diff --git a/README.md b/README.md
index d9eaf4429e5dfb7dafcd210eb8c9e8547f7b5b77..a7d1f316db418de222be5c2376c678ff82160125 100644
--- a/README.md
+++ b/README.md
@@ -8,29 +8,48 @@
GYACOMO (Gyrokinetic Advanced Collision Moment solver, 2021)
Copyright (C) 2022 EPFL
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
+This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-GNU General Public License for more details.
+This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-You should have received a copy of the GNU General Public License
-along with this program. If not, see <https://www.gnu.org/licenses/>.
+You should have received a copy of the GNU General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
Author: Antoine C.D. Hoffmann
# Citing GYACOMO
-If you use GYACOMO in your work, please cite the following paper:
+If you use GYACOMO in your work, please cite (at least) one of the following paper:
Hoffmann, A., Frei, B., & Ricci, P. (2023). Gyrokinetic simulations of plasma turbulence in a Z-pinch using a moment-based approach and advanced collision operators. Journal of Plasma Physics, 89(2), 905890214. doi:10.1017/S0022377823000284
+Hoffmann, A., Frei, B., & Ricci, P. (2023). Gyrokinetic moment-based simulations of the
+Dimits shift https://arxiv.org/abs/2308.01016
+
+# What is GYACOMO ?
+
+GYACOMO is the Gyrokinetic Advanced Collision Moment solver which solves the gyrokinetic Boltzmann equation in the delta-f flux-tube limit based on a projection of the velocity distribution function onto a Hermite-Laguerre velocity basis.
+
+It can be coupled with precomputed matrices from the code Cosolver (B.J. Frei) to incorporate advanced collision operators up to the gyro-averaged linearized exact coulomb interaction (GK Landau operator).
+
+This repository contains the solver source code (in /src) but also my personnal post-processing Matlab scripts, which are less documented. I would recommend the user to write their own post-processing scripts based on the H5 files the code outputs.
+
+#### GYACOMO can
+- evolve kinetic electrons and ions
+- use an adiabatic electrons model
+- include perpendicular magnetic perturbation (Ampere's law)
+- use Z-pinch and s-alpha geometry
+- use the Miller geometry framework (elongation, triangularity etc.)
+- use various experimental closures for the linear and nonlinear terms
+- add background ExB shear flow
+
+#### GYACOMO cannot (yet)
+- include parallel magnetic fluctuations
+- use an adiabatic ion model
+- include finite rhostar effects
+- study stellarator geometries
+
# How to compile and run GYACOMO
-A tutorial is present on the code's wiki https://gitlab.epfl.ch/ahoffman/gyacomo/-/wikis/home.
+A tutorial is now on the code's wiki https://gitlab.epfl.ch/ahoffman/gyacomo/-/wikis/home.
(older guideline)
To compile and run GYACOMO, follow these steps:
@@ -51,6 +70,8 @@ Note: For some collision operators (Sugama and Full Coulomb), you will need to r
### 4.x GYACOMO
+>4.11 background ExB shear is implemented
+
>4.1 Miller geometry is added and benchmarked for CBC adiabatic electrons
>4.0 new name and opening the code with GNU GPLv3 license
diff --git a/src/initial_mod.F90 b/src/initial_mod.F90
index d06c7cde87e44761b61523e11f198de0e7959866..a42e5d50bf4e356e5d4e665e24fccee11e421e35 100644
--- a/src/initial_mod.F90
+++ b/src/initial_mod.F90
@@ -85,7 +85,6 @@ CONTAINS
CALL apply_closure_model
CALL update_ghosts_moments
CALL solve_EM_fields ! compute phi_0=phi(N_0)
- CALL update_ghosts_EM
! through initialization
ELSE
SELECT CASE (INIT_OPT)
@@ -93,18 +92,15 @@ CONTAINS
CASE ('phi')
CALL speak('Init noisy phi')
CALL init_phi
- CALL update_ghosts_EM
CASE ('phi_ppj')
CALL speak('Init noisy phi')
CALL init_phi_ppj
- CALL update_ghosts_EM
! set moments_00 (GC density) with noise and compute phi afterwards
CASE('mom00')
CALL speak('Init noisy gyrocenter density')
CALL init_gyrodens ! init only gyrocenter density
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
! set moments_00 (GC density) with a single kx0,ky0 mode
! kx0 is setup but by the noise value and ky0 by the background value
CASE('mom00_modes','mom00_mode')
@@ -112,37 +108,32 @@ CONTAINS
CALL init_modes ! init single mode in gyrocenter density
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
! init all moments randomly (unadvised)
CASE('allmom')
CALL speak('Init noisy moments')
CALL init_moments ! init all moments
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
! init a gaussian blob in gyrodens
CASE('blob')
CALL speak('--init a blob')
CALL initialize_blob
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
! init moments 00 with a power law similarly to GENE
CASE('ppj')
CALL speak('ppj init ~ GENE')
call init_ppj
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
CASE('ricci')
CALL speak('Init Ricci')
CALL init_ricci ! init only gyrocenter density
CALL update_ghosts_moments
CALL solve_EM_fields
- CALL update_ghosts_EM
CASE DEFAULT
ERROR STOP "Initialization mode not recognized"
- END SELECT
+ END SELECT
ENDIF
! closure of j>J, p>P and j<0, p<0 moments
CALL speak('Apply closure')
@@ -460,7 +451,7 @@ CONTAINS
!******************************************************************************!
SUBROUTINE initialize_blob
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz, total_nz,&
- AA_x, AA_y, parray, jarray,&
+ AA_x, AA_y, parray, jarray, two_third_kxmax, two_third_kymax,&
ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
USE fields, ONLY: moments
USE prec_const, ONLY: xp
@@ -469,13 +460,11 @@ CONTAINS
IMPLICIT NONE
REAL(xp) ::kx, ky, z, sigma_x, sigma_y, gain
INTEGER :: ia,iky,ikx,iz,ip,ij, p, j
- sigma_y = 0.5
- sigma_x = sigma_y
+ sigma_y = two_third_kymax/2._xp
+ sigma_x = two_third_kxmax/2._xp
gain = 1.0
! One can increase the gain if we run 3D sim
IF(total_nz .GT. 1) THEN
- sigma_y = 1.0
- sigma_x = sigma_y
gain = 10.0
ENDIF
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 88f561f9c7d32fca188af440197bbeda2c581f2d..579ca824b6a9259327d46a5f1d4c3d0aa7e865de 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -133,8 +133,8 @@ CONTAINS
zjp1_cp = z_cp_stretched(z_idx_mapping(izg)+1)
zerotoone = (zi - zj_cp)/(zjp1_cp-zj_cp) ! weight for interpolation
moments(ia,ipi,iji,iky,ikx,izi,:) = &
- zerotoone *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg))&
- +(1._xp-zerotoone)*moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg+1))
+ (1._xp-zerotoone) *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg))&
+ + zerotoone *moments_cp(iacp,ipcp,ijcp,iycp,ixcp,z_idx_mapping(izg+1))
ELSE ! copy on all planes
moments(ia,ipi,iji,iky,ikx,izi,:) = moments_cp(iacp,ipcp,ijcp,iycp,ixcp,1)
ENDIF
@@ -177,35 +177,45 @@ CONTAINS
REAL(xp), DIMENSION(Nz_cp+1) , INTENT(OUT) :: z_cp_stretched
! local variables
INTEGER :: iz_new, jz_cp
- REAL(xp):: zi, zj_cp, dz_s
+ REAL(xp):: zi, zj_cp, dz_s, dz_new
LOGICAL :: in_interval
! stretched checkpoint grid interval
dz_s = (zarray_new(Nz_new)-zarray_new(1))/(Nz_cp-1)
- ! We loop over each new z grid points
- DO iz_new = 1,Nz_new
- zi = zarray_new(iz_new) ! current position
- ! Loop over the stretched checkpoint grid to find the right interval
- zj_cp = zarray_new(1) ! init cp grid position
- jz_cp = 1 ! init cp grid index
- in_interval = .FALSE. ! flag to check if we stand in the interval
- DO WHILE (.NOT. in_interval)
- in_interval = (zi .GE. zj_cp) .AND. (zi .LT. zj_cp+dz_s)
- ! Increment
- zj_cp = zj_cp + dz_s
- jz_cp = jz_cp + 1
- IF(jz_cp .GT. Nz_cp+1) ERROR STOP "STOP: could not adapt grid .."
- ENDDO ! per construction the while loop should always top
- z_idx_mapping(iz_new) = jz_cp-1 ! The last index was one too much so we store the one before
- ENDDO
- ! we build explicitly the stretched cp grid for output and double check
- DO jz_cp = 1,Nz_cp
- z_cp_stretched(jz_cp) = zarray_new(1) + (jz_cp-1)*dz_s
- ENDDO
- ! Periodicity
- z_cp_stretched(Nz_cp+1) = z_cp_stretched(1)
- z_idx_mapping (Nz_new+1) = z_idx_mapping (1)
- ! Check that the start and the end of the grids are the same
- IF(.NOT.(abs(z_cp_stretched(1)-zarray_new(1)) .LT. 1e-3).AND.(abs(z_cp_stretched(Nz_cp)-zarray_new(Nz_new)).LT.1e-3)) &
- ERROR STOP "Failed to stretch the cp grid"
+ dz_new= (zarray_new(Nz_new)-zarray_new(1))/(Nz_new-1)
+ IF(ABS(dz_s-dz_new) .LT. EPSILON(dz_s)) THEN
+ DO iz_new = 1,Nz_new
+ z_cp_stretched(iz_new) = zarray_new(iz_new)
+ z_idx_mapping(iz_new) = iz_new
+ ENDDO
+ z_cp_stretched(Nz_new+1) = zarray_new(1)
+ z_idx_mapping(Nz_new+1) = 1
+ ELSE ! Strech the grid
+ ! We loop over each new z grid points
+ DO iz_new = 1,Nz_new
+ zi = zarray_new(iz_new) ! current position
+ ! Loop over the stretched checkpoint grid to find the right interval
+ zj_cp = zarray_new(1) ! init cp grid position
+ jz_cp = 1 ! init cp grid index
+ in_interval = .FALSE. ! flag to check if we stand in the interval
+ DO WHILE (.NOT. in_interval)
+ in_interval = (zi .GE. zj_cp) .AND. (zi .LT. zj_cp+dz_s)
+ ! Increment
+ zj_cp = zj_cp + dz_s
+ jz_cp = jz_cp + 1
+ IF(jz_cp .GT. Nz_cp+1) ERROR STOP "STOP: could not adapt grid .."
+ ENDDO ! per construction the while loop should always top
+ z_idx_mapping(iz_new) = jz_cp-1 ! The last index was one too much so we store the one before
+ ENDDO
+ ! we build explicitly the stretched cp grid for output and double check
+ DO jz_cp = 1,Nz_cp
+ z_cp_stretched(jz_cp) = zarray_new(1) + (jz_cp-1)*dz_s
+ ENDDO
+ ! Periodicity
+ z_cp_stretched(Nz_cp+1) = z_cp_stretched(1)
+ z_idx_mapping (Nz_new+1) = z_idx_mapping(1)
+ ! Check that the start and the end of the grids are the same
+ IF(.NOT.(abs(z_cp_stretched(1)-zarray_new(1)) .LT. 1e-3).AND.(abs(z_cp_stretched(Nz_cp)-zarray_new(Nz_new)).LT.1e-3)) &
+ ERROR STOP "Failed to stretch the cp grid"
+ ENDIF
END SUBROUTINE z_grid_mapping
END MODULE restarts