Separate dissipation parameters for ions, electrons, neutrals (#192)

mabarnes · Apr 26, 2024 · 95f7ece · 95f7ece
1 parent 4030583
commit 95f7ece
Show file tree

Hide file tree

Showing 59 changed files with 319 additions and 174 deletions.
diff --git a/.gitignore b/.gitignore
@@ -26,3 +26,4 @@ makie_postproc.so
 plots_postproc.so
 precompile-temp
 post_processing_input.toml
+.DS_Store
diff --git a/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml b/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml
@@ -86,7 +86,7 @@ vperp_discretization = "gausslegendre_pseudospectral"
 #vzeta_bc = "periodic"
 #vzeta_discretization = "chebyshev_pseudospectral"
 
-#[numerical_dissipation]
+#[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 0.0
 #vperp_dissipation_coefficient = 0.0
 #z_dissipation_coefficient = 0.1

diff --git a/examples/geometry/1D-mirror.toml b/examples/geometry/1D-mirror.toml
@@ -71,7 +71,7 @@ vz_L = 18.0
 vz_bc = "both_zero"
 vz_discretization = "chebyshev_pseudospectral"
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 force_minimum_pdf_value = 0.0

diff --git a/examples/gk-ions/2D-periodic-gk.toml b/examples/gk-ions/2D-periodic-gk.toml
@@ -73,15 +73,19 @@ vz_L = 18.0
 vz_bc = "both_zero"
 vz_discretization = "chebyshev_pseudospectral"
 
-[numerical_dissipation]
-vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
-vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
-#r_disspipation_coefficient = 1.0e-3
-#force_minimum_pdf_value = 0.0
-
 [geometry]
 #option="1D-mirror"
 DeltaB=0.0
 option="constant-helical"
 pitch=0.1
 rhostar= 0.1
+
+[ion_numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+#r_disspipation_coefficient = 1.0e-3
+#force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+moment_dissipation_coefficient = 0.0001
+vz_dissipation_coefficient = 0.01
diff --git a/examples/numerical-dissipation/num-diss-relaxation.toml b/examples/numerical-dissipation/num-diss-relaxation.toml
@@ -58,7 +58,7 @@ nwrite = 2000
 nwrite_dfns = 2000
 split_operators = false
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 vpa_dissipation_coefficient = 0.1
 vperp_dissipation_coefficient = 0.1
 z_dissipation_coefficient = -1.0

diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml
@@ -81,8 +81,14 @@ z_width = 0.125
 source_strength = 2.0
 source_T = 2.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e-1
+#vpa_dissipation_coefficient = 1.0e-2
+#vpa_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml
@@ -81,8 +81,14 @@ z_width = 0.125
 source_strength = 2.0
 source_T = 2.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vpa_dissipation_coefficient = 1.0e-1
+#vpa_dissipation_coefficient = 1.0e-2
+#vpa_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml
@@ -87,7 +87,13 @@ z_width = 0.125
 source_strength = 2.0
 source_T = 2.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
+#vpa_dissipation_coefficient = 1.0e-1
+#vpa_dissipation_coefficient = 1.0e-2
+#vpa_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete43.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete43.toml
@@ -87,7 +87,13 @@ z_width = 0.125
 source_strength = 2.0
 source_T = 2.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
+#vpa_dissipation_coefficient = 1.0e-1
+#vpa_dissipation_coefficient = 1.0e-2
+#vpa_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/wall-bc/wall-bc_cheb.toml b/examples/wall-bc/wall-bc_cheb.toml
@@ -65,6 +65,10 @@ dt = 1.0e-5
 nwrite = 1000
 split_operators = false
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+force_minimum_pdf_value = 0.0
diff --git a/examples/wall-bc/wall-bc_cheb_split1.toml b/examples/wall-bc/wall-bc_cheb_split1.toml
@@ -66,6 +66,10 @@ dt = 1.0e-5
 nwrite = 1000
 split_operators = false
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+force_minimum_pdf_value = 0.0
diff --git a/examples/wall-bc/wall-bc_cheb_split2.toml b/examples/wall-bc/wall-bc_cheb_split2.toml
@@ -66,6 +66,10 @@ dt = 1.0e-5
 nwrite = 1000
 split_operators = false
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+force_minimum_pdf_value = 0.0
diff --git a/examples/wall-bc/wall-bc_cheb_split3.toml b/examples/wall-bc/wall-bc_cheb_split3.toml
@@ -66,6 +66,10 @@ dt = 1.0e-5
 nwrite = 1000
 split_operators = false
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
 force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+force_minimum_pdf_value = 0.0
diff --git a/examples/wall-bc/wall-bc_no-neutrals.toml b/examples/wall-bc/wall-bc_no-neutrals.toml
@@ -77,7 +77,7 @@ z_width = 0.125
 source_strength = 8.0
 source_T = 1.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/wall-bc/wall-bc_no-neutrals_split1.toml b/examples/wall-bc/wall-bc_no-neutrals_split1.toml
@@ -77,7 +77,7 @@ z_width = 0.125
 source_strength = 8.0
 source_T = 1.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/wall-bc/wall-bc_no-neutrals_split2.toml b/examples/wall-bc/wall-bc_no-neutrals_split2.toml
@@ -77,7 +77,7 @@ z_width = 0.125
 source_strength = 8.0
 source_T = 1.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/wall-bc/wall-bc_no-neutrals_split3.toml b/examples/wall-bc/wall-bc_no-neutrals_split3.toml
@@ -77,7 +77,7 @@ z_width = 0.125
 source_strength = 8.0
 source_T = 1.0
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3

diff --git a/examples/wall-bc/wall-bc_volumerecycle.toml b/examples/wall-bc/wall-bc_volumerecycle.toml
@@ -84,12 +84,18 @@ active = true
 source_type = "recycling"
 recycling_controller_fraction = 0.25
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 #force_minimum_pdf_value = 0.0
 
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-1
+#vz_dissipation_coefficient = 1.0e-2
+#vz_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
+
 [krook_collisions]
 use_krook = true
 frequency_option = "reference_parameters"
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split1.toml b/examples/wall-bc/wall-bc_volumerecycle_split1.toml
@@ -84,12 +84,18 @@ active = true
 source_type = "recycling"
 recycling_controller_fraction = 0.25
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 #force_minimum_pdf_value = 0.0
 
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-1
+#vz_dissipation_coefficient = 1.0e-2
+#vz_dissipation_coefficient = 1.0e-3
+#force_minimum_pdf_value = 0.0
+
 [krook_collisions]
 use_krook = true
 frequency_option = "reference_parameters"
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split2.toml b/examples/wall-bc/wall-bc_volumerecycle_split2.toml
@@ -84,12 +84,18 @@ active = true
 source_type = "recycling"
 recycling_controller_fraction = 0.25
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 #force_minimum_pdf_value = 0.0
 
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-1
+#vz_dissipation_coefficient = 1.0e-2
+#vz_dissipation_coefficient = 1.0e-3
+#force_minimum_pdf_value = 0.0
+
 [krook_collisions]
 use_krook = true
 frequency_option = "reference_parameters"
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split3.toml b/examples/wall-bc/wall-bc_volumerecycle_split3.toml
@@ -84,12 +84,18 @@ active = true
 source_type = "recycling"
 recycling_controller_fraction = 0.25
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #vpa_dissipation_coefficient = 1.0e-1
 #vpa_dissipation_coefficient = 1.0e-2
 #vpa_dissipation_coefficient = 1.0e-3
 force_minimum_pdf_value = 0.0
 
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-1
+#vz_dissipation_coefficient = 1.0e-2
+#vz_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
+
 [krook_collisions]
 use_krook = true
 frequency_option = "reference_parameters"
diff --git a/moment_kinetics/debug_test/restart_interpolation_inputs.jl b/moment_kinetics/debug_test/restart_interpolation_inputs.jl
@@ -59,7 +59,8 @@ base_input = Dict(
      "vzeta_nelement" => 1,
      "vr_ngrid" => 1,
      "vr_nelement" => 1,
-     "numerical_dissipation" => Dict{String,Any}("force_minimum_pdf_value" => 0.0))
+     "ion_numerical_dissipation" => Dict{String,Any}("force_minimum_pdf_value" => 0.0),
+     "neutral_numerical_dissipation" => Dict{String,Any}("force_minimum_pdf_value" => 0.0))
 
 test_input =
     merge(base_input,

diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl
@@ -638,31 +638,31 @@ using IfElse
             Si += - nuii_krook*(FMaxwellian - dfni)
         end
         include_num_diss_in_MMS = true
-        if num_diss_params.vpa_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Si += - num_diss_params.vpa_dissipation_coefficient*Dvpa(Dvpa(dfni))
+        if num_diss_params.ion.vpa_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Si += - num_diss_params.ion.vpa_dissipation_coefficient*Dvpa(Dvpa(dfni))
         end
-        if num_diss_params.vperp_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Si += - num_diss_params.vperp_dissipation_coefficient*Dvperp(Dvperp(dfni))
+        if num_diss_params.ion.vperp_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Si += - num_diss_params.ion.vperp_dissipation_coefficient*Dvperp(Dvperp(dfni))
         end
-        if num_diss_params.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Si += - rfac*num_diss_params.r_dissipation_coefficient*Dr(Dr(dfni))
+        if num_diss_params.ion.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Si += - rfac*num_diss_params.ion.r_dissipation_coefficient*Dr(Dr(dfni))
         end
-        if num_diss_params.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Si += - num_diss_params.z_dissipation_coefficient*Dz(Dz(dfni))
+        if num_diss_params.ion.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Si += - num_diss_params.ion.z_dissipation_coefficient*Dz(Dz(dfni))
         end
 
         Source_i = expand_derivatives(Si)
 
         # the neutral source to maintain the manufactured solution
         Sn = Dt(dfnn) + vz * Dz(dfnn) + rfac*vr * Dr(dfnn) + cx_frequency* (densi*dfnn - densn*vrvzvzeta_dfni) + ionization_frequency*dense*dfnn
-        if num_diss_params.vz_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Sn += - num_diss_params.vz_dissipation_coefficient*Dvz(Dvz(dfnn))
+        if num_diss_params.neutral.vz_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Sn += - num_diss_params.neutral.vz_dissipation_coefficient*Dvz(Dvz(dfnn))
         end
-        if num_diss_params.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Sn += - rfac*num_diss_params.r_dissipation_coefficient*Dr(Dr(dfnn))
+        if num_diss_params.neutral.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Sn += - rfac*num_diss_params.neutral.r_dissipation_coefficient*Dr(Dr(dfnn))
         end
-        if num_diss_params.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
-            Sn += - num_diss_params.z_dissipation_coefficient*Dz(Dz(dfnn))
+        if num_diss_params.neutral.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS
+            Sn += - num_diss_params.neutral.z_dissipation_coefficient*Dz(Dz(dfnn))
         end
 
         Source_n = expand_derivatives(Sn)

diff --git a/moment_kinetics/src/continuity.jl b/moment_kinetics/src/continuity.jl
@@ -39,7 +39,7 @@ function continuity_equation!(dens_out, fvec_in, moments, composition, dt, spect
     end
 
     # Ad-hoc diffusion to stabilise numerics...
-    diffusion_coefficient = num_diss_params.moment_dissipation_coefficient
+    diffusion_coefficient = num_diss_params.ion.moment_dissipation_coefficient
     if diffusion_coefficient > 0.0
         @loop_s_r_z is ir iz begin
             dens_out[iz,ir,is] += dt*diffusion_coefficient*moments.ion.d2dens_dz2[iz,ir,is]
@@ -80,7 +80,7 @@ function neutral_continuity_equation!(dens_out, fvec_in, moments, composition, d
     end
 
     # Ad-hoc diffusion to stabilise numerics...
-    diffusion_coefficient = num_diss_params.moment_dissipation_coefficient
+    diffusion_coefficient = num_diss_params.neutral.moment_dissipation_coefficient
     if diffusion_coefficient > 0.0
         @loop_sn_r_z isn ir iz begin
             dens_out[iz,ir,isn] += dt*diffusion_coefficient*moments.neutral.d2dens_dz2[iz,ir,isn]

diff --git a/moment_kinetics/src/energy_equation.jl b/moment_kinetics/src/energy_equation.jl
@@ -30,7 +30,7 @@ function energy_equation!(ppar, fvec, moments, collisions, dt, spectral, composi
         end
     end
 
-    diffusion_coefficient = num_diss_params.moment_dissipation_coefficient
+    diffusion_coefficient = num_diss_params.ion.moment_dissipation_coefficient
     if diffusion_coefficient > 0.0
         @loop_s_r_z is ir iz begin
             ppar[iz,ir,is] += dt*diffusion_coefficient*moments.ion.d2ppar_dz2[iz,ir,is]
@@ -82,7 +82,7 @@ function neutral_energy_equation!(pz, fvec, moments, collisions, dt, spectral,
         end
     end
 
-    diffusion_coefficient = num_diss_params.moment_dissipation_coefficient
+    diffusion_coefficient = num_diss_params.neutral.moment_dissipation_coefficient
     if diffusion_coefficient > 0.0
         @loop_sn_r_z isn ir iz begin
             pz[iz,ir,isn] += dt*diffusion_coefficient*moments.neutral.d2pz_dz2[iz,ir,isn]