From 665d0ec3526344e4968698a62680ca4cdbc899a2 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 18 Oct 2022 09:22:06 +0200
Subject: [PATCH] Miller works in linear ITG, nonlinear results are promising
---
matlab/plot/plot_metric.m | 43 ++++++++++++-------
matlab/setup.m | 3 ++
matlab/write_fort90.m | 3 ++
src/geometry_mod.F90 | 26 +++++++++--
src/miller_mod.F90 | 14 +++---
testcases/miller_example/{bug => }/fort_01.90 | 2 +-
.../{bug/fort_00.90 => fort_02.90} | 6 +--
wk/analysis_gyacomo.m | 14 +++---
wk/header_3D_results.m | 2 +-
wk/lin_ITG.m | 25 ++++++-----
10 files changed, 91 insertions(+), 47 deletions(-)
rename testcases/miller_example/{bug => }/fort_01.90 (98%)
rename testcases/miller_example/{bug/fort_00.90 => fort_02.90} (96%)
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index d548c762..377c80fc 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,4 +1,4 @@
-function [ fig ] = plot_metric( data )
+function [ fig ] = plot_metric( data, options )
names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
'gxx','gxy','gyy','gyz','gzz','hatB','hatR','hatZ'};
@@ -8,23 +8,34 @@ for i_ = 1:numel(names)
namae = names{i_};
geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end
-fig = figure;
-subplot(311)
- for i = 1:5
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
- end
- xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
+NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
+if NPLOT > 0
+ fig = figure;
+ if options.SHOW_METRICS
+ subplot(311)
+ for i = 1:5
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
-subplot(312)
- for i = 6:10
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
- end
- xlim([min(data.z),max(data.z)]); legend('show');
+ subplot(312)
+ for i = 6:10
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
-subplot(313)
- for i = 11:13
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ subplot(313)
+ for i = 11:13
+ plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
end
- xlim([min(data.z),max(data.z)]); legend('show');
+ if options.SHOW_FLUXSURF
+ R = [geo_arrays(1,:,12),geo_arrays(1,1,12)];
+ Z = [geo_arrays(1,:,13),geo_arrays(1,1,13)];
+ plot(R,Z);
+ axis equal
+ end
+end
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 26618413..3b8be6c7 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -18,6 +18,9 @@ GEOM.geom = ['''',GEOMETRY,''''];
GEOM.q0 = Q0; % q factor
GEOM.shear = SHEAR; % shear
GEOM.eps = EPS; % inverse aspect ratio
+GEOM.kappa = KAPPA; % elongation
+GEOM.delta = DELTA; % triangularity
+GEOM.zeta = ZETA; % squareness
% Model parameters
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index d42c7800..29304ef0 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -30,6 +30,9 @@ fprintf(fid,[' geom = ', GEOM.geom,'\n']);
fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']);
fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']);
fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
+fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
+fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
+fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 4921117b..6f478b7c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -78,6 +78,7 @@ CONTAINS
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
INTEGER :: fid
+ real(dp) :: G1,G2,G3,Cx,Cy
! Allocate arrays
CALL geometry_allocate_mem
@@ -121,10 +122,29 @@ CONTAINS
ENDDO
ENDDO
ENDDO
+ ! Curvature operator (Frei et al. 2022 eq 2.15)
+ DO iz = izgs,izge
+ G1 = gxy(iz,eo)*gxy(iz,eo)-gxx(iz,eo)*gyy(iz,eo)
+ G2 = gxy(iz,eo)*gxz(iz,eo)-gxx(iz,eo)*gyz(iz,eo)
+ G3 = gyy(iz,eo)*gxz(iz,eo)-gxy(iz,eo)*gyz(iz,eo)
+ Cx = (G1*gradyB(iz,eo) + G2*gradzB(iz,eo))/hatB(iz,eo)
+ Cy = (G3*gradzB(iz,eo) - G1*gradxB(iz,eo))/hatB(iz,eo)
+
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff(iz,eo) = 1._dp / jacobian(iz,eo) / hatB(iz,eo)
+ ! Factor in front of the nonlinear term
+ hatB_NL(iz,eo) = Jacobian(iz,eo)&
+ *(gxx(iz,eo)*gyy(iz,eo) - gxy(iz,eo)**2)/hatB(iz,eo)
+ ENDDO
ENDDO
- ! Factor in front of the nonlinear term
- hatB_NL(izgs:izge,0:1) = Jacobian(izgs:izge,0:1)&
- *(gxx(izgs:izge,0:1)*gyy(izgs:izge,0:1) - gxy(izgs:izge,0:1)**2)/hatB(izgs:izge,0:1)
+
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index b33e3a57..eb3a9e5f 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -542,9 +542,9 @@ CONTAINS
! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = izgs,izge
- G1 = gxx_(iz,eo)*gyy_(iz,eo)-gxy_(iz,eo)*gxy_(iz,eo)
- G2 = gxx_(iz,eo)*gyz_(iz,eo)-gxy_(iz,eo)*gxz_(iz,eo)
- G3 = gxy_(iz,eo)*gyz_(iz,eo)-gyy_(iz,eo)*gxz_(iz,eo)
+ G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
+ G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
+ G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo)
Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo)
@@ -552,13 +552,13 @@ CONTAINS
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/Bfield_(iz,eo)
+ Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
ENDDO
ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff_(iz,eo) = 1._dp / jacobian_(iz,eo) / Bfield_(iz,eo)
ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff_(iz,eo) = 1._dp / Jacobian_(iz,eo) / Bfield_(iz,eo)
- enddo
+ ENDDO
contains
diff --git a/testcases/miller_example/bug/fort_01.90 b/testcases/miller_example/fort_01.90
similarity index 98%
rename from testcases/miller_example/bug/fort_01.90
rename to testcases/miller_example/fort_01.90
index 7deb16c3..9e3c9809 100644
--- a/testcases/miller_example/bug/fort_01.90
+++ b/testcases/miller_example/fort_01.90
@@ -1,7 +1,7 @@
&BASIC
nrun = 100000000
dt = 0.01
- tmax = 400
+ tmax = 100
maxruntime = 356400
/
&GRID
diff --git a/testcases/miller_example/bug/fort_00.90 b/testcases/miller_example/fort_02.90
similarity index 96%
rename from testcases/miller_example/bug/fort_00.90
rename to testcases/miller_example/fort_02.90
index 728a3456..e0af42a2 100644
--- a/testcases/miller_example/bug/fort_00.90
+++ b/testcases/miller_example/fort_02.90
@@ -1,6 +1,6 @@
&BASIC
nrun = 100000000
- dt = 0.01
+ dt = 0.0075
tmax = 100
maxruntime = 356400
/
@@ -38,7 +38,7 @@
write_Sapj = .false.
write_dens = .true.
write_temp = .true.
- job2load = -1
+ job2load = 1
/
&MODEL_PAR
! Collisionality
@@ -49,7 +49,7 @@
mu_x = 0.1
mu_y = 0.1
N_HD = 4
- mu_z = 0.2
+ mu_z = 1.0
mu_p = 0
mu_j = 0
nu = 0.2
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 65ba29bd..54353a90 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -28,7 +28,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-i_ = 1;
+i_ = 3;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
@@ -59,14 +59,14 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
-% options.NAME = '\omega_z';
-% options.NAME = 'N_e^{00}';
+% options.NAME = '\phi';
+options.NAME = '\omega_z';
+% options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -206,7 +206,9 @@ end
if 0
%% Metric infos
-fig = plot_metric(data);
+options.SHOW_FLUXSURF = 1;
+options.SHOW_METRICS = 0;
+fig = plot_metric(data,options);
end
if 0
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 189cc0e9..779dd117 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -59,6 +59,6 @@ gyacomodir = '/home/ahoffman/gyacomo/';
% resdir = 'NL_KBM/192x64x24x5x3';
%% Linear CBC
% resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
-resdir = 'testcases/miller_example_bug';
+resdir = 'testcases/miller_example';
JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index cffb9cc6..ab514b20 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -8,13 +8,15 @@ RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.0; % Collision frequency
+NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 2.22; % ele Density '''
K_Te = 6.96; % ele Temperature '''
@@ -31,23 +33,26 @@ PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 11; % real space x-gridpoints
-NY = 32; % '' y-gridpoints
+NX = 20; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.05; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
+LY = 2*pi/0.3; % Size of the squared frequency domain
+NZ = 32; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
+KAPPA = 0.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 50; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 1; % Maximal time unit
+DT = 3e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 2D arrays
@@ -125,7 +130,7 @@ end
if 1
%% Ballooning plot
options.time_2_plot = [120];
-options.kymodes = [0.9];
+options.kymodes = [0.3];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
--
GitLab