From f4c3b659a2c2a918c9eb1b84b11381f9bc15d2ec Mon Sep 17 00:00:00 2001 From: Antoine Hoffmann <antoine.hoffmann@epfl.ch> Date: Tue, 27 Jun 2023 08:51:24 +0200 Subject: [PATCH] add circular geometry --- Makefile | 23 ++++--- src/circular_mod.F90 | 146 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 160 insertions(+), 9 deletions(-) create mode 100644 src/circular_mod.F90 diff --git a/Makefile b/Makefile index 51dc0420..48b38359 100644 --- a/Makefile +++ b/Makefile @@ -24,7 +24,7 @@ fast: compile # Debug version with all flags debug: F90 = mpif90 -debug: F90FLAGS = -C -g -traceback -ftrapuv -warn all -debug all +debug: F90FLAGS = -C -g -traceback -warn all -check all -debug all -init=zero,snan -fpe0 debug: EXEC = $(BINDIR)/gyacomo23_debug debug: dirs src/srcinfo.h debug: compile @@ -112,7 +112,7 @@ $(OBJDIR)/diagnose.o : src/srcinfo.h # Main source dependencies FOBJ=$(OBJDIR)/advance_field_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o \ -$(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o \ +$(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/circular_mod.o \ $(OBJDIR)/collision_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o \ $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \ $(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/geometry_mod.o \ @@ -155,6 +155,11 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o $(OBJDIR)/parallel_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/calculus_mod.F90 -o $@ + $(OBJDIR)/circular_mod.o : src/circular_mod.F90 \ + $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/lag_interp_mod.o\ + $(OBJDIR)/prec_const_mod.o + $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/circular_mod.F90 -o $@ + $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 \ $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o \ $(OBJDIR)/basic_mod.o @@ -207,14 +212,14 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 \ - $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/miller_mod.o \ - $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o \ - $(OBJDIR)/utility_mod.o + $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/circular_mod.o \ + $(OBJDIR)/miller_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o \ + $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@ $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 \ - $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o\ - $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o + $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o \ + $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@ $(OBJDIR)/grid_mod.o : src/grid_mod.F90 \ @@ -247,8 +252,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@ $(OBJDIR)/miller_mod.o : src/miller_mod.F90 \ - $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/lag_interp_mod.o \ - $(OBJDIR)/prec_const_mod.o + $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/lag_interp_mod.o\ + $(OBJDIR)/prec_const_mod.o $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/miller_mod.F90 -o $@ $(OBJDIR)/model_mod.o : src/model_mod.F90 \ diff --git a/src/circular_mod.F90 b/src/circular_mod.F90 new file mode 100644 index 00000000..cff23018 --- /dev/null +++ b/src/circular_mod.F90 @@ -0,0 +1,146 @@ +!! This source has been adapted from the GENE code https://genecode.org/ !! +!>Provides a circular concentric geometry +MODULE circular + USE prec_const + USE basic + USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0 + ! use coordinates,only: gcoor, get_dzprimedz + USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven + ! use discretization + USE lagrange_interpolation + USE model + + implicit none + public:: get_circ + private + +! SUBROUTINES DEFINITIONS +CONTAINS + + SUBROUTINE get_circ(trpeps,major_R,major_Z,q0,shat,& + C_y,C_xy,Cyq0_x0,gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,dBdx_,dBdy_,dBdz_,& + Bfield_,jacobian_,R_hat_,Z_hat_,dxdR_,dxdZ_) + !!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + real(xp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio) + real(xp), INTENT(INOUT) :: major_R ! major radius + real(xp), INTENT(INOUT) :: major_Z ! major Z + real(xp), INTENT(INOUT) :: q0 ! safetyfactor + real(xp), INTENT(INOUT) :: shat ! safetyfactor + real(xp), INTENT(INOUT) :: C_y, C_xy, Cyq0_x0 + real(xp), dimension(local_nz+ngz,nzgrid), INTENT(OUT) :: & + gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,& + dBdx_,dBdy_,dBdz_,Bfield_,jacobian_,& + R_hat_,Z_hat_,dxdR_,dxdZ_ + INTEGER :: iz,eo + ! No parameter in gyacomo yet + real(xp) :: SIGN_Ip=1._xp ! current sign (only normal current) + !!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + REAL(xp) :: g11,g12,g22,g33 ! normalized metric in (PSI,CHI,PHI) + REAL(xp) :: B,dBdr_c,dBdchi ! J_PCP : J_PSI,CHI,P + ! intermediate variables, t stands for theta the poloidal angle. + REAL(xp) :: qbar,dpsidr,dchidr,dchidt,dBdr,dBdt,cost,sint + REAL(xp) :: fac2, dCydr_Cy + REAL(xp) :: dxdr, dqdr, dqdx + REAL(xp) :: dqbardr, minor_r, r0, z + + minor_r = major_R + r0 = major_R*trpeps + qbar = abs(q0)*sqrt(1-trpeps**2) + C_y = abs(r0/q0)*SIGN_Ip + dxdr = 1._xp ! x=r + dpsidr = r0/qbar ! 1/(Bref*Lref)*dpsi/dr + ! dqdx = shat*q0/minor_r/dxdr ! AH: added + dqdx = shat/C_y ! AH: added + dqdr = dqdx*dxdr/minor_r + dqbardr = sqrt(1._xp-trpeps**2)*(abs(dqdr)-abs(q0)*trpeps/major_R/(1-trpeps**2)) + dCydr_Cy = 0._xp + C_xy = dpsidr/(dxdr*abs(C_y)) + Cyq0_x0 = C_y*q0/r0 + fac2 = (dCydr_Cy+dqdr/q0) + + DO eo = 1,nzgrid + do iz=1,local_nz+ngz + z = zarray(iz,eo) + ! z dep. intermediate variables + cost = (cos(z) - trpeps)/(1._xp-trpeps*cos(z)) + sint = sqrt(1._xp-trpeps**2)*sin(SIGN_Ip*z)/(1._xp-trpeps*cos(z)) + ! derivatives with respect to r,theta + dchidr = -1._xp/major_R*sin(z)/(1._xp-trpeps**2) + dchidt = SIGN_Ip * sqrt(1-trpeps**2)/(1+trpeps*cost) + + B = 1._xp/(1._xp+trpeps*cost)*sqrt(1.+(trpeps/qbar)**2) + dBdr = B/major_R*(-cost/(1._xp+trpeps*cost) + & + + trpeps/(trpeps**2+qbar**2)*(1._xp - major_R*trpeps*dqbardr/qbar)) + dBdt = trpeps*sint*B / (1._xp+trpeps*cost) + + ! metric in (r,CHI,PHI) + g11 = 1._xp + g22 = dchidr**2+1._xp/(major_R*trpeps)**2*dchidt**2 + g12 = dchidr + g33 = 1._xp/major_R**2*1._xp/(1._xp+trpeps*cost)**2 + + ! magnetic field derivatives in + dBdchi = dBdt/dchidt + dBdr_c = 1._xp/g11*(dBdr-g12*dBdt/dchidt) + + ! metric in (x,y,z) + gxx_(iz,eo) = (dxdr)**2*g11 + gyy_(iz,eo) = (C_y*q0)**2 * ((fac2*z)**2*g11 +2._xp*fac2*z*g12+g22) + C_y**2*g33 + gxy_(iz,eo) = dxdr*C_y*SIGN_Ip*q0*(fac2*z*g11+g12) + gxz_(iz,eo) = dxdr*g12 + gyz_(iz,eo) = C_y*q0*SIGN_Ip*(fac2*z*g12+g22) + gzz_(iz,eo) = g22 + + ! jacobian + jacobian_(iz,eo)=C_xy*abs(q0)*major_R*(1+trpeps*cost)**2 + + ! Bfield + Bfield_(iz,eo) = B + dBdx_(iz,eo) = dBdr_c/dxdr + dBdy_(iz,eo) = 0._xp + dBdz_(iz,eo) = dBdchi + + !derivatives with respect to cylindrical coordinates + dxdR_(iz,eo) = cost + dxdZ_(iz,eo) = sint + R_hat_(iz,eo) = major_R*(1.+trpeps*cost) + Z_hat_(iz,eo) = major_Z+major_R*trpeps*sint + + !!!!!!! TEST salpha + ! ! metric + ! gxx_(iz,eo) = 1._xp + ! gxy_(iz,eo) = shat*z - amhd*SIN(z) + ! gxz_(iz,eo) = 0._xp + ! gyy_(iz,eo) = 1._xp + (shat*z - amhd*SIN(z))**2 + ! gyz_(iz,eo) = 1._xp/trpeps + ! gzz_(iz,eo) = 1._xp/trpeps**2 + ! dxdR_(iz,eo)= COS(z) + ! dxdZ_(iz,eo)= SIN(z) + + ! ! ! Poloidal plane coordinates + ! R_hat_(iz,eo) = 1._xp + trpeps*COS(z) + ! Z_hat_(iz,eo) = trpeps*SIN(z) + + ! ! Relative strengh of modulus of B + ! Bfield_(iz,eo) = 1._xp/(1._xp + trpeps*COS(z)) + + ! ! Jacobian + ! Jacobian_(iz,eo) = q0/Bfield_(iz,eo) + + ! ! Derivative of the magnetic field strenght + ! dBdx_(iz,eo) = -COS(z)*Bfield_(iz,eo)**2 ! LB = 1 + ! dBdy_(iz,eo) = 0._xp + ! dBdz_(iz,eo) = trpeps*SIN(z)*Bfield_(iz,eo)**2 + + ! ! Curvature factor + ! C_xy = 1._xp + !!!!! END TEST + end do + END DO + + END SUBROUTINE get_circ + + !----------------------------------------------------------- +END MODULE circular -- GitLab