Skip to content
Snippets Groups Projects
Commit f4c3b659 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

add circular geometry

parent cc60ce6c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Makefile
+ 14
− 9
View file @ f4c3b659
......@@ -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 \
......
src/circular_mod.F90 0 → 100644
+ 146
− 0
View file @ f4c3b659
!! 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment