From 182902d16cef899a172c9e3a644043e30e2317cd Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Thu, 23 Jan 2025 17:18:30 -0500 Subject: [PATCH 1/6] WIP: deltaprime symmetrization outline --- stride/sym.f | 156 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 stride/sym.f diff --git a/stride/sym.f b/stride/sym.f new file mode 100644 index 00000000..42483e63 --- /dev/null +++ b/stride/sym.f @@ -0,0 +1,156 @@ +c----------------------------------------------------------------------- +c file stride_netcdf.f +c writes stride.out information to a netcdf file +c----------------------------------------------------------------------- +c code organization. +c----------------------------------------------------------------------- +c 0. stride_netcdf_mod +c 1. check +c 2. stride_netcdf_out +c----------------------------------------------------------------------- +c subprogram 0. stride_netcdf_mod +c module declarations. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + MODULE sym_mod + + USE dcon_mod + + IMPLICIT NONE + + CONTAINS +c----------------------------------------------------------------------- +c subprogram 1. symmetrize. +c Generate and symmetrize Delta matrices. +c----------------------------------------------------------------------- + SUBROUTINE symmetrize(msing,dp) + + COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp + + INTEGER, INTENT(IN) :: msing, cmpsi + COMPLEX(r8), INTENT(IN) :: dp + sq + + + INTEGER :: ising,jsing + LOGICAL :: is_hermitian + REAL(r8) :: tolerance + COMPLEX(r8), DIMENSION(msing,msing) :: A_prime,B_prime, + $ Gamma_prime,Delta_prime + COMPLEX(r8), DIMENSION(msing,msing) :: A_prime_tmp,B_prime_tmp, + $ Gamma_prime_tmp,Delta_prime_tmp + COMPLEX(r8), DIMENSION(msing,msing) :: A_prime_sym,B_prime_sym, + $ Gamma_prime_sym,Delta_prime_sym + + ! Construct PEST3 matrices + ! construct PEST3 matching data (keep synced with RDCON!) + A_prime=0.0 + B_prime=0.0 + Gamma_prime=0.0 + Delta_prime=0.0 + DO ising=1,msing + DO jsing=1,msing + A_prime(ising,jsing)=dp(2*ising,2*jsing) + $ +dp(2*ising,2*jsing-1) + $ +dp(2*ising-1,2*jsing) + $ +dp(2*ising-1,2*jsing-1) + B_prime(ising,jsing)=dp(2*ising,2*jsing) + $ -dp(2*ising,2*jsing-1) + $ +dp(2*ising-1,2*jsing) + $ -dp(2*ising-1,2*jsing-1) + Gamma_prime(ising,jsing)=dp(2*ising,2*jsing) + $ +dp(2*ising,2*jsing-1) + $ -dp(2*ising-1,2*jsing) + $ -dp(2*ising-1,2*jsing-1) + Delta_prime(ising,jsing)=dp(2*ising,2*jsing) + $ -dp(2*ising,2*jsing-1) + $ -dp(2*ising-1,2*jsing) + $ +dp(2*ising-1,2*jsing-1) + ENDDO + ENDDO + + ! STILL NEED G(ipsi), psi_a + I = 0 + DO ipsi=1,cmpsi + CALL spline_eval(sq,psi(ipsi),0) + I = I + 2.0*(sq%f(4) / G(ipsi)) !DB !!! + + ENDDO + + DO ising=1,msing + resnum=NINT(singtype(ising)%q*nn)-mlow+1 + respsi=singtype(ising)%psifac + CALL spline_eval(sq,respsi,1) + + DO itheta=0,mthsurf + CALL bicube_eval(rzphi,respsi,theta(itheta),1) + jac=rzphi%f(4) + w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac + w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) + delpsi(itheta)=SQRT(w(1,1)**2+w(1,2)**2) + + J(ising)=J(ising)+(1.0/delpsi(itheta)**2)*jac!*(1.0/twopi) !DB !!! + ENDDO + + rho(ising) = (J(ising) * G(respsi))/sq%f(4) !DB + shear(ising)=mfac(resnum)*sq%f1(4)/sq%f(4)**2 + + L(ising) = I * (J(ising) * (( resm * G(ising) )/sq%f(4))**2+ + $ (nn * psi_a)**2 ) !DB + + ! CALL mercier_mod + DI(ising) = locstab%fs(respsi,1)/sq%xs(respsi) !DB + + nu_L(ising) = 0.5 - SQRT(-DI(ising)) !DB + nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB + + f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising) - + $ nu_L(ising)) / L(ising))**0.5)*shear(ising)*resm !DB + f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising) - + $ nu_L(ising)) / L(ising))**0.5)*shear(ising)*resm !DB + + A_prime_tmp = MATMUL(A_prime,f_L) !DB + B_prime_tmp = MATMUL(B_prime,f_L) !DB + Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB + Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB + + A_prime_sym = MATMUL(RESHAPE([f_S], [1,n]),A_prime_tmp) !DB + B_prime_sym = MATMUL(RESHAPE([f_S], [1,n]),B_prime_tmp) !DB + Gamma_prime_sym = MATMUL(RESHAPE([f_S], + $ [1,n]),Gamma_prime_tmp) !DB + Delta_prime_sym = MATMUL(RESHAPE([f_S], + $ [1,n]),Delta_prime_tmp) !DB + ENDDO + + tolerance = 1.0e-2 + + is_hermitian = all(abs(A_prime_sym - + $ transpose(conjg(A_prime_sym))) < tolerance) + IF (is_hermitian) then + PRINT(*,*), "A_prime is Hermitian" + ELSE + PRINT(*,*), "A_prime is not Hermitian" + END IF + is_hermitian = all(abs(B_prime_sym - + $ transpose(conjg(B_prime_sym))) < tolerance) + IF (is_hermitian) then + PRINT(*,*), "B_prime is Hermitian" + ELSE + PRINT(*,*), "B_prime is not Hermitian" + END IF + is_hermitian = all(abs(Gamma_prime_sym - + $ transpose(conjg(Gamma_prime_sym))) < tolerance) + IF (is_hermitian) then + PRINT(*,*), "Gamma_prime_sym is Hermitian" + ELSE + PRINT(*,*), "Gamma_prime_sym is not Hermitian" + END IF + is_hermitian = all(abs(Delta_prime_sym - + $ transpose(conjg(Delta_prime_sym))) < tolerance) + IF (is_hermitian) then + PRINT(*,*), "Delta_prime_sym is Hermitian" + ELSE + PRINT(*,*), "Delta_prime_sym is not Hermitian" + END IF \ No newline at end of file From f2713f5f273796f33fac94989f4dfbb5e81a1cb4 Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Tue, 28 Jan 2025 18:45:03 -0500 Subject: [PATCH 2/6] WIP: STRIDE sym dependency updates --- stride/free.f | 2 + stride/makefile | 4 +- stride/sym.f | 101 +++++++++++++++++++++++++++++------------------- 3 files changed, 67 insertions(+), 40 deletions(-) diff --git a/stride/free.f b/stride/free.f index 23ffb84b..57f8c92e 100755 --- a/stride/free.f +++ b/stride/free.f @@ -21,6 +21,7 @@ c----------------------------------------------------------------------- MODULE free_mod USE ode_mod + USE sym_mod USE stride_netcdf_mod IMPLICIT NONE @@ -292,6 +293,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) c optionally write netcdf file. c----------------------------------------------------------------------- IF(present(op_netcdf_out))THEN + IF(op_netcdf_out) CALL symmetrize(delta_prime_mat) IF(op_netcdf_out) CALL stride_netcdf_out(wp,wv,wt,ep,ev,et, $ delta_prime_mat) ENDIF diff --git a/stride/makefile b/stride/makefile index 22a8b775..32e63187 100755 --- a/stride/makefile +++ b/stride/makefile @@ -29,6 +29,7 @@ OBJS = \ sing.o \ zvode1.o \ ode.o \ + sym.o \ stride_netcdf.o \ free.o \ riccati.o \ @@ -85,6 +86,7 @@ sing.o: fourfit.o zvode1.o: ../zvode/zvode.o ../zlange/zlange.o riccati.o: sing.o ode.o: sing.o debug.o zvode1.o riccati.o +sym.o: dcon_mod.o stride_netcdf.o: dcon_mod.o -free.o: ode.o stride_netcdf.o +free.o: ode.o stride_netcdf.o sym.o stride.o: bal.o free.o mercier.o diff --git a/stride/sym.f b/stride/sym.f index 42483e63..43e8b112 100644 --- a/stride/sym.f +++ b/stride/sym.f @@ -25,26 +25,36 @@ MODULE sym_mod c subprogram 1. symmetrize. c Generate and symmetrize Delta matrices. c----------------------------------------------------------------------- - SUBROUTINE symmetrize(msing,dp) + SUBROUTINE symmetrize(dp) COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp - INTEGER, INTENT(IN) :: msing, cmpsi - COMPLEX(r8), INTENT(IN) :: dp - sq - - - INTEGER :: ising,jsing + INTEGER :: m,ising,jsing,itheta,mthsurf,npsi,ipsi + INTEGER, DIMENSION(mpert) :: mvec + REAL(r8) :: resnum,resm_sing, I + INTEGER, DIMENSION(msing) :: resm LOGICAL :: is_hermitian REAL(r8) :: tolerance + REAL(r8), DIMENSION(0:mtheta) :: r,z,theta + + REAL(r8),DIMENSION(msing) :: L,f_L,f_S,nu_L,nu_S,DI,J,rho,shear + COMPLEX(r8), DIMENSION(msing,msing) :: A_prime,B_prime, $ Gamma_prime,Delta_prime - COMPLEX(r8), DIMENSION(msing,msing) :: A_prime_tmp,B_prime_tmp, + COMPLEX(r8), DIMENSION(msing) :: A_prime_tmp,B_prime_tmp, $ Gamma_prime_tmp,Delta_prime_tmp COMPLEX(r8), DIMENSION(msing,msing) :: A_prime_sym,B_prime_sym, $ Gamma_prime_sym,Delta_prime_sym - ! Construct PEST3 matrices + REAL(r8) :: delpsi,jac,psifac,q,q1,rfac + REAL(r8) :: respsi,psi_a + REAL(r8), DIMENSION(2,2) :: w + TYPE(spline_type) :: spl + + npsi = sq%mx+1 + + psi_a=0.1 !!!!!!!!!!!!!! + ! construct PEST3 matching data (keep synced with RDCON!) A_prime=0.0 B_prime=0.0 @@ -71,17 +81,24 @@ SUBROUTINE symmetrize(msing,dp) ENDDO ENDDO - ! STILL NEED G(ipsi), psi_a - I = 0 - DO ipsi=1,cmpsi - CALL spline_eval(sq,psi(ipsi),0) - I = I + 2.0*(sq%f(4) / G(ipsi)) !DB !!! + mvec=(/(m,m=mlow,mhigh)/) + DO ising=1,msing + respsi=sing(ising)%psifac + resnum=NINT(sing(ising)%q*nn)-mlow+1 + resm_sing=mvec(resnum) + resm(ising)=resm_sing + ENDDO + ! STILL NEED psi_a + I = 0 + DO ipsi=1,npsi + !CALL spline_eval(sq,psi(ipsi),0) + I = I + 2.0*(sq%f(4) / (sq%fs(ipsi,1)/twopi)) !DB this is F! ENDDO DO ising=1,msing - resnum=NINT(singtype(ising)%q*nn)-mlow+1 - respsi=singtype(ising)%psifac + !resnum=NINT(sing(ising)%q*nn)-mlow+1 + respsi=sing(ising)%psifac CALL spline_eval(sq,respsi,1) DO itheta=0,mthsurf @@ -89,16 +106,16 @@ SUBROUTINE symmetrize(msing,dp) jac=rzphi%f(4) w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - delpsi(itheta)=SQRT(w(1,1)**2+w(1,2)**2) + delpsi=SQRT(w(1,1)**2+w(1,2)**2) - J(ising)=J(ising)+(1.0/delpsi(itheta)**2)*jac!*(1.0/twopi) !DB !!! + J(ising)=J(ising)+(1.0/delpsi**2)*jac!*(1.0/twopi) !DB !!! ENDDO - rho(ising) = (J(ising) * G(respsi))/sq%f(4) !DB - shear(ising)=mfac(resnum)*sq%f1(4)/sq%f(4)**2 + rho(ising) = (J(ising) * (sq%fs(respsi,1)/twopi))/sq%f(4) !DB + shear(ising)=sing(ising)%q1 - L(ising) = I * (J(ising) * (( resm * G(ising) )/sq%f(4))**2+ - $ (nn * psi_a)**2 ) !DB + L(ising)=I*(J(ising)*((resm(ising)*(sq%fs(respsi,1)/twopi) + $ )/sq%f(4))**2+(nn * psi_a)**2 ) !DB ! CALL mercier_mod DI(ising) = locstab%fs(respsi,1)/sq%xs(respsi) !DB @@ -107,21 +124,23 @@ SUBROUTINE symmetrize(msing,dp) nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising) - - $ nu_L(ising)) / L(ising))**0.5)*shear(ising)*resm !DB + $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising) - - $ nu_L(ising)) / L(ising))**0.5)*shear(ising)*resm !DB + $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB A_prime_tmp = MATMUL(A_prime,f_L) !DB B_prime_tmp = MATMUL(B_prime,f_L) !DB Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB - A_prime_sym = MATMUL(RESHAPE([f_S], [1,n]),A_prime_tmp) !DB - B_prime_sym = MATMUL(RESHAPE([f_S], [1,n]),B_prime_tmp) !DB - Gamma_prime_sym = MATMUL(RESHAPE([f_S], - $ [1,n]),Gamma_prime_tmp) !DB - Delta_prime_sym = MATMUL(RESHAPE([f_S], - $ [1,n]),Delta_prime_tmp) !DB + A_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + $ RESHAPE([A_prime_tmp], [1,msing])) !DB + B_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + $ RESHAPE([B_prime_tmp], [1,msing])) !DB + Gamma_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + $ RESHAPE([Gamma_prime_tmp], [1,msing])) !DB + Delta_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + $ RESHAPE([Delta_prime_tmp], [1,msing])) !DB ENDDO tolerance = 1.0e-2 @@ -129,28 +148,32 @@ SUBROUTINE symmetrize(msing,dp) is_hermitian = all(abs(A_prime_sym - $ transpose(conjg(A_prime_sym))) < tolerance) IF (is_hermitian) then - PRINT(*,*), "A_prime is Hermitian" + WRITE(*,*), "A_prime is Hermitian" ELSE - PRINT(*,*), "A_prime is not Hermitian" + WRITE(*,*), "A_prime is not Hermitian" END IF is_hermitian = all(abs(B_prime_sym - $ transpose(conjg(B_prime_sym))) < tolerance) IF (is_hermitian) then - PRINT(*,*), "B_prime is Hermitian" + WRITE(*,*), "B_prime is Hermitian" ELSE - PRINT(*,*), "B_prime is not Hermitian" + WRITE(*,*), "B_prime is not Hermitian" END IF is_hermitian = all(abs(Gamma_prime_sym - $ transpose(conjg(Gamma_prime_sym))) < tolerance) IF (is_hermitian) then - PRINT(*,*), "Gamma_prime_sym is Hermitian" + WRITE(*,*), "Gamma_prime_sym is Hermitian" ELSE - PRINT(*,*), "Gamma_prime_sym is not Hermitian" + WRITE(*,*), "Gamma_prime_sym is not Hermitian" END IF is_hermitian = all(abs(Delta_prime_sym - $ transpose(conjg(Delta_prime_sym))) < tolerance) IF (is_hermitian) then - PRINT(*,*), "Delta_prime_sym is Hermitian" + WRITE(*,*), "Delta_prime_sym is Hermitian" ELSE - PRINT(*,*), "Delta_prime_sym is not Hermitian" - END IF \ No newline at end of file + WRITE(*,*), "Delta_prime_sym is not Hermitian" + END IF + + RETURN + END SUBROUTINE symmetrize + END MODULE sym_mod \ No newline at end of file From dce3e95a1fa60f7de446591c1af305fef0aee6a4 Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Fri, 31 Jan 2025 14:56:42 -0500 Subject: [PATCH 3/6] WIP: broken, but integral progress --- stride/sym.f | 107 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 81 insertions(+), 26 deletions(-) diff --git a/stride/sym.f b/stride/sym.f index 43e8b112..de67b3f9 100644 --- a/stride/sym.f +++ b/stride/sym.f @@ -29,14 +29,14 @@ SUBROUTINE symmetrize(dp) COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp - INTEGER :: m,ising,jsing,itheta,mthsurf,npsi,ipsi + INTEGER :: m,ising,jsing,itheta,mthsurf,ipsi,jpsi,npsi,idx INTEGER, DIMENSION(mpert) :: mvec - REAL(r8) :: resnum,resm_sing, I + REAL(r8) :: resnum,resm_sing,min_diff,dtheta INTEGER, DIMENSION(msing) :: resm LOGICAL :: is_hermitian - REAL(r8) :: tolerance + REAL(r8) :: tolerance, dx REAL(r8), DIMENSION(0:mtheta) :: r,z,theta - + REAL(r8), DIMENSION(sq%mx+1) :: I_psin REAL(r8),DIMENSION(msing) :: L,f_L,f_S,nu_L,nu_S,DI,J,rho,shear COMPLEX(r8), DIMENSION(msing,msing) :: A_prime,B_prime, @@ -51,10 +51,15 @@ SUBROUTINE symmetrize(dp) REAL(r8), DIMENSION(2,2) :: w TYPE(spline_type) :: spl - npsi = sq%mx+1 + npsi = sq%mx psi_a=0.1 !!!!!!!!!!!!!! + WRITE(*,*)"sq%mx=",sq%mx + WRITE(*,*)"npsi=",npsi + !WRITE(*,*)"sq%xs=",sq%xs + + ! construct PEST3 matching data (keep synced with RDCON!) A_prime=0.0 B_prime=0.0 @@ -89,36 +94,86 @@ SUBROUTINE symmetrize(dp) resm(ising)=resm_sing ENDDO + !WRITE(*,*)"respsi=",respsi + WRITE(*,*)"resnum=",resnum + WRITE(*,*)"resm=",resm + ! STILL NEED psi_a - I = 0 - DO ipsi=1,npsi - !CALL spline_eval(sq,psi(ipsi),0) - I = I + 2.0*(sq%f(4) / (sq%fs(ipsi,1)/twopi)) !DB this is F! - ENDDO + !I = 0 + !DO ipsi=1,npsi + ! !CALL spline_eval(sq,psi(ipsi),0) + ! I = I + 2.0*(sq%f(4) / (sq%fs(ipsi,1)/twopi)) !DB this is F! + !ENDDO - DO ising=1,msing - !resnum=NINT(sing(ising)%q*nn)-mlow+1 - respsi=sing(ising)%psifac - CALL spline_eval(sq,respsi,1) + I_psin(1) = 0.0 ! integral at x=0 is 0 + + !WRITE(*,*)"sq%fs(:,4)=",sq%fs(:,4) + !WRITE(*,*)"SIZE(sq%fs(:,4))=",SIZE(sq%fs(:,4)) - DO itheta=0,mthsurf - CALL bicube_eval(rzphi,respsi,theta(itheta),1) - jac=rzphi%f(4) - w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac - w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - delpsi=SQRT(w(1,1)**2+w(1,2)**2) + !WRITE(*,*)"sq%fs(npsi,4)=",sq%fs(npsi,4) - J(ising)=J(ising)+(1.0/delpsi**2)*jac!*(1.0/twopi) !DB !!! - ENDDO + !WRITE(*,*)"sq%fs(5,4)=",sq%fs(5,4) + + DO ipsi = 2, npsi + I_psin(ipsi) = 0.0 + DO jpsi = 1, ipsi-1 + dx = sq%xs(jpsi+1) - sq%xs(jpsi) + I_psin(ipsi) = I_psin(ipsi) + + $ (2.0*sq%fs(jpsi,4)/(sq%fs(jpsi,1)/twopi))*dx ! check q + END DO + END DO - rho(ising) = (J(ising) * (sq%fs(respsi,1)/twopi))/sq%f(4) !DB + !WRITE(*,*)"I_psin=",I_psin + + DO ising=1,msing + !resnum=NINT(sing(ising)%q*nn)-mlow+1 + respsi=sing(ising)%psifac + WRITE(*,*)"respsi=",respsi + + !!! FIND CORRECT I_PSIN INDEX + min_diff = abs(sq%xs(1) - respsi) + idx = 1 + + do ipsi = 2, npsi + if (abs(sq%xs(ipsi) - respsi) < min_diff) then + min_diff = abs(sq%xs(ipsi) - respsi) + idx = ipsi + end if + end do + + WRITE(*,*)"idx=",idx + + !CALL spline_eval(sq,respsi,1) + dtheta = twopi / real(mthsurf) + + !DO itheta=0,mthsurf + ! CALL bicube_eval(rzphi,respsi,theta(itheta),1) + ! jac=rzphi%f(4) + ! w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac + ! w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) + ! delpsi=SQRT(w(1,1)**2+w(1,2)**2) + ! WRITE(*,*)"delpsi=",delpsi + + ! J(ising)=J(ising)+(1.0/delpsi**2)*0.5*dtheta!jac!*(1.0/twopi) !DB !!! + !ENDDO + + WRITE(*,*)"J=",J + + rho(ising) = (J(ising) * (sq%fs(idx,1)/twopi))/sq%f(4) !DB shear(ising)=sing(ising)%q1 + WRITE(*,*)"rho=",rho + WRITE(*,*)"sq%f(4)=",sq%f(4) + + L(ising)=I_psin(idx)*(J(ising)*((resm(ising)* + $ (sq%fs(idx,1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB - L(ising)=I*(J(ising)*((resm(ising)*(sq%fs(respsi,1)/twopi) - $ )/sq%f(4))**2+(nn * psi_a)**2 ) !DB + WRITE(*,*)"L=",L + WRITE(*,*)"sq%fs(respsi,1)=",sq%fs(respsi,1) + WRITE(*,*)"sq%fs(idx,1)=",sq%fs(idx,1) ! CALL mercier_mod - DI(ising) = locstab%fs(respsi,1)/sq%xs(respsi) !DB + DI(ising) = locstab%fs(idx,1)/sq%xs(idx) !DB + WRITE(*,*)"DI=",DI nu_L(ising) = 0.5 - SQRT(-DI(ising)) !DB nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB From 4568a5c06b7b0655827dd11d95ae3d081cbb7fcd Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Tue, 4 Feb 2025 14:54:56 -0500 Subject: [PATCH 4/6] WIP: looks roughly correct but yields wrong answers --- stride/sym.f | 107 +++++++++++++++++++++++++++++---------------------- 1 file changed, 62 insertions(+), 45 deletions(-) diff --git a/stride/sym.f b/stride/sym.f index de67b3f9..fffbc952 100644 --- a/stride/sym.f +++ b/stride/sym.f @@ -34,8 +34,8 @@ SUBROUTINE symmetrize(dp) REAL(r8) :: resnum,resm_sing,min_diff,dtheta INTEGER, DIMENSION(msing) :: resm LOGICAL :: is_hermitian - REAL(r8) :: tolerance, dx - REAL(r8), DIMENSION(0:mtheta) :: r,z,theta + REAL(r8) :: tolerance,dx,angle,chi1 + REAL(r8), DIMENSION(:), ALLOCATABLE :: r,z,theta REAL(r8), DIMENSION(sq%mx+1) :: I_psin REAL(r8),DIMENSION(msing) :: L,f_L,f_S,nu_L,nu_S,DI,J,rho,shear @@ -51,14 +51,19 @@ SUBROUTINE symmetrize(dp) REAL(r8), DIMENSION(2,2) :: w TYPE(spline_type) :: spl + WRITE(*,*)"mthsurf=",mthsurf + + mthsurf=40*mthsurf0*MAX(ABS(mlow),ABS(mhigh)) + ALLOCATE(r(0:mthsurf),z(0:mthsurf),theta(0:mthsurf)) + npsi = sq%mx - psi_a=0.1 !!!!!!!!!!!!!! + !chi1=twopi*psio + psi_a=psio + WRITE(*,*)"psi_a=",psi_a WRITE(*,*)"sq%mx=",sq%mx WRITE(*,*)"npsi=",npsi - !WRITE(*,*)"sq%xs=",sq%xs - ! construct PEST3 matching data (keep synced with RDCON!) A_prime=0.0 @@ -94,26 +99,11 @@ SUBROUTINE symmetrize(dp) resm(ising)=resm_sing ENDDO - !WRITE(*,*)"respsi=",respsi WRITE(*,*)"resnum=",resnum WRITE(*,*)"resm=",resm - ! STILL NEED psi_a - !I = 0 - !DO ipsi=1,npsi - ! !CALL spline_eval(sq,psi(ipsi),0) - ! I = I + 2.0*(sq%f(4) / (sq%fs(ipsi,1)/twopi)) !DB this is F! - !ENDDO - + ! Calculate I(psi_N) integral I_psin(1) = 0.0 ! integral at x=0 is 0 - - !WRITE(*,*)"sq%fs(:,4)=",sq%fs(:,4) - !WRITE(*,*)"SIZE(sq%fs(:,4))=",SIZE(sq%fs(:,4)) - - !WRITE(*,*)"sq%fs(npsi,4)=",sq%fs(npsi,4) - - !WRITE(*,*)"sq%fs(5,4)=",sq%fs(5,4) - DO ipsi = 2, npsi I_psin(ipsi) = 0.0 DO jpsi = 1, ipsi-1 @@ -123,14 +113,19 @@ SUBROUTINE symmetrize(dp) END DO END DO - !WRITE(*,*)"I_psin=",I_psin + ! Prepare theta quantities for J(psi_N) integral + CALL spline_alloc(spl,mthsurf,4) + theta=(/(itheta,itheta=0,mthsurf)/)/REAL(mthsurf,r8) + spl%xs=theta + psifac=psilim + ! Loop across rational surfaces to evaluate remaining quantities DO ising=1,msing !resnum=NINT(sing(ising)%q*nn)-mlow+1 respsi=sing(ising)%psifac WRITE(*,*)"respsi=",respsi - !!! FIND CORRECT I_PSIN INDEX + ! Find correct I(psiN) index min_diff = abs(sq%xs(1) - respsi) idx = 1 @@ -143,51 +138,65 @@ SUBROUTINE symmetrize(dp) WRITE(*,*)"idx=",idx - !CALL spline_eval(sq,respsi,1) + ! Evaluate splines on rational surface + CALL spline_eval(sq,respsi,1) + CALL spline_eval(locstab,respsi,1) + + WRITE(*,*)"SIZE(theta)=",SIZE(theta) + WRITE(*,*)"mthsurf=",mthsurf + + ! Calculate J(psi_N) contour integral dtheta = twopi / real(mthsurf) + DO itheta=0,mthsurf + CALL bicube_eval(rzphi,psilim,theta(itheta),1) + rfac=SQRT(rzphi%f(1)) + angle=twopi*(theta(itheta)+rzphi%f(2)) + r(itheta)=ro+rfac*COS(angle) + z(itheta)=zo+rfac*SIN(angle) + jac=rzphi%f(4) - !DO itheta=0,mthsurf - ! CALL bicube_eval(rzphi,respsi,theta(itheta),1) - ! jac=rzphi%f(4) - ! w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac - ! w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - ! delpsi=SQRT(w(1,1)**2+w(1,2)**2) - ! WRITE(*,*)"delpsi=",delpsi + w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac + w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - ! J(ising)=J(ising)+(1.0/delpsi**2)*0.5*dtheta!jac!*(1.0/twopi) !DB !!! - !ENDDO + delpsi=SQRT(w(1,1)**2+w(1,2)**2) + + J(ising)=J(ising)+ + $ (1.0/(delpsi**2))*(dtheta/twopi)!/jac !DB !! + ENDDO WRITE(*,*)"J=",J - rho(ising) = (J(ising) * (sq%fs(idx,1)/twopi))/sq%f(4) !DB - shear(ising)=sing(ising)%q1 + rho(ising) = (J(ising) * (sq%f(1)/twopi))/sq%f(4) !DB + shear(ising) = sing(ising)%q1 WRITE(*,*)"rho=",rho + WRITE(*,*)"shear(ising)=",shear(ising) + WRITE(*,*)"sq%f(4)=",sq%f(4) + ! Combine integrated quantities into L(psi_N) L(ising)=I_psin(idx)*(J(ising)*((resm(ising)* - $ (sq%fs(idx,1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB + $ (sq%f(1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB WRITE(*,*)"L=",L - WRITE(*,*)"sq%fs(respsi,1)=",sq%fs(respsi,1) - WRITE(*,*)"sq%fs(idx,1)=",sq%fs(idx,1) - ! CALL mercier_mod - DI(ising) = locstab%fs(idx,1)/sq%xs(idx) !DB + DI(ising) = locstab%f(1)/respsi !DB WRITE(*,*)"DI=",DI nu_L(ising) = 0.5 - SQRT(-DI(ising)) !DB nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB - f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising) - + f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising)- $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB - f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising) - + f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising)- $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB + ! First component of vector*matrix*vector multiply A_prime_tmp = MATMUL(A_prime,f_L) !DB B_prime_tmp = MATMUL(B_prime,f_L) !DB Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB + ! Second component of vector*matrix*vector multiply A_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), $ RESHAPE([A_prime_tmp], [1,msing])) !DB B_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), @@ -198,7 +207,15 @@ SUBROUTINE symmetrize(dp) $ RESHAPE([Delta_prime_tmp], [1,msing])) !DB ENDDO - tolerance = 1.0e-2 + WRITE(*,*)"Delta_prime=",Delta_prime + WRITE(*,*)"Delta_prime_sym=",Delta_prime_sym + WRITE(*,*)"Delta_diff=",abs(Delta_prime_sym - + $ transpose(conjg(Delta_prime_sym))) + !WRITE(*,*)"B_prime=",B_prime + !WRITE(*,*)"Gamma_prime_sym=",Gamma_prime_sym + WRITE(*,*)"B_Gamma_diff=",(B_prime_sym - + $ transpose(conjg(Gamma_prime_sym))) + tolerance = 1.0e-01 is_hermitian = all(abs(A_prime_sym - $ transpose(conjg(A_prime_sym))) < tolerance) @@ -207,14 +224,14 @@ SUBROUTINE symmetrize(dp) ELSE WRITE(*,*), "A_prime is not Hermitian" END IF - is_hermitian = all(abs(B_prime_sym - + is_hermitian = all(abs(Gamma_prime_sym - $ transpose(conjg(B_prime_sym))) < tolerance) IF (is_hermitian) then WRITE(*,*), "B_prime is Hermitian" ELSE WRITE(*,*), "B_prime is not Hermitian" END IF - is_hermitian = all(abs(Gamma_prime_sym - + is_hermitian = all(abs(B_prime_sym - $ transpose(conjg(Gamma_prime_sym))) < tolerance) IF (is_hermitian) then WRITE(*,*), "Gamma_prime_sym is Hermitian" From 897cda8236438ca8478929706b644da8d2e25e3c Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Mon, 17 Feb 2025 15:53:48 -0500 Subject: [PATCH 5/6] WIP: check original delta matrices for symmetry --- stride/sym.f | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/stride/sym.f b/stride/sym.f index fffbc952..4f612472 100644 --- a/stride/sym.f +++ b/stride/sym.f @@ -207,6 +207,13 @@ SUBROUTINE symmetrize(dp) $ RESHAPE([Delta_prime_tmp], [1,msing])) !DB ENDDO + WRITE(*,*)"Original Delta diff=",(Delta_prime - + $ transpose(conjg(Delta_prime))) + WRITE(*,*)"Original A diff=",(A_prime - + $ transpose(conjg(A_prime))) + WRITE(*,*)"Original B_Gamma diff=",(B_prime - + $ transpose(conjg(Gamma_prime))) + WRITE(*,*)"Delta_prime=",Delta_prime WRITE(*,*)"Delta_prime_sym=",Delta_prime_sym WRITE(*,*)"Delta_diff=",abs(Delta_prime_sym - From 6b46ca4aed5ff34fc042012d37f93df1c64d030d Mon Sep 17 00:00:00 2001 From: Daniel Burgess Date: Thu, 26 Jun 2025 17:26:19 -0400 Subject: [PATCH 6/6] new splining and correct integral(J) --- stride/stride.F | 2 +- stride/sym.f | 311 ++++++++++++++++++++++++++++++------------------ 2 files changed, 195 insertions(+), 118 deletions(-) diff --git a/stride/stride.F b/stride/stride.F index 3104f397..77dad123 100644 --- a/stride/stride.F +++ b/stride/stride.F @@ -157,7 +157,7 @@ PROGRAM stride PRINT *, " > Forcing reform_eq_with_psilim=t" ENDIF IF(psilim /= psihigh .OR. psilow /= sq%xs(0))THEN - psilow_tmp = psilow ! if we feed psilow directly, it get's overwritten by namelist read + psilow_tmp = psilow ! if we feed psilow directly, it gets overwritten by namelist read psilim_tmp = psilim CALL equil_read(out_unit, psilim_tmp, psilow_tmp) CALL equil_out_global diff --git a/stride/sym.f b/stride/sym.f index 4f612472..3b4b954a 100644 --- a/stride/sym.f +++ b/stride/sym.f @@ -16,7 +16,7 @@ c----------------------------------------------------------------------- MODULE sym_mod - USE dcon_mod + USE stride_dcon_mod IMPLICIT NONE @@ -29,14 +29,17 @@ SUBROUTINE symmetrize(dp) COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp - INTEGER :: m,ising,jsing,itheta,mthsurf,ipsi,jpsi,npsi,idx + INTEGER :: m,ising,jsing,itheta,mthsurf,ipsi,jpsi,npsi,idx, + $ si,sj + REAL(r8) :: bsq,chi1,dpsisq,eta,jac,p1,psifac,q,q1,r, + $ rfac,theta,twopif,v1,v2,v21,v22,v23,v33 INTEGER, DIMENSION(mpert) :: mvec REAL(r8) :: resnum,resm_sing,min_diff,dtheta INTEGER, DIMENSION(msing) :: resm LOGICAL :: is_hermitian - REAL(r8) :: tolerance,dx,angle,chi1 - REAL(r8), DIMENSION(:), ALLOCATABLE :: r,z,theta - REAL(r8), DIMENSION(sq%mx+1) :: I_psin + REAL(r8) :: tolerance,dx,angle + !REAL(r8), DIMENSION(:), ALLOCATABLE :: r,z,theta + REAL(r8), DIMENSION(sq%mx+1) :: I_psin,ln_q REAL(r8),DIMENSION(msing) :: L,f_L,f_S,nu_L,nu_S,DI,J,rho,shear COMPLEX(r8), DIMENSION(msing,msing) :: A_prime,B_prime, @@ -44,22 +47,25 @@ SUBROUTINE symmetrize(dp) COMPLEX(r8), DIMENSION(msing) :: A_prime_tmp,B_prime_tmp, $ Gamma_prime_tmp,Delta_prime_tmp COMPLEX(r8), DIMENSION(msing,msing) :: A_prime_sym,B_prime_sym, - $ Gamma_prime_sym,Delta_prime_sym + $ Gamma_prime_sym,Delta_prime_sym,tmp_arr - REAL(r8) :: delpsi,jac,psifac,q,q1,rfac + REAL(r8) :: delpsi REAL(r8) :: respsi,psi_a REAL(r8), DIMENSION(2,2) :: w - TYPE(spline_type) :: spl + TYPE(spline_type) :: spl,psi_t,I_spl,shr_spl,J_spl - WRITE(*,*)"mthsurf=",mthsurf + REAL(r8), DIMENSION(:), POINTER :: avg + TYPE(spline_type), TARGET :: fspl - mthsurf=40*mthsurf0*MAX(ABS(mlow),ABS(mhigh)) - ALLOCATE(r(0:mthsurf),z(0:mthsurf),theta(0:mthsurf)) + !WRITE(*,*)"mthsurf=",mthsurf + + !mthsurf=40*mthsurf0*MAX(ABS(mlow),ABS(mhigh)) + !ALLOCATE(r(0:mthsurf),z(0:mthsurf),theta(0:mthsurf)) npsi = sq%mx !chi1=twopi*psio - psi_a=psio + psi_a=twopi*psio WRITE(*,*)"psi_a=",psi_a WRITE(*,*)"sq%mx=",sq%mx @@ -102,149 +108,220 @@ SUBROUTINE symmetrize(dp) WRITE(*,*)"resnum=",resnum WRITE(*,*)"resm=",resm + ! Prepare toroidal flux spline + CALL spline_alloc(psi_t,SIZE(sq%fsi(:, 4)),1) + psi_t%xs=sq%xs(:) + psi_t%fs(:,1)=sq%fsi(:,4)*twopi*psio ! Un-normalize toroidal flux + CALL spline_fit(psi_t,"extrap") + + ! Prepare shear spline + CALL spline_alloc(shr_spl,SIZE(sq%xs(:)),1) + shr_spl%xs=sq%xs(:) + ln_q=LOG(sq%fs(:,4)) + shr_spl%fs(:,1)=ln_q! log(q) + CALL spline_fit(shr_spl,"extrap") + ! Calculate I(psi_N) integral I_psin(1) = 0.0 ! integral at x=0 is 0 DO ipsi = 2, npsi I_psin(ipsi) = 0.0 DO jpsi = 1, ipsi-1 dx = sq%xs(jpsi+1) - sq%xs(jpsi) + CALL spline_eval(psi_t,sq%xs(jpsi),1) I_psin(ipsi) = I_psin(ipsi) + - $ (2.0*sq%fs(jpsi,4)/(sq%fs(jpsi,1)/twopi))*dx ! check q - END DO - END DO + $ (2.0*sq%fs(jpsi,4)/(psi_t%f(1)))*dx + ENDDO + ENDDO + + ! Prepare I(psiN) spline + CALL spline_alloc(I_spl,SIZE(sq%xs(:)),1) + I_spl%xs=sq%xs(:) + I_spl%fs(:,1)=I_psin + CALL spline_fit(I_spl,"extrap") + + ! Compute J integral + CALL spline_alloc(fspl,mtheta,2) + fspl%xs=rzphi%ys + CALL spline_alloc(J_spl,mpsi,1) + J_spl%xs=sq%xs(:) + WRITE(*,*)'(mtheta)=',mtheta + + DO ipsi=0,mpsi + psifac=sq%xs(ipsi) + twopif=sq%fs(ipsi,1) + p1=sq%fs1(ipsi,2) + v1=sq%fs(ipsi,3) + v2=sq%fs1(ipsi,3) + q=sq%fs(ipsi,4) + q1=sq%fs1(ipsi,4) + chi1=twopi*psio +c----------------------------------------------------------------------- +c evaluate coordinates and jacobian. +c----------------------------------------------------------------------- + DO itheta=0,mtheta + CALL bicube_eval(rzphi,rzphi%xs(ipsi),rzphi%ys(itheta),1) + theta=rzphi%ys(itheta) + rfac=SQRT(rzphi%f(1)) + eta=twopi*(theta+rzphi%f(2)) + r=ro+rfac*COS(eta) + jac=rzphi%f(4) +c----------------------------------------------------------------------- +c evaluate other local quantities. +c----------------------------------------------------------------------- + v21=rzphi%fy(1)/(2*rfac*jac) + v22=(1+rzphi%fy(2))*twopi*rfac/jac + v23=rzphi%fy(3)*r/jac + v33=twopi*r/jac + bsq=chi1**2*(v21**2+v22**2+(v23+q*v33)**2) + dpsisq=(twopi*r)**2*(v21**2+v22**2) +c----------------------------------------------------------------------- +c evaluate integrands. +c----------------------------------------------------------------------- + fspl%fs(itheta,1)=1/dpsisq!*(chi1**2.0)) dpsisq is grad(psi_N)**2 ?? + fspl%fs(itheta,2)=bsq + fspl%fs(itheta,:)=fspl%fs(itheta,:)*(jac/v1) + ENDDO +c----------------------------------------------------------------------- +c integrate quantities with respect to theta. +c----------------------------------------------------------------------- + CALL spline_fit(fspl,"periodic") + CALL spline_int(fspl) + avg => fspl%fsi(mtheta,:) + J_spl%fs(ipsi,1)=avg(1) - ! Prepare theta quantities for J(psi_N) integral - CALL spline_alloc(spl,mthsurf,4) - theta=(/(itheta,itheta=0,mthsurf)/)/REAL(mthsurf,r8) - spl%xs=theta - psifac=psilim + ENDDO + CALL spline_dealloc(fspl) + CALL spline_fit(J_spl,"extrap") ! Loop across rational surfaces to evaluate remaining quantities DO ising=1,msing - !resnum=NINT(sing(ising)%q*nn)-mlow+1 respsi=sing(ising)%psifac WRITE(*,*)"respsi=",respsi - ! Find correct I(psiN) index - min_diff = abs(sq%xs(1) - respsi) - idx = 1 - - do ipsi = 2, npsi - if (abs(sq%xs(ipsi) - respsi) < min_diff) then - min_diff = abs(sq%xs(ipsi) - respsi) - idx = ipsi - end if - end do - - WRITE(*,*)"idx=",idx - ! Evaluate splines on rational surface CALL spline_eval(sq,respsi,1) CALL spline_eval(locstab,respsi,1) + CALL spline_eval(psi_t,respsi,1) + CALL spline_eval(I_spl,respsi,1) + CALL spline_eval(shr_spl,respsi,1) + CALL spline_eval(J_spl,respsi,1) - WRITE(*,*)"SIZE(theta)=",SIZE(theta) - WRITE(*,*)"mthsurf=",mthsurf - - ! Calculate J(psi_N) contour integral - dtheta = twopi / real(mthsurf) - DO itheta=0,mthsurf - CALL bicube_eval(rzphi,psilim,theta(itheta),1) - rfac=SQRT(rzphi%f(1)) - angle=twopi*(theta(itheta)+rzphi%f(2)) - r(itheta)=ro+rfac*COS(angle) - z(itheta)=zo+rfac*SIN(angle) - jac=rzphi%f(4) - - w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac - w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac) - - delpsi=SQRT(w(1,1)**2+w(1,2)**2) - - J(ising)=J(ising)+ - $ (1.0/(delpsi**2))*(dtheta/twopi)!/jac !DB !! - ENDDO - - WRITE(*,*)"J=",J - - rho(ising) = (J(ising) * (sq%f(1)/twopi))/sq%f(4) !DB - shear(ising) = sing(ising)%q1 - WRITE(*,*)"rho=",rho - WRITE(*,*)"shear(ising)=",shear(ising) + J(ising)=J_spl%f(1) + rho(ising) = (J(ising) * (psi_t%f(1)))/sq%f(4) !DB - WRITE(*,*)"sq%f(4)=",sq%f(4) + WRITE(*,*)"GPEC shear=",sing(ising)%q1 + WRITE(*,*)"d(ln(q))/dpsi(ising)=",shr_spl%f1(1) ! Combine integrated quantities into L(psi_N) - L(ising)=I_psin(idx)*(J(ising)*((resm(ising)* - $ (sq%f(1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB + L(ising)=I_spl%f(1)*(J(ising)*(((resm(ising)* + $ psi_t%f(1))/sq%f(4))**2.0) + (nn * psi_a)**2.0 ) !DB - WRITE(*,*)"L=",L - ! CALL mercier_mod DI(ising) = locstab%f(1)/respsi !DB - WRITE(*,*)"DI=",DI nu_L(ising) = 0.5 - SQRT(-DI(ising)) !DB nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising)- - $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB + $ nu_L(ising)) /L(ising))**0.5)*shr_spl%f1(1)*resm(ising) !DB f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising)- - $ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB + $ nu_L(ising)) /L(ising))**0.5)*shr_spl%f1(1)*resm(ising) !DB ! First component of vector*matrix*vector multiply - A_prime_tmp = MATMUL(A_prime,f_L) !DB - B_prime_tmp = MATMUL(B_prime,f_L) !DB - Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB - Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB + ! A_prime_tmp = MATMUL(A_prime,f_L) !DB + ! B_prime_tmp = MATMUL(B_prime,f_L) !DB + ! Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB + +c Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB ! Second component of vector*matrix*vector multiply - A_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), - $ RESHAPE([A_prime_tmp], [1,msing])) !DB - B_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), - $ RESHAPE([B_prime_tmp], [1,msing])) !DB - Gamma_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), - $ RESHAPE([Gamma_prime_tmp], [1,msing])) !DB - Delta_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), - $ RESHAPE([Delta_prime_tmp], [1,msing])) !DB + ! A_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + !$ RESHAPE([A_prime_tmp], [1,msing])) !DB + ! B_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + !$ RESHAPE([B_prime_tmp], [1,msing])) !DB + ! Gamma_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), + !$! RESHAPE([Gamma_prime_tmp], [1,msing])) !DB + +c Delta_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]), +c $ RESHAPE([Delta_prime_tmp], [1,msing])) !DB + ENDDO + + ! New cosine method from Richard's paper + DO si=1,msing + DO sj=1,msing + Delta_prime_sym(si,sj) = COS((si + sj)*pi)* + $ (f_S(si)/f_L(sj))*Delta_prime(si,sj) + ENDDO + ENDDO + + tmp_arr=(Delta_prime - transpose(conjg(Delta_prime))) + WRITE(*,*)"Original delta diff=" + DO si = 1, msing + DO sj = 1, msing + WRITE(*,'(F8.3)',advance='no') REAL(tmp_arr(si,sj)) + IF (sj < msing) WRITE(*,'(A)', advance='no') ' ' + ENDDO + WRITE(*,*) ! New line after each row + ENDDO + + WRITE(*,*)"Delta_prime=" + DO si = 1, msing + DO sj = 1, msing + WRITE(*,'(F8.3)', advance='no') REAL(Delta_prime(si,sj)) + IF (sj < msing) WRITE(*,'(A)', advance='no') ' ' + ENDDO + WRITE(*,*) ! New line after each row + ENDDO + WRITE(*,*)"Delta_prime_sym=" + DO si = 1, msing + DO sj = 1, msing + WRITE(*,'(F8.3)',advance='no') REAL(Delta_prime_sym(si,sj)) + IF (sj < msing) WRITE(*,'(A)', advance='no') ' ' + ENDDO + WRITE(*,*) ! New line after each row ENDDO - WRITE(*,*)"Original Delta diff=",(Delta_prime - - $ transpose(conjg(Delta_prime))) - WRITE(*,*)"Original A diff=",(A_prime - - $ transpose(conjg(A_prime))) - WRITE(*,*)"Original B_Gamma diff=",(B_prime - - $ transpose(conjg(Gamma_prime))) - - WRITE(*,*)"Delta_prime=",Delta_prime - WRITE(*,*)"Delta_prime_sym=",Delta_prime_sym - WRITE(*,*)"Delta_diff=",abs(Delta_prime_sym - - $ transpose(conjg(Delta_prime_sym))) +c WRITE(*, '(3(F5.3, 1X))') ((abs(Delta_prime_sym - +c $ transpose(conjg(Delta_prime_sym)))(i, j), j = 1, +c $ msing), i = 1, msing) + !WRITE(*,*)"B_prime=",B_prime !WRITE(*,*)"Gamma_prime_sym=",Gamma_prime_sym - WRITE(*,*)"B_Gamma_diff=",(B_prime_sym - - $ transpose(conjg(Gamma_prime_sym))) - tolerance = 1.0e-01 + !WRITE(*,*)"B_Gamma_diff=",(B_prime_sym - + !$ transpose(conjg(Gamma_prime_sym))) + tolerance = 1.0!e-01 + +c is_hermitian = all(abs(A_prime_sym - +c $ transpose(conjg(A_prime_sym))) < tolerance) +c IF (is_hermitian) then +c WRITE(*,*), "A_prime is Hermitian" +c ELSE +c WRITE(*,*), "A_prime is not Hermitian" +c END IF +c is_hermitian = all(abs(Gamma_prime_sym - +c $ transpose(conjg(B_prime_sym))) < tolerance) +c IF (is_hermitian) then +c WRITE(*,*), "B_prime is Hermitian" +c ELSE +c WRITE(*,*), "B_prime is not Hermitian" +c END IF +c is_hermitian = all(abs(B_prime_sym - +c $ transpose(conjg(Gamma_prime_sym))) < tolerance) +c IF (is_hermitian) then +c WRITE(*,*), "Gamma_prime_sym is Hermitian" +c ELSE +c WRITE(*,*), "Gamma_prime_sym is not Hermitian" +c END IF + + tmp_arr=(Delta_prime_sym - transpose(conjg(Delta_prime_sym))) + WRITE(*,*)"New delta diff=" + DO si = 1, msing + DO sj = 1, msing + WRITE(*,'(F8.3)',advance='no') REAL(tmp_arr(si,sj)) + IF (sj < msing) WRITE(*,'(A)', advance='no') ' ' + ENDDO + WRITE(*,*) ! New line after each row + ENDDO - is_hermitian = all(abs(A_prime_sym - - $ transpose(conjg(A_prime_sym))) < tolerance) - IF (is_hermitian) then - WRITE(*,*), "A_prime is Hermitian" - ELSE - WRITE(*,*), "A_prime is not Hermitian" - END IF - is_hermitian = all(abs(Gamma_prime_sym - - $ transpose(conjg(B_prime_sym))) < tolerance) - IF (is_hermitian) then - WRITE(*,*), "B_prime is Hermitian" - ELSE - WRITE(*,*), "B_prime is not Hermitian" - END IF - is_hermitian = all(abs(B_prime_sym - - $ transpose(conjg(Gamma_prime_sym))) < tolerance) - IF (is_hermitian) then - WRITE(*,*), "Gamma_prime_sym is Hermitian" - ELSE - WRITE(*,*), "Gamma_prime_sym is not Hermitian" - END IF is_hermitian = all(abs(Delta_prime_sym - $ transpose(conjg(Delta_prime_sym))) < tolerance) IF (is_hermitian) then