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

release almost ready

parent f69bcc08
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/stepon.F90
+ 151
− 128
View file @ 2489ac29
SUBROUTINE stepon SUBROUTINE stepon
! Advance one time step, (num_step=4 for Runge Kutta 4 scheme) ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
USE advance_field_routine, ONLY: advance_time_level, advance_moments USE advance_field_routine, ONLY: advance_time_level, advance_moments
USE basic, ONLY: nlend USE basic, ONLY: nlend, chrono_advf, chrono_pois,&
USE closure, ONLY: apply_closure_model chrono_chck, chrono_clos, chrono_ghst, start_chrono, stop_chrono
USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM USE closure, ONLY: apply_closure_model
use mpi, ONLY: MPI_COMM_WORLD USE ghosts, ONLY: update_ghosts_moments, update_ghosts_EM
USE time_integration, ONLY: ntimelevel use mpi, ONLY: MPI_COMM_WORLD
USE prec_const, ONLY: dp USE time_integration, ONLY: ntimelevel
IMPLICIT NONE USE prec_const, ONLY: dp
IMPLICIT NONE
INTEGER :: num_step, ierr
LOGICAL :: mlend INTEGER :: num_step, ierr
LOGICAL :: mlend
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4 DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- BEFORE: All fields+ghosts are updated for step = n !----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields ! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n) ! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
CALL assemble_RHS CALL assemble_RHS
...@@ -25,134 +26,156 @@ SUBROUTINE stepon ...@@ -25,134 +26,156 @@ SUBROUTINE stepon
! Update moments with the hierarchy RHS (step by step) ! Update moments with the hierarchy RHS (step by step)
! N_n+1 = N_n + N_rhs(n) ! N_n+1 = N_n + N_rhs(n)
CALL advance_moments CALL start_chrono(chrono_advf)
CALL advance_moments
CALL stop_chrono(chrono_advf)
! Closure enforcement of moments ! Closure enforcement of moments
CALL apply_closure_model CALL start_chrono(chrono_clos)
CALL apply_closure_model
CALL stop_chrono(chrono_clos)
! Exchanges the ghosts values of N_n+1 ! Exchanges the ghosts values of N_n+1
CALL update_ghosts_moments CALL start_chrono(chrono_ghst)
CALL update_ghosts_moments
CALL stop_chrono(chrono_ghst)
! Update electrostatic potential phi_n = phi(N_n+1) and potential vect psi ! Update electrostatic potential phi_n = phi(N_n+1) and potential vect psi
CALL solve_EM_fields CALL start_chrono(chrono_pois)
CALL update_ghosts_EM CALL solve_EM_fields
CALL stop_chrono(chrono_pois)
! Update EM ghosts
CALL start_chrono(chrono_ghst)
CALL update_ghosts_EM
CALL stop_chrono(chrono_ghst)
!- Check before next step !- Check before next step
CALL checkfield_all() CALL start_chrono(chrono_chck)
CALL checkfield_all()
CALL stop_chrono(chrono_chck)
IF( nlend ) EXIT ! exit do loop IF( nlend ) EXIT ! exit do loop
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
!----- AFTER: All fields are updated for step = n+1 !----- AFTER: All fields are updated for step = n+1
END DO END DO
CONTAINS CONTAINS
!!!! Basic structure to simplify stepon !!!! Basic structure to simplify stepon
SUBROUTINE assemble_RHS SUBROUTINE assemble_RHS
USE moments_eq_rhs, ONLY: compute_moments_eq_rhs USE basic, ONLY: chrono_mrhs, chrono_sapj, chrono_coll, chrono_grad, chrono_napj, start_chrono, stop_chrono
USE collision, ONLY: compute_Capj USE moments_eq_rhs, ONLY: compute_moments_eq_rhs
USE nonlinear, ONLY: compute_Sapj USE collision, ONLY: compute_Capj
USE processing, ONLY: compute_nadiab_moments_z_gradients_and_interp USE nonlinear, ONLY: compute_Sapj
IMPLICIT NONE USE processing, ONLY: compute_nadiab_moments, compute_interp_z, compute_gradients_z
! compute auxiliary non adiabatic moments array, gradients and interp IMPLICIT NONE
CALL compute_nadiab_moments_z_gradients_and_interp ! compute auxiliary non adiabatic moments array
! compute nonlinear term ("if linear" is included inside) CALL start_chrono(chrono_napj)
CALL compute_Sapj CALL compute_nadiab_moments
! compute collision term ("if coll, if nu >0" is included inside) CALL stop_chrono(chrono_napj)
CALL compute_Capj
! compute the moments equation rhs ! compute gradients and interp
CALL compute_moments_eq_rhs CALL start_chrono(chrono_grad)
END SUBROUTINE assemble_RHS CALL compute_gradients_z
CALL compute_interp_z
SUBROUTINE checkfield_all ! Check all the fields for inf or nan CALL stop_chrono(chrono_grad)
USE utility,ONLY: is_nan, is_inf
USE basic, ONLY: t0_checkfield, t1_checkfield, tc_checkfield ! compute nonlinear term ("if linear" is included inside)
USE fields, ONLY: phi CALL start_chrono(chrono_sapj)
USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,& CALL compute_Sapj
ngp,ngj,ngz CALL stop_chrono(chrono_sapj)
USE MPI
USE time_integration, ONLY: updatetlevel ! compute collision term ("if coll, if nu >0" is included inside)
USE model, ONLY: LINEARITY CALL start_chrono(chrono_coll)
IMPLICIT NONE CALL compute_Capj
LOGICAL :: checkf_ CALL stop_chrono(chrono_coll)
INTEGER :: ia, ip, ij, iky, ikx, iz
REAL :: sum_ ! compute the moments equation rhs
! Execution time start CALL start_chrono(chrono_mrhs)
CALL cpu_time(t0_checkfield) CALL compute_moments_eq_rhs
CALL stop_chrono(chrono_mrhs)
IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule END SUBROUTINE assemble_RHS
IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0
SUBROUTINE checkfield_all ! Check all the fields for inf or nan
mlend=.FALSE. USE utility,ONLY: is_nan, is_inf
IF(.NOT.nlend) THEN USE fields, ONLY: phi
sum_ = SUM(ABS(phi)) USE MPI
checkf_ = (is_nan(sum_,'phi') .OR. is_inf(sum_,'phi')) USE model, ONLY: LINEARITY
mlend = (mlend .or. checkf_) IMPLICIT NONE
CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) LOGICAL :: checkf_
ENDIF REAL :: sum_
! Execution time end IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule
CALL cpu_time(t1_checkfield) IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0
tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield)
END SUBROUTINE checkfield_all mlend=.FALSE.
IF(.NOT.nlend) THEN
SUBROUTINE anti_aliasing sum_ = SUM(ABS(phi))
USE fields, ONLY: moments checkf_ = (is_nan(sum_,'phi') .OR. is_inf(sum_,'phi'))
USE time_integration, ONLY: updatetlevel mlend = (mlend .or. checkf_)
USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,& CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
ngp,ngj,ngz, AA_x, AA_y ENDIF
IMPLICIT NONE END SUBROUTINE checkfield_all
INTEGER :: ia, ip, ij, iky, ikx, iz
z: DO iz = 1,local_nz+ngz SUBROUTINE anti_aliasing
kx:DO ikx= 1,local_nkx USE fields, ONLY: moments
ky:DO iky=1,local_nky USE time_integration, ONLY: updatetlevel
j: DO ij=1,local_nj+ngj USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
p: DO ip=1,local_np+ngp ngp,ngj,ngz, AA_x, AA_y
a: DO ia=1,local_na IMPLICIT NONE
moments(ia,ip,ij,iky,ikx,iz,updatetlevel) =& INTEGER :: ia, ip, ij, iky, ikx, iz
AA_x(ikx)* AA_y(iky) * moments(ia,ip,ij,iky,ikx,iz,updatetlevel) z: DO iz = 1,local_nz+ngz
ENDDO a kx:DO ikx= 1,local_nkx
ENDDO p ky:DO iky=1,local_nky
ENDDO j j: DO ij=1,local_nj+ngj
ENDDO ky p: DO ip=1,local_np+ngp
ENDDO kx a: DO ia=1,local_na
ENDDO z moments(ia,ip,ij,iky,ikx,iz,updatetlevel) =&
END SUBROUTINE anti_aliasing AA_x(ikx)* AA_y(iky) * moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
ENDDO a
SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry ENDDO p
USE fields, ONLY: phi, psi, moments ENDDO j
USE time_integration, ONLY: updatetlevel ENDDO ky
USE grid, ONLY: local_na,local_np,local_nj,total_nkx,local_nz,& ENDDO kx
ngp,ngj,ngz, ikx0,iky0, contains_ky0 ENDDO z
IMPLICIT NONE END SUBROUTINE anti_aliasing
INTEGER :: ia, ip, ij, ikx, iz
IF ( contains_ky0 ) THEN SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
! moments USE fields, ONLY: phi, psi, moments
! z: DO iz = 1,local_nz+ngz USE time_integration, ONLY: updatetlevel
! j: DO ij=1,local_nj+ngj USE grid, ONLY: total_nkx, ikx0,iky0, contains_ky0
! p: DO ip=1,local_np+ngp IMPLICIT NONE
! a: DO ia=1,local_na INTEGER :: ikx
DO ikx=2,total_nkx/2 !symmetry at ky = 0 IF ( contains_ky0 ) THEN
moments(:,:,:,iky0,ikx,:,updatetlevel) = & ! moments
CONJG(moments(:,:,:,iky0,total_nkx+2-ikx,:,updatetlevel)) ! z: DO iz = 1,local_nz+ngz
END DO ! j: DO ij=1,local_nj+ngj
! must be real at origin and top right ! p: DO ip=1,local_np+ngp
moments(:,:,:, iky0,ikx0,:,updatetlevel) = & ! a: DO ia=1,local_na
REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),dp) DO ikx=2,total_nkx/2 !symmetry at ky = 0
! ENDDO a moments(:,:,:,iky0,ikx,:,updatetlevel) = &
! ENDDO p CONJG(moments(:,:,:,iky0,total_nkx+2-ikx,:,updatetlevel))
! ENDDO j END DO
! ENDDO z ! must be real at origin and top right
! Phi moments(:,:,:, iky0,ikx0,:,updatetlevel) = &
DO ikx=2,total_nkx/2 !symmetry at ky = 0 REAL(moments(:,:,:, iky0,ikx0,:,updatetlevel),dp)
! ENDDO a
! ENDDO p
! ENDDO j
! ENDDO z
! Phi
DO ikx=2,total_nkx/2 !symmetry at ky = 0
phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:) phi(iky0,ikx,:) = phi(iky0,total_nkx+2-ikx,:)
END DO END DO
! must be real at origin ! must be real at origin
phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp) phi(iky0,ikx0,:) = REAL(phi(iky0,ikx0,:),dp)
! Psi ! Psi
DO ikx=2,total_nkx/2 !symmetry at ky = 0 DO ikx=2,total_nkx/2 !symmetry at ky = 0
psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:) psi(iky0,ikx,:) = psi(iky0,total_nkx+2-ikx,:)
END DO END DO
! must be real at origin ! must be real at origin
psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp) psi(iky0,ikx0,:) = REAL(psi(iky0,ikx0,:),dp)
ENDIF ENDIF
END SUBROUTINE enforce_symmetry END SUBROUTINE enforce_symmetry
END SUBROUTINE stepon END SUBROUTINE stepon
src/tesend.F90
+ 1
− 1
View file @ 2489ac29
...@@ -44,7 +44,7 @@ SUBROUTINE tesend ...@@ -44,7 +44,7 @@ SUBROUTINE tesend
!________________________________________________________________________________ !________________________________________________________________________________
! 4. Test on run time ! 4. Test on run time
CALL cpu_time(tnow) CALL cpu_time(tnow)
mlend = (1.1*(tnow-start)) .GT. maxruntime mlend = (1.1*(tnow-chrono_runt%tstart)) .GT. maxruntime
CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
......
testcases/cyclone_example/fort.90 deleted 100644 → 0
+ 0
− 92
View file @ f69bcc08
&BASIC
nrun = 99999999
dt = 0.01
tmax = 25
maxruntime = 356400
job2load = -1
/
&GRID
pmax = 4
jmax = 2
Nx = 128
Lx = 120
Ny = 48
Ly = 120
Nz = 16
SG = .f.
Nexc = 0
/
&GEOMETRY
geom = 's-alpha'
q0 = 1.4
shear = 0.8
eps = 0.18
kappa = 1.0
s_kappa= 0.0
delta = 0.0
s_delta= 0.0
zeta = 0.0
s_zeta = 0.0
parallel_bc = 'dirichlet'
shift_y= 0.0
/
&OUTPUT_PAR
dtsave_0d = 0.1
dtsave_1d = -1
dtsave_2d = -1
dtsave_3d = 1
dtsave_5d = 10
write_doubleprecision = .f.
write_gamma = .t.
write_hf = .t.
write_phi = .t.
write_Na00 = .t.
write_Napj = .t.
write_dens = .t.
write_fvel = .t.
write_temp = .t.
/
&MODEL_PAR
! Collisionality
CLOS = 0
NL_CLOS = 0
LINEARITY = 'nonlinear'
Na = 1 ! number of species
mu_x = 1.0
mu_y = 1.0
N_HD = 4
mu_z = 2.0
HYP_V = 'hypcoll'
mu_p = 0.0
mu_j = 0.0
nu = 0.001
beta = 0.0
ADIAB_E = .t.
tau_e = 1.0
/
&SPECIES
! ions
name_ = 'ions'
tau_ = 1.0
sigma_= 1.0
q_ = 1.0
k_N_ = 2.22
k_T_ = 6.96
/
&COLLISION_PAR
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
GK_CO = .f.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 0.0
init_noiselvl = 0.005
iseed = 42
/
&TIME_INTEGRATION_PAR
numerical_scheme = 'RK4'
/
wk/analysis_gyacomo.m
+ 9
− 8
View file @ 2489ac29
%% UNCOMMENT FOR TUTORIAL %% UNCOMMENT FOR TUTORIAL
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
resdir = '../testcases/zpinch_example/'; %Name of the directory where the results are located resdir = '../testcases/cyclone_example/'; %Name of the directory where the results are located
PARTITION =''; PARTITION ='';
JOBNUMMIN = 00; JOBNUMMAX = 10; JOBNUMMIN = 03; JOBNUMMAX = 03;
%% %%
addpath(genpath([gyacomodir,'matlab'])) % ... add addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
...@@ -80,7 +80,7 @@ options.NAME = '\phi'; ...@@ -80,7 +80,7 @@ options.NAME = '\phi';
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
% options.NAME = 'n_i-n_e'; % options.NAME = 'n_i-n_e';
options.PLAN = 'xy'; options.PLAN = 'kxky';
% options.NAME = 'f_i'; % options.NAME = 'f_i';
% options.PLAN = 'sx'; % options.PLAN = 'sx';
options.COMP = 'avg'; options.COMP = 'avg';
...@@ -93,27 +93,28 @@ options.RESOLUTION = 256; ...@@ -93,27 +93,28 @@ options.RESOLUTION = 256;
create_film(data,options,'.gif') create_film(data,options,'.gif')
end end
if 0 if 1
%% fields snapshots %% fields snapshots
% Options % Options
options.INTERP = 0; options.INTERP = 0;
options.POLARPLOT = 0; options.POLARPLOT = 0;
options.AXISEQUAL = 0; options.AXISEQUAL = 0;
options.NORMALIZE = 0; options.NORMALIZE = 0;
options.NAME = '\phi'; % options.NAME = '\phi';
% options.NAME = '\psi'; % options.NAME = '\psi';
% options.NAME = '\omega_z'; % options.NAME = '\omega_z';
% options.NAME = 'T_i'; % options.NAME = 'T_i';
% options.NAME = 'n_i'; % options.NAME = 'n_i';
% options.NAME = '\phi^{NZ}'; % options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}'; options.NAME = 'N_i^{00}';
% options.NAME = 'N_i^{00}-N_e^{00}'; % options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}'; % options.NAME = 's_{Ex}';
% options.NAME = 'Q_x'; % options.NAME = 'Q_x';
% options.NAME = 'k^2n_e'; % options.NAME = 'k^2n_e';
options.PLAN = 'xy'; options.PLAN = 'kxky';
options.COMP = 'avg'; options.COMP = 'avg';
options.TIME = [15 30 50]; % options.TIME = [49 50 51];
options.TIME = data.Ts3D(51:55);
options.RESOLUTION = 256; options.RESOLUTION = 256;
data.a = data.EPS * 2e3; data.a = data.EPS * 2e3;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment