From 59e481d6659bc3b2b96d77f165bdf2338ba123e9 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Tue, 26 Nov 2024 17:38:16 -0800 Subject: [PATCH 1/8] new LE3 GK branch --- cgyro/src/cgyro_equilibrium.F90 | 16 +++++ le3/src/le3_geometry.f90 | 8 +++ le3/src/le3_geometry_matrix.f90 | 112 +++++++++++++++++++++++++++++--- le3/src/le3_geometry_rho.f90 | 12 ++++ le3/src/le3_globals.f90 | 3 + 5 files changed, 141 insertions(+), 10 deletions(-) diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index 453ea1681..993b828d5 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/cgyro/src/cgyro_equilibrium.F90 @@ -203,6 +203,22 @@ subroutine cgyro_equilibrium ! 2. Compute terms required for O(M) rotation, and no rotation. ! + ! EAB LE3 test + open(unit=1,file='out.cgyro.geogk',status='replace') + write (1,'(i2)') n_theta + write (1,'(e16.8)') theta(:) + write (1,'(e16.8)') geo_b(:) + write (1,'(e16.8)') 1.0/(q*rmaj*geo_g_theta(:)) + write (1,'(e16.8)') (geo_dbdt(:)/geo_b(:))/(q*rmaj*geo_g_theta(:)) + write (1,'(e16.8)') -1.0/geo_b(:)*geo_grad_r(:)/rmaj*geo_gsin(:) + write (1,'(e16.8)') -1.0/geo_b(:)*geo_gq(:)/rmaj& + *(geo_gcos1(:)+geo_gcos2(:)+geo_captheta(:)*geo_gsin(:)) + write (1,'(e16.8)') 1.0/geo_b(:)*geo_gq(:)/rmaj*geo_gcos2(:) + write (1,'(e16.8)') geo_grad_r(:)**2 + write (1,'(e16.8)') geo_gq(:)**2 * (1.0+geo_captheta(:)**2) + write (1,'(e16.8)') 2.0*geo_grad_r(:)*geo_gq(:)*geo_captheta(:) + close(1) + do it=1,n_theta do is=1,n_species diff --git a/le3/src/le3_geometry.f90 b/le3/src/le3_geometry.f90 index ce1eaefd1..e6045975d 100644 --- a/le3/src/le3_geometry.f90 +++ b/le3/src/le3_geometry.f90 @@ -46,6 +46,14 @@ subroutine le3_geometry allocate(vdrift_x(nt,np)) allocate(vdrift_dt(nt,np)) allocate(vdrift_dp(nt,np)) + allocate(vdrift_gk(nt,np,3)) + allocate(grad_perpsq_gk(nt,np,3)) + allocate(grad_cc(nt,np)) + allocate(grad_tt(nt,np)) + allocate(grad_pp(nt,np)) + allocate(grad_pt(nt,np)) + allocate(grad_ct(nt,np)) + allocate(grad_cp(nt,np)) allocate(vexb_dt(nt,np)) allocate(vexb_dp(nt,np)) allocate(dgdp(nt,np)) diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90 index f8efac05e..9d4564ab3 100644 --- a/le3/src/le3_geometry_matrix.f90 +++ b/le3/src/le3_geometry_matrix.f90 @@ -20,26 +20,70 @@ subroutine le3_geometry_matrix real :: xsum integer, dimension(:), allocatable :: i_piv - ! bhat dot grad = bdotgrad * (iota d/dt + d/dp) + ! anorm* bhat dot grad = bdotgrad * (iota d/dt + d/dp) bdotgrad(:,:) = 1.0/(bmag * g) - ! (bhat dot grad B)/B + ! anorm * (bhat dot grad B)/B bdotgradB_overB(:,:) = bdotgrad * (iota * dbdt + dbdp) / bmag - ! bhat cross grad B dot grad r / B^2 - vdrift_x(:,:) = 1/(rmin*bmag * g**2) & + ! B cross grad B dot grad chi / (B^2 * r * B) + vdrift_x(:,:) = 1/(rmin*bmag* g**2) & * (-dbdt * (gpp + iota * gpt) + dbdp * (gpt + iota*gtt)) / bmag**2 - ! bhat cross grad B dot grad theta / B^2 + ! B cross grad B dot grad theta / B^2 + ! EAB: is there (1/bmag) missing here?? vdrift_dt(:,:) = 1/g**2 & * ((b1*bmag/chi1 - dtheta*dbdt)*(gpp + iota * gpt) & - dbdp*gc) / bmag**2 - ! bhat cross grad B dot grad phi / B^2 + ! B cross grad B dot grad phi / B^2 + ! EAB: is there (1/bmag) missing here?? vdrift_dp(:,:) = 1/g**2 & * (-(b1*bmag/chi1 - dtheta*dbdt)*(gpt + iota * gtt) & + dbdt*gc) / bmag**2 - + + ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)] + vdrift_gk(:,:,1) = vdrift_x(:,:) + vdrift_gk(:,:,2) = (-iota*vdrift_dp(:,:) + vdrift_dt(:,:))/bmag(:,:) + do i=1,nt + do j=1,np + vdrift_gk(i,j,2) = vdrift_gk(i,j,2) & + - vdrift_x(i,j) * rmin*(iota_p/iota)*t(i) + enddo + enddo + vdrift_gk(:,:,2) = vdrift_gk(:,:,2)*rmin + i=1 + j=1 + print *, vdrift_gk(i,j,2), b1(i,j)*bmag(i,j) + print *, gpp(i,j), gpt(i,j), gtt(i,j), dbdp(i,j), dbdt(i,j), gc(i,j) + + ! (a/cs)*vdrift_gradp = 1/(cs*anorm*Omega_ca_unit)*(vpar^2) [(3) (i k_theta a)] + vdrift_gk(:,:,3) = (-0.5*beta_star)/bmag(:,:)**2 + + ! grad_perpsq: (1) d^2/d(r/a)^2 + (2) (k_theta a)^2 + (3) (i k_theta a) d/d(r/a) + grad_cc = 1.0/g**2 * (gtt * gpp - gpt**2) + grad_tt = 1.0/g**2 * (gpp * gcc - gcp**2) + grad_pp = 1.0/g**2 * (gtt * gcc - gct**2) + grad_pt = 1.0/g**2 * (gcp * gct - gpt*gcc) + grad_ct = 1.0/g**2 * (gpt * gcp - gpp*gct) + grad_cp = 1.0/g**2 * (gpt * gct - gtt*gcp) + grad_perpsq_gk(:,:,1) = 1.0/rmin**2 * grad_cc(:,:) + grad_perpsq_gk(:,:,2) = grad_pp + 1.0/iota**2*grad_tt - 2.0/iota*grad_pt + do i=1,nt + do j=1,np + grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + (iota_p/iota**2*t(i))**2 * grad_cc(i,j) & + + 2.0*iota_p/iota**2*t(i)*grad_cp(i,j) - 2.0*iota_p/iota**3*t(i)*grad_ct(i,j) + enddo + enddo + grad_perpsq_gk(:,:,2) = grad_perpsq_gk(:,:,2)*(iota*rmin)**2 + grad_perpsq_gk(:,:,3) = grad_cp - 1.0/iota*grad_ct + do i=1,nt + do j=1,np + grad_perpsq_gk(i,j,3) = grad_perpsq_gk(i,j,3) + iota_p/iota**2*t(i)* grad_cc(i,j) + enddo + enddo + grad_perpsq_gk(:,:,3) = grad_perpsq_gk(:,:,3)*(-2.0*iota) + ! flux-surface d volume / dr vprime = 0.0 do i=1,nt @@ -64,7 +108,7 @@ subroutine le3_geometry_matrix vexb_dp(:,:) = -1/(rmin*bmag * g**2) & * (gpt + iota*gtt) * bmag/bsq_avg - ! construct the geo collocation matices + ! construct the geo collocation matices for NEO-3D allocate(mat_stream_dt(matsize,matsize)) allocate(mat_stream_dp(matsize,matsize)) @@ -174,7 +218,7 @@ subroutine le3_geometry_matrix enddo close(1) - ! Construct the geo vectors + ! Construct the geo vectors for NEO-3D vec_thetabar(:) = vec_thetabar(:) / (nt*np) vec_vdriftx(:) = vec_vdriftx(:) / (nt*np) @@ -212,7 +256,7 @@ subroutine le3_geometry_matrix write (1,'(i2,i2,i2)') itype(i),m_indx(i), n_indx(i) enddo close(1) - + ! Map to Boozer coordinates ! J: B ~ J(psi) grad phi @@ -412,6 +456,43 @@ subroutine le3_geometry_matrix call le3_compute_theory + ! Write out quantitites for GK-3D + + open(unit=1,file='out.le3.geogk',status='replace') + write (1,'(i2)') nt + write (1,'(i2)') np + write (1,'(e16.8)') t + write (1,'(e16.8)') p + + ! theta_bar(theta,phi) and derivatives + write (1,'(e16.8)') tb(:,:) + write (1,'(e16.8)') dtbdt(:,:) + write (1,'(e16.8)') dtbdp(:,:) + + ! B/Bunit + write (1,'(e16.8)') bmag(:,:) + + ! anorm*(bhat dot grad) = bdotgrad * (iota d/dt + d/dp) = (1) d/dt + (2) d/dp + write (1,'(e16.8)') bdotgrad(:,:)*iota + write (1,'(e16.8)') bdotgrad(:,:) + + ! anorm*(bhat dot grad B)/B + write (1,'(e16.8)') bdotgradB_overB(:,:) + + ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)] + write (1,'(e16.8)') vdrift_gk(:,:,1) + write (1,'(e16.8)') vdrift_gk(:,:,2) + + ! (a/cs)*vdrift_gradp = 1/(cs*anorm*Omega_ca_unit)*(par^2)*(i k_theta a) [(3) (i ktheta a)] + write (1,'(e16.8)') vdrift_gk(:,:,3) + + ! grad_perpsq: (1) d^2/d(r/a)^2 + (2) (k_theta a)^2 + (3) (i k_theta a) d/d(r/a) + write (1,'(e16.8)') grad_perpsq_gk(:,:,1) + write (1,'(e16.8)') grad_perpsq_gk(:,:,2) + write (1,'(e16.8)') grad_perpsq_gk(:,:,3) + + close(1) + deallocate(itype) deallocate(m_indx) deallocate(n_indx) @@ -423,6 +504,9 @@ subroutine le3_geometry_matrix deallocate(gpp) deallocate(gtt) deallocate(gpt) + deallocate(gct) + deallocate(gcp) + deallocate(gcc) deallocate(cosu) deallocate(btor) deallocate(bpol) @@ -440,6 +524,14 @@ subroutine le3_geometry_matrix deallocate(vexb_dp) deallocate(vdrift_dt) deallocate(vdrift_dp) + deallocate(vdrift_gk) + deallocate(grad_perpsq_gk) + deallocate(grad_cc) + deallocate(grad_tt) + deallocate(grad_pp) + deallocate(grad_pt) + deallocate(grad_ct) + deallocate(grad_cp) deallocate(mat_stream_dt) deallocate(mat_stream_dp) deallocate(mat_trap) diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index 0a2db3f3a..79e6001e2 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -45,6 +45,9 @@ subroutine le3_geometry_rho allocate(dthetap(nt,np)) allocate(b1(nt,np)) allocate(gc(nt,np)) + allocate(gct(nt,np)) + allocate(gcp(nt,np)) + allocate(gcc(nt,np)) do i=1,nt do j=1,np @@ -234,6 +237,15 @@ subroutine le3_geometry_rho ! g_cp + i g_ct gc(i,j) = ysinuv/chi1(i,j)-c_m0*dtheta(i,j) + ! g_ct + gct(i,j) = -dtheta(i,j) * x**2 + + ! g_cp + gcp(i,j) = ysinuv/chi1(i,j) - dtheta(i,j)*x*ycosuv + + ! g_cc + gcc(i,j) = 1.0/chi1(i,j)**2 + dtheta(i,j)**2 * x**2 + enddo enddo diff --git a/le3/src/le3_globals.f90 b/le3/src/le3_globals.f90 index 52e1340fa..e7465424d 100644 --- a/le3/src/le3_globals.f90 +++ b/le3/src/le3_globals.f90 @@ -38,6 +38,7 @@ module le3_globals real, dimension(:,:), allocatable :: sinn,cosn real, dimension(:,:), allocatable :: gpp,gtt,gpt + real, dimension(:,:), allocatable :: gct, gcp, gcc real, dimension(:,:), allocatable :: cosu,rc real, dimension(:,:), allocatable :: bpol,btor real, dimension(:,:), allocatable :: chi1 @@ -72,6 +73,8 @@ module le3_globals real, dimension(:,:), allocatable :: bdotgradB_overB real, dimension(:,:), allocatable :: dgdp,vexb_dt,vexb_dp real, dimension(:,:), allocatable :: vdrift_x, vdrift_dt, vdrift_dp + real, dimension(:,:,:), allocatable :: grad_perpsq_gk, vdrift_gk + real, dimension(:,:), allocatable :: grad_cc, grad_tt, grad_pp, grad_pt, grad_ct, grad_cp real, dimension(:,:), allocatable :: g real, dimension(:,:), allocatable :: mat_stream_dt real, dimension(:,:), allocatable :: mat_stream_dp From 63da77c9594a97b31b97bad41fbd4f269efa64f1 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Mon, 9 Dec 2024 10:22:14 -0800 Subject: [PATCH 2/8] le3 b1 updates --- le3/src/le3_geometry_matrix.f90 | 27 +++++++++++++-------------- le3/src/le3_geometry_rho.f90 | 27 ++++++++++++++++++++++----- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90 index 9d4564ab3..2d03e6260 100644 --- a/le3/src/le3_geometry_matrix.f90 +++ b/le3/src/le3_geometry_matrix.f90 @@ -26,25 +26,23 @@ subroutine le3_geometry_matrix ! anorm * (bhat dot grad B)/B bdotgradB_overB(:,:) = bdotgrad * (iota * dbdt + dbdp) / bmag - ! B cross grad B dot grad chi / (B^2 * r * B) + ! bhat cross grad B dot grad chi / (B^2 * r) vdrift_x(:,:) = 1/(rmin*bmag* g**2) & * (-dbdt * (gpp + iota * gpt) + dbdp * (gpt + iota*gtt)) / bmag**2 - ! B cross grad B dot grad theta / B^2 - ! EAB: is there (1/bmag) missing here?? + ! bhat cross grad B dot grad theta / B^2 vdrift_dt(:,:) = 1/g**2 & * ((b1*bmag/chi1 - dtheta*dbdt)*(gpp + iota * gpt) & - - dbdp*gc) / bmag**2 + - dbdp*gc) / bmag**3 - ! B cross grad B dot grad phi / B^2 - ! EAB: is there (1/bmag) missing here?? + ! bhat cross grad B dot grad phi / B^2 vdrift_dp(:,:) = 1/g**2 & * (-(b1*bmag/chi1 - dtheta*dbdt)*(gpt + iota * gtt) & - + dbdt*gc) / bmag**2 + + dbdt*gc) / bmag**3 ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)] vdrift_gk(:,:,1) = vdrift_x(:,:) - vdrift_gk(:,:,2) = (-iota*vdrift_dp(:,:) + vdrift_dt(:,:))/bmag(:,:) + vdrift_gk(:,:,2) = -iota*vdrift_dp(:,:) + vdrift_dt(:,:) do i=1,nt do j=1,np vdrift_gk(i,j,2) = vdrift_gk(i,j,2) & @@ -52,10 +50,11 @@ subroutine le3_geometry_matrix enddo enddo vdrift_gk(:,:,2) = vdrift_gk(:,:,2)*rmin - i=1 - j=1 - print *, vdrift_gk(i,j,2), b1(i,j)*bmag(i,j) - print *, gpp(i,j), gpt(i,j), gtt(i,j), dbdp(i,j), dbdt(i,j), gc(i,j) + !i=1 + !j=1 + !print *, vdrift_gk(i,j,2), -iota*vdrift_dp(i,j)*rmin, vdrift_dt(i,j)*rmin + !print *, rmin*b1(i,j)/chi1(i,j) + !print *, g(i,j), b1(i,j)*bmag(i,j), chi1(i,j), bmag(i,j) ! (a/cs)*vdrift_gradp = 1/(cs*anorm*Omega_ca_unit)*(vpar^2) [(3) (i k_theta a)] vdrift_gk(:,:,3) = (-0.5*beta_star)/bmag(:,:)**2 @@ -459,8 +458,8 @@ subroutine le3_geometry_matrix ! Write out quantitites for GK-3D open(unit=1,file='out.le3.geogk',status='replace') - write (1,'(i2)') nt - write (1,'(i2)') np + write (1,'(i4)') nt + write (1,'(i4)') np write (1,'(e16.8)') t write (1,'(e16.8)') p diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index 79e6001e2..1b89dcc77 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -15,7 +15,7 @@ subroutine le3_geometry_rho real :: up,pp,pt real :: d0,eta - real, dimension(:,:), allocatable :: s1a,s1b + real, dimension(:,:), allocatable :: s1a,s1b, s1c real, dimension(:,:), allocatable :: s2a,s2b real, dimension(:,:), allocatable :: a11a,a11b real, dimension(:,:), allocatable :: a12a @@ -31,6 +31,7 @@ subroutine le3_geometry_rho ! Step 1: Generate required functions on (theta,phi) mesh allocate(s1a(nt,np)) allocate(s1b(nt,np)) + allocate(s1c(nt,np)) allocate(s2a(nt,np)) allocate(s2b(nt,np)) allocate(a11a(nt,np)) @@ -77,6 +78,7 @@ subroutine le3_geometry_rho s1a(i,j) = -0.5*d0*beta_star/rmin+(c_dn+iota*c_dm)/(chi1(i,j)*d0) s1b(i,j) = ysinuv/(chi1(i,j)*d0) + s1c(i,j) = (c_dn+iota*c_dm)/d0 s2a(i,j) = c_dm/(chi1(i,j)*d0) s2b(i,j) = c_dn/(chi1(i,j)*d0) @@ -231,9 +233,24 @@ subroutine le3_geometry_rho !------------------------------------------------------------- ! B1/Bs - b1(i,j) = 0.5*chi1(i,j)/a12a(i,j)*(iota_p*c_m0/d0+s1a(i,j)) & - -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j)) - + !b1(i,j) = 0.5*chi1(i,j)/a12a(i,j)*(iota_p*c_m0/d0+s1a(i,j)) & + ! -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j)) + + b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & + -0.5*(eta-2*dchi(i,j)*chi1(i,j)-2*chi1(i,j)*dthetap(i,j)) + + if(j == 1) then + if (i == 1) then + print *, t(i), tb(i,j) + print *, b1(i,j) + print *, 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0) + print *, 0.5/a12a(i,j)*s1c(i,j) + print *, eta + print *, d0 + print *, chi1(i,j), dchi(i,j), dthetap(i,j) + endif + endif + ! g_cp + i g_ct gc(i,j) = ysinuv/chi1(i,j)-c_m0*dtheta(i,j) @@ -282,7 +299,7 @@ subroutine le3_geometry_rho endif !-------------------------------------------------------------------- - + deallocate(sys_m) deallocate(sys_b) deallocate(i_piv) From 81bb08f4fc718b0fa3e11a8dfbcfdeabcb0faa69 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Tue, 10 Dec 2024 15:52:11 -0800 Subject: [PATCH 3/8] new b1 and delta_chi params in le3 for checking --- cgyro/src/cgyro_equilibrium.F90 | 2 ++ le3/src/le3_geometry_matrix.f90 | 4 ++++ le3/src/le3_geometry_rho.f90 | 33 +++++++++++++++++++-------------- 3 files changed, 25 insertions(+), 14 deletions(-) diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index 993b828d5..ac7351b31 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/cgyro/src/cgyro_equilibrium.F90 @@ -217,6 +217,8 @@ subroutine cgyro_equilibrium write (1,'(e16.8)') geo_grad_r(:)**2 write (1,'(e16.8)') geo_gq(:)**2 * (1.0+geo_captheta(:)**2) write (1,'(e16.8)') 2.0*geo_grad_r(:)*geo_gq(:)*geo_captheta(:) + write (1,'(e16.8)') geo_chi2(:) + write (1,'(e16.8)') -(geo_gcos1(:) + geo_gcos2(:))/rmaj close(1) do it=1,n_theta diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90 index 2d03e6260..f263cc434 100644 --- a/le3/src/le3_geometry_matrix.f90 +++ b/le3/src/le3_geometry_matrix.f90 @@ -489,6 +489,10 @@ subroutine le3_geometry_matrix write (1,'(e16.8)') grad_perpsq_gk(:,:,1) write (1,'(e16.8)') grad_perpsq_gk(:,:,2) write (1,'(e16.8)') grad_perpsq_gk(:,:,3) + + write (1,'(e16.8)') dchi(:,:) + + write (1,'(e16.8)') b1(:,:) close(1) diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index 1b89dcc77..b09c22bad 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -91,7 +91,6 @@ subroutine le3_geometry_rho a21a(i,j) = (-2*c_m0+iota*x**2)/d0 a21b(i,j) = x**2/d0 - enddo enddo @@ -237,19 +236,25 @@ subroutine le3_geometry_rho ! -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j)) b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & - -0.5*(eta-2*dchi(i,j)*chi1(i,j)-2*chi1(i,j)*dthetap(i,j)) - - if(j == 1) then - if (i == 1) then - print *, t(i), tb(i,j) - print *, b1(i,j) - print *, 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0) - print *, 0.5/a12a(i,j)*s1c(i,j) - print *, eta - print *, d0 - print *, chi1(i,j), dchi(i,j), dthetap(i,j) - endif - endif + -0.5*(eta-2.0*dchi(i,j)*chi1(i,j)-2.0*chi1(i,j)*dthetap(i,j)) + + print *, dchi(i,j)+dthetap(i,j), & + -iota*(c_m0/c_n0)*dchi(i,j) - 1.0/(chi1(i,j)*c_n0)*(c_dn + iota*c_dm), & + chi1(i,j), eta + + !if(j == 1) then + ! if (i == 1) then + ! print *, t(i), tb(i,j) + ! print *, b1(i,j) + ! print *, (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j))/(r(i,j)**2 + iota**2 * x**2) + ! print *, rc(i,j), rmin + !print *, b1(i,j), & + ! -eta + 0.5*chi1(i,j)*(dchi(i,j)+dthetap(i,j)) & + ! + (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j)& + ! + chi1(i,j)*iota_p*iota*x**2) & + ! /(r(i,j)**2 + iota**2 * x**2) + ! endif + !endif ! g_cp + i g_ct gc(i,j) = ysinuv/chi1(i,j)-c_m0*dtheta(i,j) From e0041c5ace38eababe079fdf813850e8587ecb70 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Thu, 12 Dec 2024 17:05:00 -0800 Subject: [PATCH 4/8] le3 updates --- cgyro/src/cgyro_equilibrium.F90 | 1 + le3/src/le3_geometry_matrix.f90 | 14 +++++++++++++- le3/src/le3_geometry_rho.f90 | 23 ++++++++++++++++++----- le3/src/le3_globals.f90 | 2 +- 4 files changed, 33 insertions(+), 7 deletions(-) diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index ac7351b31..a08b83741 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/cgyro/src/cgyro_equilibrium.F90 @@ -204,6 +204,7 @@ subroutine cgyro_equilibrium ! ! EAB LE3 test + print *, geo_f open(unit=1,file='out.cgyro.geogk',status='replace') write (1,'(i2)') n_theta write (1,'(e16.8)') theta(:) diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90 index f263cc434..08abb9957 100644 --- a/le3/src/le3_geometry_matrix.f90 +++ b/le3/src/le3_geometry_matrix.f90 @@ -490,9 +490,21 @@ subroutine le3_geometry_matrix write (1,'(e16.8)') grad_perpsq_gk(:,:,2) write (1,'(e16.8)') grad_perpsq_gk(:,:,3) + write (1,'(e16.8)') b1(:,:) + write (1,'(e16.8)') dchi(:,:) - write (1,'(e16.8)') b1(:,:) + write (1,'(e16.8)') dthetap(:,:) + + write (1,'(e16.8)') chi1(:,:) + + write (1,'(e16.8)') b1_temp1(:,:) + + write (1,'(e16.8)') b1_temp2(:,:) + + write (1,'(e16.8)') b1_temp3(:,:) + + write (1,'(e16.8)') dtheta(:,:) close(1) diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index b09c22bad..d8be4ef80 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -40,6 +40,10 @@ subroutine le3_geometry_rho allocate(a21a(nt,np)) allocate(a21b(nt,np)) + allocate(b1_temp1(nt,np)) + allocate(b1_temp2(nt,np)) + allocate(b1_temp3(nt,np)) + ! Global second order quantities allocate(dtheta(nt,np)) allocate(dchi(nt,np)) @@ -236,14 +240,23 @@ subroutine le3_geometry_rho ! -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j)) b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & - -0.5*(eta-2.0*dchi(i,j)*chi1(i,j)-2.0*chi1(i,j)*dthetap(i,j)) - - print *, dchi(i,j)+dthetap(i,j), & - -iota*(c_m0/c_n0)*dchi(i,j) - 1.0/(chi1(i,j)*c_n0)*(c_dn + iota*c_dm), & - chi1(i,j), eta + -0.5*(eta-2*dchi(i,j)*chi1(i,j)-2*chi1(i,j)*dthetap(i,j)) + + !print *, chi1(i,j)*(dchi(i,j)+dthetap(i,j)), & + ! -chi1(i,j)*iota**2 *x**2/r(i,j)**2 * dchi(i,j) & + ! - 1.0/(r(i,j)**2)*(c_dn + iota*c_dm), & + ! chi1(i,j)*(dchi(i,j)+dthetap(i,j))+c_dn/c_n0, & + ! rc(i,j),x**2 + ! chi1(i,j), eta + print *, t(i), c_n0/d0 + + b1_temp1(i,j) = eta + b1_temp2(i,j) = 1.0/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0) + b1_temp3(i,j) = 1.0/a12a(i,j)*s1c(i,j) !if(j == 1) then ! if (i == 1) then + ! print *, chi1(i,j), d0 ! print *, t(i), tb(i,j) ! print *, b1(i,j) ! print *, (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j))/(r(i,j)**2 + iota**2 * x**2) diff --git a/le3/src/le3_globals.f90 b/le3/src/le3_globals.f90 index e7465424d..417906f41 100644 --- a/le3/src/le3_globals.f90 +++ b/le3/src/le3_globals.f90 @@ -53,7 +53,7 @@ module le3_globals real, dimension(:,:), allocatable :: dtheta real, dimension(:,:), allocatable :: dchi real, dimension(:,:), allocatable :: dthetap - real, dimension(:,:), allocatable :: b1 + real, dimension(:,:), allocatable :: b1, b1_temp1, b1_temp2, b1_temp3 real, dimension(:,:), allocatable :: gc !-------------------------------------------------------- From fa33724d024bf2fba200eae7fb9d33176bb72a10 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Wed, 18 Dec 2024 18:36:58 -0800 Subject: [PATCH 5/8] fixed LE3 B1 for 2D limit; still need 3D fix --- cgyro/src/cgyro_equilibrium.F90 | 2 +- le3/src/le3_geometry_rho.f90 | 49 +++++++++++++++++++++++++-------- 2 files changed, 39 insertions(+), 12 deletions(-) diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index a08b83741..db50eb49a 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/cgyro/src/cgyro_equilibrium.F90 @@ -204,7 +204,7 @@ subroutine cgyro_equilibrium ! ! EAB LE3 test - print *, geo_f + print *, geo_f, geo_ffprime/geo_f open(unit=1,file='out.cgyro.geogk',status='replace') write (1,'(i2)') n_theta write (1,'(e16.8)') theta(:) diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index d8be4ef80..fa003b33f 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -14,6 +14,7 @@ subroutine le3_geometry_rho real :: ysinuv real :: up,pp,pt real :: d0,eta + real :: sum1, sum2, sum3, sum4 real, dimension(:,:), allocatable :: s1a,s1b, s1c real, dimension(:,:), allocatable :: s2a,s2b @@ -236,19 +237,15 @@ subroutine le3_geometry_rho !------------------------------------------------------------- ! B1/Bs - !b1(i,j) = 0.5*chi1(i,j)/a12a(i,j)*(iota_p*c_m0/d0+s1a(i,j)) & - ! -0.5*(eta-dchi(i,j)-chi1(i,j)*dthetap(i,j)) - b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & - -0.5*(eta-2*dchi(i,j)*chi1(i,j)-2*chi1(i,j)*dthetap(i,j)) + !b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & + ! -0.5*(eta-2.0*dchi(i,j)*chi1(i,j)-2.0*chi1(i,j)*dthetap(i,j)) - !print *, chi1(i,j)*(dchi(i,j)+dthetap(i,j)), & - ! -chi1(i,j)*iota**2 *x**2/r(i,j)**2 * dchi(i,j) & - ! - 1.0/(r(i,j)**2)*(c_dn + iota*c_dm), & - ! chi1(i,j)*(dchi(i,j)+dthetap(i,j))+c_dn/c_n0, & - ! rc(i,j),x**2 - ! chi1(i,j), eta - print *, t(i), c_n0/d0 + ! EAB temp + b1(i,j) = -eta+chi1(i,j)*(dchi(i,j)+dthetap(i,j)) & + + 1.0/(c_n0 + iota*c_m0) * (r(i,j)*cosu(i,j) & + + iota**2 * x**2 * (1.0/rc(i,j)- chi1(i,j)*dthetap(i,j)) & + + chi1(i,j)*iota*iota_p*x**2) b1_temp1(i,j) = eta b1_temp2(i,j) = 1.0/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0) @@ -284,6 +281,36 @@ subroutine le3_geometry_rho enddo enddo + j=1 + sum1=0.0 + sum2=0.0 + sum3=0.0 + sum4=0.0 + do i=1,nt + sum4 = sum4 + (cos(tb(i,j)))/r(i,j)-1/rmin + sum1 = sum1 + cos(tb(i,j))/r(i,j) + sum2 = sum2 + 1/(r(i,j))*dtbdt(i,j) + sum3 = sum3 + dchi(i,j) + enddo + sum1 = sum1/nt + sum2 = sum2/nt + sum3 = sum3/nt + sum4 = sum4/nt + print *, (-rmin**2 * sum2**2 * chi1(i,1)*iota_p*iota & + -2*rmin**2 * sum2**2 *iota**2 * (-sum4)) & + / (1.0 + iota**2 * rmin**2 * sum2**2) + do i=1,nt + x = sqrt(gtt(i,j)) + print *, b1(i,j), -1.0/(r(i,j)**2 + iota**2 * x**2) & + * (r(i,j)*cos(tb(i,j)) + iota**2*x**2/rmin & + + 0.0*2*iota**2 * x**2*(cos(tb(i,j))/r(i,j) - sum1)) + !print *, chi1(i,j)*r(i,j)/x, r(i,j)/x/iota & + ! *(-1/rc(i,j)+cosu(i,j)/r(i,j)+chi1(i,j)*(dchi(i,j)+dthetap(i,j))) + enddo + !print *, chi1(i,1)*(1+iota**2 * gtt(i,1)/r(i,1)**2)*sum3, & + ! - sum1 + !print *, sum4 - sum1/(1.0 + iota**2 * rmin**2 * sum2**2) + !-------------------------------------------------------------------- ! Check with GEO result ! From 8d2c685e23f88d1eded67280eb5aab1f0fbb4aa4 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Thu, 2 Jan 2025 16:11:25 -0800 Subject: [PATCH 6/8] re-mapping theta in le3 --- le3/src/le3_geometry_matrix.f90 | 126 +++++++++++++++++++++++--------- le3/src/le3_geometry_rho.f90 | 106 +++++++++------------------ le3/src/le3_globals.f90 | 4 +- 3 files changed, 129 insertions(+), 107 deletions(-) diff --git a/le3/src/le3_geometry_matrix.f90 b/le3/src/le3_geometry_matrix.f90 index 08abb9957..15d612553 100644 --- a/le3/src/le3_geometry_matrix.f90 +++ b/le3/src/le3_geometry_matrix.f90 @@ -2,7 +2,7 @@ subroutine le3_geometry_matrix use le3_globals implicit none - integer :: i,j,its,ips,kt,kp + integer :: i,j,k,its,ips,kt,kp real, dimension(:), allocatable :: vec_vdriftx real, dimension(:), allocatable :: vec_flux @@ -17,6 +17,8 @@ subroutine le3_geometry_matrix real :: J_boozer, I_boozer real, dimension(:), allocatable :: vec_boozer_p, vec_boozer_t, vec_boozer real, dimension(:,:), allocatable :: mat_boozer, t_boozer, p_boozer + real, dimension(:,:), allocatable :: xtemp + real, dimension(:), allocatable :: tmap real :: xsum integer, dimension(:), allocatable :: i_piv @@ -41,24 +43,35 @@ subroutine le3_geometry_matrix + dbdt*gc) / bmag**3 ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)] + ! Remap theta: 0..2pi -> -pi..pi + + allocate(xtemp(nt,np)) + allocate(tmap(nt)) + do i=1,nt/2 + k=nt/2+i + tmap(i) = t(k)-2*pi + tmap(k) = t(i) + enddo + vdrift_gk(:,:,1) = vdrift_x(:,:) - vdrift_gk(:,:,2) = -iota*vdrift_dp(:,:) + vdrift_dt(:,:) + call remap_theta(vdrift_gk(:,:,1),xtemp(:,:),0.0) + vdrift_gk(:,:,1) = xtemp(:,:) + + vdrift_gk(:,:,2) = (-iota*vdrift_dp(:,:) + vdrift_dt(:,:))*rmin + call remap_theta(vdrift_gk(:,:,2),xtemp(:,:),0.0) + vdrift_gk(:,:,2) = xtemp(:,:) do i=1,nt do j=1,np vdrift_gk(i,j,2) = vdrift_gk(i,j,2) & - - vdrift_x(i,j) * rmin*(iota_p/iota)*t(i) + - vdrift_gk(i,j,1) * rmin**2 *(iota_p/iota)*tmap(i) enddo enddo - vdrift_gk(:,:,2) = vdrift_gk(:,:,2)*rmin - !i=1 - !j=1 - !print *, vdrift_gk(i,j,2), -iota*vdrift_dp(i,j)*rmin, vdrift_dt(i,j)*rmin - !print *, rmin*b1(i,j)/chi1(i,j) - !print *, g(i,j), b1(i,j)*bmag(i,j), chi1(i,j), bmag(i,j) ! (a/cs)*vdrift_gradp = 1/(cs*anorm*Omega_ca_unit)*(vpar^2) [(3) (i k_theta a)] vdrift_gk(:,:,3) = (-0.5*beta_star)/bmag(:,:)**2 - + call remap_theta(vdrift_gk(:,:,3),xtemp(:,:),0.0) + vdrift_gk(:,:,3) = xtemp(:,:) + ! grad_perpsq: (1) d^2/d(r/a)^2 + (2) (k_theta a)^2 + (3) (i k_theta a) d/d(r/a) grad_cc = 1.0/g**2 * (gtt * gpp - gpt**2) grad_tt = 1.0/g**2 * (gpp * gcc - gcp**2) @@ -66,19 +79,41 @@ subroutine le3_geometry_matrix grad_pt = 1.0/g**2 * (gcp * gct - gpt*gcc) grad_ct = 1.0/g**2 * (gpt * gcp - gpp*gct) grad_cp = 1.0/g**2 * (gpt * gct - gtt*gcp) + grad_perpsq_gk(:,:,1) = 1.0/rmin**2 * grad_cc(:,:) - grad_perpsq_gk(:,:,2) = grad_pp + 1.0/iota**2*grad_tt - 2.0/iota*grad_pt + call remap_theta(grad_perpsq_gk(:,:,1),xtemp(:,:),0.0) + grad_perpsq_gk(:,:,1) = xtemp(:,:) + + grad_perpsq_gk(:,:,2) = grad_pp + 1.0/iota**2*grad_tt - 2.0/iota*grad_pt + call remap_theta(grad_perpsq_gk(:,:,2),xtemp(:,:),0.0) + grad_perpsq_gk(:,:,2) = xtemp(:,:) + call remap_theta(grad_cc(:,:),xtemp(:,:),0.0) do i=1,nt do j=1,np - grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + (iota_p/iota**2*t(i))**2 * grad_cc(i,j) & - + 2.0*iota_p/iota**2*t(i)*grad_cp(i,j) - 2.0*iota_p/iota**3*t(i)*grad_ct(i,j) + grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + (iota_p/iota**2*tmap(i))**2 * xtemp(i,j) + enddo + enddo + call remap_theta(grad_cp(:,:),xtemp(:,:),0.0) + do i=1,nt + do j=1,np + grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) + 2.0*iota_p/iota**2*tmap(i)*xtemp(i,j) + enddo + enddo + call remap_theta(grad_ct(:,:),xtemp(:,:),0.0) + do i=1,nt + do j=1,np + grad_perpsq_gk(i,j,2) = grad_perpsq_gk(i,j,2) - 2.0*iota_p/iota**3*tmap(i)*xtemp(i,j) enddo enddo grad_perpsq_gk(:,:,2) = grad_perpsq_gk(:,:,2)*(iota*rmin)**2 + grad_perpsq_gk(:,:,3) = grad_cp - 1.0/iota*grad_ct + call remap_theta(grad_perpsq_gk(:,:,3),xtemp(:,:),0.0) + grad_perpsq_gk(:,:,3) = xtemp(:,:) + call remap_theta(grad_cc(:,:),xtemp(:,:),0.0) do i=1,nt do j=1,np - grad_perpsq_gk(i,j,3) = grad_perpsq_gk(i,j,3) + iota_p/iota**2*t(i)* grad_cc(i,j) + grad_perpsq_gk(i,j,3) = grad_perpsq_gk(i,j,3) + iota_p/iota**2*tmap(i)* xtemp(i,j) enddo enddo grad_perpsq_gk(:,:,3) = grad_perpsq_gk(:,:,3)*(-2.0*iota) @@ -456,27 +491,34 @@ subroutine le3_geometry_matrix call le3_compute_theory ! Write out quantitites for GK-3D - + ! remap theta: 0..2pi -> -pi..pi + open(unit=1,file='out.le3.geogk',status='replace') write (1,'(i4)') nt write (1,'(i4)') np - write (1,'(e16.8)') t + write (1,'(e16.8)') tmap write (1,'(e16.8)') p ! theta_bar(theta,phi) and derivatives - write (1,'(e16.8)') tb(:,:) - write (1,'(e16.8)') dtbdt(:,:) - write (1,'(e16.8)') dtbdp(:,:) + call remap_theta(tb(:,:),xtemp(:,:),1.0) + write (1,'(e16.8)') xtemp(:,:) + call remap_theta(dtbdt(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) + call remap_theta(dtbdp(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) ! B/Bunit - write (1,'(e16.8)') bmag(:,:) + call remap_theta(bmag(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) ! anorm*(bhat dot grad) = bdotgrad * (iota d/dt + d/dp) = (1) d/dt + (2) d/dp - write (1,'(e16.8)') bdotgrad(:,:)*iota - write (1,'(e16.8)') bdotgrad(:,:) + call remap_theta(bdotgrad(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:)*iota + write (1,'(e16.8)') xtemp(:,:) ! anorm*(bhat dot grad B)/B - write (1,'(e16.8)') bdotgradB_overB(:,:) + call remap_theta(bdotgradB_overB(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) ! (a/cs)*vdrift_gradB = 1/(cs*anorm*Omega_ca_unit)*(vperp^2/2+vpar^2)*[(1) d/d(r/a) + (2) (i k_theta a)] write (1,'(e16.8)') vdrift_gk(:,:,1) @@ -490,23 +532,24 @@ subroutine le3_geometry_matrix write (1,'(e16.8)') grad_perpsq_gk(:,:,2) write (1,'(e16.8)') grad_perpsq_gk(:,:,3) - write (1,'(e16.8)') b1(:,:) - - write (1,'(e16.8)') dchi(:,:) - - write (1,'(e16.8)') dthetap(:,:) - - write (1,'(e16.8)') chi1(:,:) + call remap_theta(b1(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) - write (1,'(e16.8)') b1_temp1(:,:) + call remap_theta(chi1(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) - write (1,'(e16.8)') b1_temp2(:,:) + call remap_theta(dchi(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) - write (1,'(e16.8)') b1_temp3(:,:) + call remap_theta(dtheta(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) - write (1,'(e16.8)') dtheta(:,:) + call remap_theta(dtheta_t(:,:),xtemp(:,:),0.0) + write (1,'(e16.8)') xtemp(:,:) close(1) + deallocate(xtemp) + deallocate(tmap) deallocate(itype) deallocate(m_indx) @@ -565,3 +608,18 @@ subroutine le3_geometry_matrix deallocate(vec_ntv) end subroutine le3_geometry_matrix + +subroutine remap_theta(xold,xnew,fac) + use le3_globals + implicit none + real, intent(inout), dimension(nt,np) :: xold,xnew + real, intent(in) :: fac + integer :: i,j,k + do i=1,nt/2 + k=nt/2+i + do j=1,np + xnew(i,j) = xold(k,j) - fac*2*pi + xnew(k,j) = xold(i,j) + enddo + enddo +end subroutine remap_theta diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index fa003b33f..3d44c62d7 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -13,10 +13,10 @@ subroutine le3_geometry_rho real :: ycosuv real :: ysinuv real :: up,pp,pt - real :: d0,eta + real :: d0, eta, d1, gpp_1, gtt_1, gpt_1 real :: sum1, sum2, sum3, sum4 - real, dimension(:,:), allocatable :: s1a,s1b, s1c + real, dimension(:,:), allocatable :: s1a,s1b real, dimension(:,:), allocatable :: s2a,s2b real, dimension(:,:), allocatable :: a11a,a11b real, dimension(:,:), allocatable :: a12a @@ -32,7 +32,6 @@ subroutine le3_geometry_rho ! Step 1: Generate required functions on (theta,phi) mesh allocate(s1a(nt,np)) allocate(s1b(nt,np)) - allocate(s1c(nt,np)) allocate(s2a(nt,np)) allocate(s2b(nt,np)) allocate(a11a(nt,np)) @@ -41,14 +40,11 @@ subroutine le3_geometry_rho allocate(a21a(nt,np)) allocate(a21b(nt,np)) - allocate(b1_temp1(nt,np)) - allocate(b1_temp2(nt,np)) - allocate(b1_temp3(nt,np)) - ! Global second order quantities allocate(dtheta(nt,np)) allocate(dchi(nt,np)) - allocate(dthetap(nt,np)) + allocate(dtheta_t(nt,np)) + allocate(dtheta_p(nt,np)) allocate(b1(nt,np)) allocate(gc(nt,np)) allocate(gct(nt,np)) @@ -71,7 +67,7 @@ subroutine le3_geometry_rho pt = (x/rc(i,j))*ycosuv-chi1t(i,j)/chi1(i,j)*ysinuv eta = 1.0/rc(i,j)+cosu(i,j)/r(i,j) - + c_m0 = iota*x**2+x*ycosuv c_n0 = r(i,j)**2+y**2+iota*x*ycosuv @@ -83,7 +79,6 @@ subroutine le3_geometry_rho s1a(i,j) = -0.5*d0*beta_star/rmin+(c_dn+iota*c_dm)/(chi1(i,j)*d0) s1b(i,j) = ysinuv/(chi1(i,j)*d0) - s1c(i,j) = (c_dn+iota*c_dm)/d0 s2a(i,j) = c_dm/(chi1(i,j)*d0) s2b(i,j) = c_dn/(chi1(i,j)*d0) @@ -189,17 +184,18 @@ subroutine le3_geometry_rho dtheta = 0.0 dchi = 0.0 - dthetap = 0.0 + dtheta_t = 0.0 + dtheta_p = 0.0 do k=1,matsize call le3_basis(itype(k),m_indx(k),n_indx(k),bk(:,:),'d0') call le3_basis(itype(k),m_indx(k),n_indx(k),bk_t(:,:),'dt') + call le3_basis(itype(k),m_indx(k),n_indx(k),bk_p(:,:),'dp') do j=1,np do i=1,nt - - dtheta(i,j) = dtheta(i,j) + sys_b(k)*bk(i,j) - dchi(i,j) = dchi(i,j) + sys_b(k+matsize)*bk(i,j) - - dthetap(i,j) = dthetap(i,j) + sys_b(k)*bk_t(i,j) + dchi(i,j) = dchi(i,j) + sys_b(k+matsize)*bk(i,j) + dtheta(i,j) = dtheta(i,j) + sys_b(k)*bk(i,j) + dtheta_t(i,j) = dtheta_t(i,j) + sys_b(k)*bk_t(i,j) + dtheta_p(i,j) = dtheta_p(i,j) + sys_b(k)*bk_p(i,j) enddo enddo @@ -233,38 +229,26 @@ subroutine le3_geometry_rho c_dm = x*up+pt+iota*x**2*2.0/rc(i,j)+chi1(i,j)*iota_p*x**2-c_m0*eta c_dn = 2*r(i,j)*cosu(i,j)+iota*x*up+2*pp+iota*pt+chi1(i,j)*iota_p*x*ycosuv-c_n0*eta + ! D0 = sqrt(g_s) d0 = g(i,j) - !------------------------------------------------------------- - ! B1/Bs + ! D1/D0 + d1 = eta - chi1(i,j)*(dchi(i,j)+dtheta_t(i,j)) - !b1(i,j) = 0.5/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0+s1c(i,j)) & - ! -0.5*(eta-2.0*dchi(i,j)*chi1(i,j)-2.0*chi1(i,j)*dthetap(i,j)) + gpp_1 = 2.0*(r(i,j)*cosu(i,j) - chi1p(i,j)/chi1(i,j)*ysinuv & + + (up - chi1(i,j)*x*dtheta_p(i,j))*ycosuv) - ! EAB temp - b1(i,j) = -eta+chi1(i,j)*(dchi(i,j)+dthetap(i,j)) & - + 1.0/(c_n0 + iota*c_m0) * (r(i,j)*cosu(i,j) & - + iota**2 * x**2 * (1.0/rc(i,j)- chi1(i,j)*dthetap(i,j)) & - + chi1(i,j)*iota*iota_p*x**2) + gtt_1 = 2.0 * x**2 * (1.0/rc(i,j)- chi1(i,j)*dtheta_t(i,j)) - b1_temp1(i,j) = eta - b1_temp2(i,j) = 1.0/a12a(i,j)*(chi1(i,j)*iota_p*c_m0/d0) - b1_temp3(i,j) = 1.0/a12a(i,j)*s1c(i,j) + ! Not done + gpt_1 = 0.0 - !if(j == 1) then - ! if (i == 1) then - ! print *, chi1(i,j), d0 - ! print *, t(i), tb(i,j) - ! print *, b1(i,j) - ! print *, (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j))/(r(i,j)**2 + iota**2 * x**2) - ! print *, rc(i,j), rmin - !print *, b1(i,j), & - ! -eta + 0.5*chi1(i,j)*(dchi(i,j)+dthetap(i,j)) & - ! + (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j)& - ! + chi1(i,j)*iota_p*iota*x**2) & - ! /(r(i,j)**2 + iota**2 * x**2) - ! endif - !endif + !------------------------------------------------------------- + + ! B1/Bs + b1(i,j) = -d1 & + + 0.5/(c_n0 + iota*c_m0) * (gpp_1 + iota**2 * gtt_1 + 2.0*iota*gpt_1 & + + chi1(i,j) * iota_p * (2*iota*gtt(i,j) + gpt(i,j)) ) ! g_cp + i g_ct gc(i,j) = ysinuv/chi1(i,j)-c_m0*dtheta(i,j) @@ -281,35 +265,15 @@ subroutine le3_geometry_rho enddo enddo - j=1 - sum1=0.0 - sum2=0.0 - sum3=0.0 - sum4=0.0 - do i=1,nt - sum4 = sum4 + (cos(tb(i,j)))/r(i,j)-1/rmin - sum1 = sum1 + cos(tb(i,j))/r(i,j) - sum2 = sum2 + 1/(r(i,j))*dtbdt(i,j) - sum3 = sum3 + dchi(i,j) - enddo - sum1 = sum1/nt - sum2 = sum2/nt - sum3 = sum3/nt - sum4 = sum4/nt - print *, (-rmin**2 * sum2**2 * chi1(i,1)*iota_p*iota & - -2*rmin**2 * sum2**2 *iota**2 * (-sum4)) & - / (1.0 + iota**2 * rmin**2 * sum2**2) - do i=1,nt - x = sqrt(gtt(i,j)) - print *, b1(i,j), -1.0/(r(i,j)**2 + iota**2 * x**2) & - * (r(i,j)*cos(tb(i,j)) + iota**2*x**2/rmin & - + 0.0*2*iota**2 * x**2*(cos(tb(i,j))/r(i,j) - sum1)) - !print *, chi1(i,j)*r(i,j)/x, r(i,j)/x/iota & - ! *(-1/rc(i,j)+cosu(i,j)/r(i,j)+chi1(i,j)*(dchi(i,j)+dthetap(i,j))) - enddo - !print *, chi1(i,1)*(1+iota**2 * gtt(i,1)/r(i,1)**2)*sum3, & - ! - sum1 - !print *, sum4 - sum1/(1.0 + iota**2 * rmin**2 * sum2**2) + !do i=1,nt + ! j=1 + ! x = sqrt(gtt(i,j)) + !print *, b1(i,j), -1.0/(r(i,j)**2 + iota**2 * x**2) & + ! * (r(i,j)*cosu(i,j) + iota**2*x**2/rc(i,j)), & + ! chi1(i,j)*r(i,j)/sqrt(gtt(i,j)), & + ! r(i,j)/sqrt(gtt(i,j))& + ! *(-1/rc(i,j)+cosu(i,j)/r(i,j)+chi1(i,j)*(dchi(i,j)+dtheta_t(i,j)))/iota + !enddo !-------------------------------------------------------------------- ! Check with GEO result diff --git a/le3/src/le3_globals.f90 b/le3/src/le3_globals.f90 index 417906f41..033ad57be 100644 --- a/le3/src/le3_globals.f90 +++ b/le3/src/le3_globals.f90 @@ -52,8 +52,8 @@ module le3_globals ! "Second order" quantities real, dimension(:,:), allocatable :: dtheta real, dimension(:,:), allocatable :: dchi - real, dimension(:,:), allocatable :: dthetap - real, dimension(:,:), allocatable :: b1, b1_temp1, b1_temp2, b1_temp3 + real, dimension(:,:), allocatable :: dtheta_t, dtheta_p + real, dimension(:,:), allocatable :: b1 real, dimension(:,:), allocatable :: gc !-------------------------------------------------------- From 8a91478ad44c949881784363396bfc8a8362b23e Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Thu, 2 Jan 2025 17:20:02 -0800 Subject: [PATCH 7/8] last part of b1 in le3 --- le3/src/le3_geometry_rho.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/le3/src/le3_geometry_rho.f90 b/le3/src/le3_geometry_rho.f90 index 3d44c62d7..87f6b044c 100644 --- a/le3/src/le3_geometry_rho.f90 +++ b/le3/src/le3_geometry_rho.f90 @@ -240,8 +240,9 @@ subroutine le3_geometry_rho gtt_1 = 2.0 * x**2 * (1.0/rc(i,j)- chi1(i,j)*dtheta_t(i,j)) - ! Not done - gpt_1 = 0.0 + gpt_1 = x*(up - chi1(i,j)*x*dtheta_p(i,j)) & + + (x/rc(i,j) - chi1(i,j)*x*dtheta_t(i,j))*ycosuv & + - chi1t(i,j)/chi1(i,j)*ysinuv !------------------------------------------------------------- From 6fb223e73752ed61c9f8769b989cca94360f4563 Mon Sep 17 00:00:00 2001 From: Emily Belli Date: Mon, 6 Jan 2025 15:35:31 -0800 Subject: [PATCH 8/8] CGYRO LE3 test file in comments --- cgyro/src/cgyro_equilibrium.F90 | 34 ++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/cgyro/src/cgyro_equilibrium.F90 b/cgyro/src/cgyro_equilibrium.F90 index db50eb49a..55b2493bf 100644 --- a/cgyro/src/cgyro_equilibrium.F90 +++ b/cgyro/src/cgyro_equilibrium.F90 @@ -204,23 +204,23 @@ subroutine cgyro_equilibrium ! ! EAB LE3 test - print *, geo_f, geo_ffprime/geo_f - open(unit=1,file='out.cgyro.geogk',status='replace') - write (1,'(i2)') n_theta - write (1,'(e16.8)') theta(:) - write (1,'(e16.8)') geo_b(:) - write (1,'(e16.8)') 1.0/(q*rmaj*geo_g_theta(:)) - write (1,'(e16.8)') (geo_dbdt(:)/geo_b(:))/(q*rmaj*geo_g_theta(:)) - write (1,'(e16.8)') -1.0/geo_b(:)*geo_grad_r(:)/rmaj*geo_gsin(:) - write (1,'(e16.8)') -1.0/geo_b(:)*geo_gq(:)/rmaj& - *(geo_gcos1(:)+geo_gcos2(:)+geo_captheta(:)*geo_gsin(:)) - write (1,'(e16.8)') 1.0/geo_b(:)*geo_gq(:)/rmaj*geo_gcos2(:) - write (1,'(e16.8)') geo_grad_r(:)**2 - write (1,'(e16.8)') geo_gq(:)**2 * (1.0+geo_captheta(:)**2) - write (1,'(e16.8)') 2.0*geo_grad_r(:)*geo_gq(:)*geo_captheta(:) - write (1,'(e16.8)') geo_chi2(:) - write (1,'(e16.8)') -(geo_gcos1(:) + geo_gcos2(:))/rmaj - close(1) + !print *, geo_f, geo_ffprime/geo_f + !open(unit=1,file='out.cgyro.geogk',status='replace') + !write (1,'(i2)') n_theta + !write (1,'(e16.8)') theta(:) + !write (1,'(e16.8)') geo_b(:) + !write (1,'(e16.8)') 1.0/(q*rmaj*geo_g_theta(:)) + !write (1,'(e16.8)') (geo_dbdt(:)/geo_b(:))/(q*rmaj*geo_g_theta(:)) + !write (1,'(e16.8)') -1.0/geo_b(:)*geo_grad_r(:)/rmaj*geo_gsin(:) + !write (1,'(e16.8)') -1.0/geo_b(:)*geo_gq(:)/rmaj& + ! *(geo_gcos1(:)+geo_gcos2(:)+geo_captheta(:)*geo_gsin(:)) + !write (1,'(e16.8)') 1.0/geo_b(:)*geo_gq(:)/rmaj*geo_gcos2(:) + !write (1,'(e16.8)') geo_grad_r(:)**2 + !write (1,'(e16.8)') geo_gq(:)**2 * (1.0+geo_captheta(:)**2) + !write (1,'(e16.8)') 2.0*geo_grad_r(:)*geo_gq(:)*geo_captheta(:) + !write (1,'(e16.8)') geo_chi2(:) + !write (1,'(e16.8)') -(geo_gcos1(:) + geo_gcos2(:))/rmaj + !close(1) do it=1,n_theta