From 789906ab465e5a49069cb96f5fd0817acdccf106 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 25 Apr 2022 13:55:31 +0200
Subject: [PATCH] created a geometry namelist for input and moved variables to
geom module
---
src/array_mod.F90 | 17 +---
src/diagnose.F90 | 7 +-
src/geometry_mod.F90 | 157 +++++++++++++++++++++++++++++++++++--
src/grid_mod.F90 | 29 ++++---
src/memory.F90 | 23 ------
src/moments_eq_rhs_mod.F90 | 1 +
src/processing_mod.F90 | 9 ++-
src/readinputs.F90 | 4 +
8 files changed, 183 insertions(+), 64 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 841f80d2..5ac91e6e 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -48,22 +48,7 @@ MODULE array
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1 ! mirror lin coeff for adiab mom
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_e, xphijp1_e, xphijm1_e
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_i, xphijp1_i, xphijm1_i
- ! Geometrical operators
- ! Curvature
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
- ! Jacobian
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
- ! Metric
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
- ! derivatives of magnetic field strength
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
- ! Relative magnetic field strength
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: hatB
- ! Relative strength of major radius
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
- ! Some geometrical coefficients
- REAL(dp), DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+
! Kernel function evaluation (ij,ikx,iky,iz,odd/even p)
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index f3efd5ff..50c0e171 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -13,6 +13,7 @@ SUBROUTINE diagnose(kstep)
USE utility
USE prec_const
USE collision, ONLY: coll_outputinputs
+ USE geometry
IMPLICIT NONE
INCLUDE 'srcinfo.h'
@@ -227,15 +228,11 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, TRIM(str), "write_temp", write_temp)
CALL grid_outputinputs(fidres, str)
-
+ CALL geometry_outputinputs(fidres, str)
CALL diag_par_outputinputs(fidres, str)
-
CALL model_outputinputs(fidres, str)
-
CALL coll_outputinputs(fidres, str)
-
CALL initial_outputinputs(fidres, str)
-
CALL time_integration_outputinputs(fidres, str)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index ebf1a642..9c8ff809 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -11,10 +11,45 @@ module geometry
use calculus, ONLY: simpson_rule_z
implicit none
- public
- COMPLEX(dp), PROTECTED :: iInt_Jacobian
+ PRIVATE
+ ! Geometry parameters
+ CHARACTER(len=16), &
+ PUBLIC, PROTECTED :: geom
+ REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
+ REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
+ REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
-contains
+ ! Geometrical operators
+ ! Curvature
+ REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
+ ! Jacobian
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
+ COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
+ ! Metric
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
+ ! derivatives of magnetic field strength
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
+ ! Relative magnetic field strength
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
+ ! Relative strength of major radius
+ REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
+ ! Some geometrical coefficients
+ REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+
+ ! Functions
+ PUBLIC :: geometry_readinputs, geometry_outputinputs, eval_magnetic_geometry
+
+CONTAINS
+
+
+ SUBROUTINE geometry_readinputs
+ ! Read the input parameters
+ IMPLICIT NONE
+ NAMELIST /GEOMETRY/ geom, q0, shear, eps
+ READ(lu_in,geometry)
+
+ END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
@@ -22,13 +57,24 @@ contains
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
INTEGER :: fid
+
+ ! Allocate arrays
+ CALL geometry_allocate_mem
!
IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
call eval_1D_geometry
ELSE
- IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
- call eval_salphaB_geometry
+ SELECT CASE(geom)
+ CASE('s-alpha')
+ IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
+ call eval_salphaB_geometry
+ CASE('Z-pinch')
+ IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
+ call eval_zpinch_geometry
+ CASE DEFAULT
+ ERROR STOP 'Error stop: geometry not recognized!!'
+ END SELECT
ENDIF
!
! Evaluate perpendicular wavenumber
@@ -118,6 +164,64 @@ contains
!--------------------------------------------------------------------------------
!
+ SUBROUTINE eval_zpinch_geometry
+ ! evaluate s-alpha geometry model
+ implicit none
+ REAL(dp) :: z, kx, ky, alpha_MHD
+ alpha_MHD = 0._dp
+
+ parity: DO eo = 0,1
+ zloop: DO iz = izgs,izge
+ z = zarray(iz,eo)
+
+ ! metric
+ gxx(iz,eo) = 1._dp
+ gxy(iz,eo) = 0._dp
+ gxz(iz,eo) = 0._dp
+ gyy(iz,eo) = 1._dp
+ gyz(iz,eo) = 0._dp
+ gzz(iz,eo) = 1._dp
+ dxdR(iz,eo)= COS(z)
+ dxdZ(iz,eo)= SIN(z)
+
+ ! Relative strengh of radius
+ hatR(iz,eo) = 1._dp
+ hatZ(iz,eo) = 1._dp
+
+ ! toroidal coordinates
+ Rc (iz,eo) = hatR(iz,eo)
+ phic(iz,eo) = z
+ Zc (iz,eo) = hatZ(iz,eo)
+
+ ! Jacobian
+ Jacobian(iz,eo) = 1._dp
+
+ ! Relative strengh of modulus of B
+ hatB(iz,eo) = 1._dp
+
+ ! Derivative of the magnetic field strenght
+ gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
+
+ ! Curvature operator
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky(ikx, iky, iz,eo) = - ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
+
+ ENDDO zloop
+ ENDDO parity
+
+ END SUBROUTINE eval_zpinch_geometry
+ !
+ !--------------------------------------------------------------------------------
+ !
subroutine eval_1D_geometry
! evaluate 1D perp geometry model
implicit none
@@ -158,4 +262,47 @@ contains
!
!--------------------------------------------------------------------------------
!
+
+ SUBROUTINE geometry_allocate_mem
+
+ ! Curvature and geometry
+ CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+ CALL allocate_array( Jacobian,izgs,izge, 0,1)
+ CALL allocate_array( gxx,izgs,izge, 0,1)
+ CALL allocate_array( gxy,izgs,izge, 0,1)
+ CALL allocate_array( gxz,izgs,izge, 0,1)
+ CALL allocate_array( gyy,izgs,izge, 0,1)
+ CALL allocate_array( gyz,izgs,izge, 0,1)
+ CALL allocate_array( gzz,izgs,izge, 0,1)
+ CALL allocate_array( gradxB,izgs,izge, 0,1)
+ CALL allocate_array( gradyB,izgs,izge, 0,1)
+ CALL allocate_array( gradzB,izgs,izge, 0,1)
+ CALL allocate_array( hatB,izgs,izge, 0,1)
+ CALL allocate_array( hatR,izgs,izge, 0,1)
+ CALL allocate_array( hatZ,izgs,izge, 0,1)
+ CALL allocate_array( Rc,izgs,izge, 0,1)
+ CALL allocate_array( phic,izgs,izge, 0,1)
+ CALL allocate_array( Zc,izgs,izge, 0,1)
+ CALL allocate_array( dxdR,izgs,izge, 0,1)
+ CALL allocate_array( dxdZ,izgs,izge, 0,1)
+ call allocate_array(gradz_coeff,izgs,izge, 0,1)
+ CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+
+ END SUBROUTINE geometry_allocate_mem
+
+ SUBROUTINE geometry_outputinputs(fidres, str)
+ ! Write the input parameters to the results_xx.h5 file
+
+ USE futils, ONLY: attach
+
+ USE prec_const
+ IMPLICIT NONE
+
+ INTEGER, INTENT(in) :: fidres
+ CHARACTER(len=256), INTENT(in) :: str
+ CALL attach(fidres, TRIM(str),"geometry", geom)
+ CALL attach(fidres, TRIM(str), "q0", q0)
+ CALL attach(fidres, TRIM(str), "shear", shear)
+ CALL attach(fidres, TRIM(str), "eps", eps)
+ END SUBROUTINE geometry_outputinputs
end module geometry
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 782699c6..52889c10 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -21,9 +21,6 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: Nz = 1 ! Number of total perpendicular planes
REAL(dp), PUBLIC, PROTECTED :: Npol = 1._dp ! number of poloidal turns
INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes
- REAL(dp), PUBLIC, PROTECTED :: q0 = 1._dp ! safety factor
- REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
- REAL(dp), PUBLIC, PROTECTED :: eps = 0._dp ! inverse aspect ratio
INTEGER, PUBLIC, PROTECTED :: Nkx = 8 ! Number of total internal grid points in kx
REAL(dp), PUBLIC, PROTECTED :: Lkx = 1._dp ! horizontal length of the fourier box
INTEGER, PUBLIC, PROTECTED :: Nky = 16 ! Number of total internal grid points in ky
@@ -120,7 +117,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nx, Lx, Ny, Ly, Nz, Npol, q0, shear, eps, SG
+ Nx, Lx, Ny, Ly, Nz, Npol, SG
READ(lu_in,grid)
!! Compute the maximal degree of full GF moments set
@@ -136,6 +133,7 @@ CONTAINS
IF(Nz .EQ. 1) THEN
deltape = 2; deltapi = 2;
pp2 = 1; ! index p+2 is ip+1
+ SG = .FALSE.
ELSE
deltape = 1; deltapi = 1;
pp2 = 2; ! index p+2 is ip+1
@@ -336,13 +334,13 @@ CONTAINS
INTEGER :: i_, counter
Nky = Ny;
- ALLOCATE(kyarray_full(1:Nky))
! Local data
! Start and END indices of grid
ikys = 1
ikye = Nky
local_nky = ikye - ikys + 1
ALLOCATE(kyarray(ikys:ikye))
+ ALLOCATE(kyarray_full(1:Nky))
IF (Ny .EQ. 1) THEN ! "cancel" y dimension
deltaky = 1._dp
kyarray(1) = 0._dp
@@ -360,7 +358,12 @@ CONTAINS
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
local_kymax = 0._dp
DO iky = ikys,ikye
- kyarray(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
+ SELECT CASE (LINEARITY)
+ CASE ('linear') ! only positive freq for linear runs dk*(0 1 2 3 4 5)
+ kyarray(iky) = deltaky*REAL(iky-1,dp)
+ CASE DEFAULT !
+ kyarray(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
+ END SELECT
if (iky .EQ. Ny/2+1) kyarray(iky) = -kyarray(iky)
! Finding ky=0
IF (kyarray(iky) .EQ. 0) THEN
@@ -376,9 +379,14 @@ CONTAINS
END DO
! Build the full grids on process 0 to diagnose it without comm
! ky
- DO iky = ikys,ikye
- kyarray_full(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
- IF (iky .EQ. Ny/2+1) kyarray_full(iky) = -kyarray_full(iky)
+ DO iky = 1,Nky
+ SELECT CASE (LINEARITY)
+ CASE ('linear') ! only positive freq for linear runs dk*(0 1 2 3 4 5)
+ kyarray_full(iky) = deltaky*REAL(iky-1,dp)
+ CASE DEFAULT !
+ kyarray_full(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
+ IF (iky .EQ. Ny/2+1) kyarray_full(iky) = -kyarray_full(iky)
+ END SELECT
END DO
ENDIF
! Orszag 2/3 filter
@@ -491,9 +499,6 @@ CONTAINS
CALL attach(fidres, TRIM(str), "Ny", Ny)
CALL attach(fidres, TRIM(str), "Ly", Ly)
CALL attach(fidres, TRIM(str), "Nz", Nz)
- CALL attach(fidres, TRIM(str), "q0", q0)
- CALL attach(fidres, TRIM(str),"shear",shear)
- CALL attach(fidres, TRIM(str), "eps", eps)
CALL attach(fidres, TRIM(str), "Nkx", Nkx)
CALL attach(fidres, TRIM(str), "Lkx", Lkx)
CALL attach(fidres, TRIM(str), "Nky", Nky)
diff --git a/src/memory.F90 b/src/memory.F90
index 5c94ee7e..3f51d7aa 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -95,29 +95,6 @@ SUBROUTINE memory
CALL allocate_array( xphijp1_i, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijm1_i, ips_i,ipe_i, ijs_i,ije_i)
- ! Curvature and geometry
- CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
- CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
- CALL allocate_array( Jacobian,izgs,izge, 0,1)
- CALL allocate_array( gxx,izgs,izge, 0,1)
- CALL allocate_array( gxy,izgs,izge, 0,1)
- CALL allocate_array( gxz,izgs,izge, 0,1)
- CALL allocate_array( gyy,izgs,izge, 0,1)
- CALL allocate_array( gyz,izgs,izge, 0,1)
- CALL allocate_array( gzz,izgs,izge, 0,1)
- CALL allocate_array( gradxB,izgs,izge, 0,1)
- CALL allocate_array( gradyB,izgs,izge, 0,1)
- CALL allocate_array( gradzB,izgs,izge, 0,1)
- CALL allocate_array( hatB,izgs,izge, 0,1)
- CALL allocate_array( hatR,izgs,izge, 0,1)
- CALL allocate_array( hatZ,izgs,izge, 0,1)
- CALL allocate_array( Rc,izgs,izge, 0,1)
- CALL allocate_array( phic,izgs,izge, 0,1)
- CALL allocate_array( Zc,izgs,izge, 0,1)
- CALL allocate_array( dxdR,izgs,izge, 0,1)
- CALL allocate_array( dxdZ,izgs,izge, 0,1)
- call allocate_array(gradz_coeff,izgs,izge, 0,1)
-
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
IF (gyrokin_CO) THEN !GK collision matrices (one for each kperp)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 0934b533..ac081a88 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -152,6 +152,7 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
+ USE geometry
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 0a4878f4..7139fe5c 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -5,6 +5,7 @@ MODULE processing
USE grid
USE geometry
USE utility
+ USE calculus
implicit none
REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
@@ -19,7 +20,8 @@ CONTAINS
! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
SUBROUTINE compute_radial_ion_transport
USE fields, ONLY : moments_i, phi
- USE array, ONLY : kernel_i, Jacobian
+ USE array, ONLY : kernel_i
+ USE geometry, ONLY : Jacobian
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
COMPLEX(dp) :: pflux_local, gflux_local, integral
@@ -92,7 +94,8 @@ END SUBROUTINE compute_radial_ion_transport
! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
SUBROUTINE compute_radial_heatflux
USE fields, ONLY : moments_i, moments_e, phi
- USE array, ONLY : kernel_e, kernel_i, Jacobian
+ USE array, ONLY : kernel_e, kernel_i
+ USE geometry, ONLY : Jacobian
USE time_integration, ONLY : updatetlevel
USE model, ONLY : qe_taue, qi_taui, KIN_E
IMPLICIT NONE
@@ -404,7 +407,7 @@ END SUBROUTINE compute_Tpar
! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre)
SUBROUTINE compute_fluid_moments
- USE array, ONLY : dens_e, Tpar_e, Tper_e, dens_i, Tpar_i, Tper_i
+ USE array, ONLY : dens_e, Tpar_e, Tper_e, dens_i, Tpar_i, Tper_i, temp_e, temp_i
USE model, ONLY : KIN_E
IMPLICIT NONE
CALL compute_density
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 1322ddeb..56cdca06 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -7,6 +7,7 @@ SUBROUTINE readinputs
USE model, ONLY: model_readinputs
USE initial_par, ONLY: initial_readinputs
USE time_integration, ONLY: time_integration_readinputs
+ USE geometry, ONLY: geometry_readinputs
USE prec_const
IMPLICIT NONE
@@ -18,6 +19,9 @@ SUBROUTINE readinputs
! Load grid data from input file
CALL grid_readinputs
+ ! Load geometry parameters
+ CALL geometry_readinputs
+
! Load diagnostic options from input file
CALL diag_par_readinputs
--
GitLab