From 5473945172db3b80f8eadf6b87b7deab37a73664 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 8 May 2023 18:23:15 +0200
Subject: [PATCH] Npol is a real variable now (useful for 3D zpinch)
---
fort.90 | 7 ++++++-
src/geometry_mod.F90 | 8 +++++---
src/grid_mod.F90 | 48 +++++++++++++++++++++-----------------------
3 files changed, 34 insertions(+), 29 deletions(-)
diff --git a/fort.90 b/fort.90
index 88f8dbfb..3638ba08 100644
--- a/fort.90
+++ b/fort.90
@@ -19,8 +19,12 @@
&GEOMETRY
geom = 'Z-pinch'
q0 = 0.0
- shear = 0.0
+ shear = 0.0
eps = 0.0
+ !geom = 'miller'
+ !q0 = 1.4
+ !shear = 0.8
+ !eps = 0.18
kappa = 1.0
s_kappa= 0.0
delta = 0.0
@@ -29,6 +33,7 @@
s_zeta = 0.0
parallel_bc = 'shearless'
shift_y= 0.0
+ Npol = 1.0
/
&OUTPUT_PAR
dtsave_0d = 1
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index ffcedfdc..552a48c9 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -20,7 +20,7 @@ implicit none
REAL(xp), PUBLIC, PROTECTED :: s_zeta = 0._xp ! '' szeta = r dzeta/dr
! to apply shift in the parallel z-BC if shearless
REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3
- INTEGER, PUBLIC, PROTECTED :: Npol = 1 ! number of poloidal turns
+ REAL(xp), PUBLIC, PROTECTED :: Npol = 1._xp ! number of poloidal turns (real for 3D zpinch studies)
! Chooses the type of parallel BC we use for the unconnected kx modes (active for non-zero shear only)
! 'periodic' : Connect a disconnected kx to a mode on the other cadran
! 'dirichlet' : Connect a disconnected kx to 0
@@ -103,7 +103,7 @@ CONTAINS
REAL(xp) :: kx,ky
COMPLEX(xp), DIMENSION(local_nz) :: integrant
real(xp) :: G1,G2,G3,Cx,Cy
- INTEGER :: eo,iz,iky,ikx
+ INTEGER :: eo,iz,iky,ikx, Npol_i
! Allocate arrays
CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
@@ -114,6 +114,7 @@ CONTAINS
ELSE
SELECT CASE(geom)
CASE('s-alpha')
+ IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for s-alpha geometry"
CALL speak('s-alpha geometry')
call eval_salpha_geometry
CASE('Z-pinch','z-pinch','Zpinch','zpinch')
@@ -126,8 +127,9 @@ CONTAINS
kappa = 1._xp
CASE('miller')
CALL speak('Miller geometry')
+ IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
- call get_miller(eps,major_R,major_Z,q0,shear,Npol,alpha_MHD,edge_opt,&
+ call get_miller(eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
C_y,C_xy,xpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
dBdx,dBdy,hatB,jacobian,dBdz,hatR,hatZ,dxdR,dxdZ,&
Ckxky,gradz_coeff)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 4817c7b8..89f53324 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -146,8 +146,7 @@ CONTAINS
! 1 2 3 4
SUBROUTINE set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
USE fourier, ONLY: init_grid_distr_and_plans
- REAL(xp), INTENT(IN) :: shear
- INTEGER, INTENT(IN) :: Npol
+ REAL(xp), INTENT(IN) :: shear, Npol
CHARACTER(len=*), INTENT(IN) :: LINEARITY
INTEGER, INTENT(IN) :: N_HD
LOGICAL, INTENT(IN) :: EM
@@ -398,8 +397,7 @@ CONTAINS
SUBROUTINE set_kxgrid(shear,Npol,LINEARITY,N_HD)
USE prec_const
IMPLICIT NONE
- REAL(xp), INTENT(IN) :: shear
- INTEGER, INTENT(IN) :: Npol
+ REAL(xp), INTENT(IN) :: shear, Npol
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx
@@ -490,11 +488,11 @@ CONTAINS
USE prec_const
USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE
- REAL(xp):: grid_shift, Lz, zmax, zmin
- INTEGER :: istart, iend, in, Npol, iz, ig, eo, iglob
+ REAL(xp):: grid_shift, Lz, zmax, zmin, Npol
+ INTEGER :: istart, iend, in, iz, ig, eo, iglob
total_nz = Nz
! Length of the flux tube (in ballooning angle)
- Lz = 2._xp*pi*REAL(Npol,xp)
+ Lz = 2._xp*pi*Npol
! Z stepping (#interval = #points since periodic)
deltaz = Lz/REAL(Nz,xp)
inv_deltaz = 1._xp/deltaz
@@ -515,24 +513,24 @@ CONTAINS
ENDIF
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(zarray_full(total_nz))
- IF (Nz .EQ. 1) Npol = 0
- zmax = 0; zmin = 0;
- DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol
- zarray_full(iz) = REAL(iz-1,xp)*deltaz - Lz/2._xp
- IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz)
- IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
- END DO
- !! Parallel data distribution
- IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) &
- ERROR STOP '>> ERROR << Cannot have multiple core in z-direction (Nz = 1)'
- ! Local data distribution
- CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
- local_nz = ize - izs + 1
- local_nz_offset = izs - 1
- ! Ghosts boundaries (depend on the order of z operators)
- IF(Nz .EQ. 1) THEN
- ngz = 0
- zarray_full(izs) = 0
+ IF (Nz .EQ. 1) Npol = 0._xp
+ zmax = 0; zmin = 0;
+ DO iz = 1,total_nz ! z in [-pi pi-dz] x Npol
+ zarray_full(iz) = REAL(iz-1,xp)*deltaz - Lz/2._xp
+ IF(zarray_full(iz) .GT. zmax) zmax = zarray_full(iz)
+ IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
+ END DO
+ !! Parallel data distribution
+ IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) &
+ ERROR STOP '>> ERROR << Cannot have multiple core in z-direction (Nz = 1)'
+ ! Local data distribution
+ CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
+ local_nz = ize - izs + 1
+ local_nz_offset = izs - 1
+ ! Ghosts boundaries (depend on the order of z operators)
+ IF(Nz .EQ. 1) THEN
+ ngz = 0
+ zarray_full(izs) = 0
ELSEIF(Nz .GE. 4) THEN
ngz =4
IF(mod(Nz,2) .NE. 0 ) THEN
--
GitLab