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

--------- Co-authored-by: Lucas <[email protected]> Co-authored-by: John Omotani <[email protected]>
mabarnes · Mar 28, 2024 · 44c914e · 44c914e
1 parent 79e7cf8
commit 44c914e
Show file tree

Hide file tree

Showing 64 changed files with 354 additions and 201 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
@@ -88,7 +88,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/kinetic-electrons/wall+sheath-bc_boltzmann_loworder.toml b/examples/kinetic-electrons/wall+sheath-bc_boltzmann_loworder.toml
@@ -82,11 +82,15 @@ vz_discretization = "chebyshev_pseudospectral"
 [output]
 ascii_output = true
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 #moment_dissipation_coefficient = 0.0001
 #moment_dissipation_coefficient = 1.0
-vpa_dissipation_coefficient = 0.002
-#vpa_dissipation_coefficient = 0.1
+#vpa_dissipation_coefficient = 0.002
+vpa_dissipation_coefficient = 0.1
 #vpa_dissipation_coefficient = 0.2
 #vpa_dissipation_coefficient = 2.0
 #vpa_dissipation_coefficient = 20.0
+
+[electron_numerical_dissipation]
+
+[neutral_numerical_dissipation]
diff --git a/examples/kinetic-electrons/wall+sheath-bc_kinetic.toml b/examples/kinetic-electrons/wall+sheath-bc_kinetic.toml
@@ -79,7 +79,7 @@ vz_discretization = "chebyshev_pseudospectral"
 [output]
 ascii_output = true
 
-[numerical_dissipation]
+[electron_numerical_dissipation]
 #moment_dissipation_coefficient = 0.0001
 #moment_dissipation_coefficient = 1.0
 #vpa_dissipation_coefficient = 0.002

diff --git a/examples/kinetic-electrons/wall+sheath-bc_kinetic_krook_loworder.toml b/examples/kinetic-electrons/wall+sheath-bc_kinetic_krook_loworder.toml
@@ -80,7 +80,7 @@ vz_discretization = "chebyshev_pseudospectral"
 [output]
 ascii_output = true
 
-[numerical_dissipation]
+[electron_numerical_dissipation]
 #moment_dissipation_coefficient = 0.0001
 #moment_dissipation_coefficient = 1.0
 #vpa_dissipation_coefficient = 0.002

diff --git a/examples/kinetic-electrons/wall+sheath-bc_kinetic_loworder.toml b/examples/kinetic-electrons/wall+sheath-bc_kinetic_loworder.toml
@@ -82,7 +82,7 @@ vz_discretization = "chebyshev_pseudospectral"
 [output]
 ascii_output = true
 
-[numerical_dissipation]
+[electron_numerical_dissipation]
 #moment_dissipation_coefficient = 0.0001
 #moment_dissipation_coefficient = 1.0
 #vpa_dissipation_coefficient = 0.002

diff --git a/examples/numerical-dissipation/num-diss-relaxation.toml b/examples/numerical-dissipation/num-diss-relaxation.toml
@@ -58,7 +58,7 @@ vperp_L = 3.0
 vperp_bc = "zero"
 vperp_discretization = "gausslegendre_pseudospectral"
 
-[numerical_dissipation]
+[ion_numerical_dissipation]
 vpa_dissipation_coefficient = 0.1
 vperp_dissipation_coefficient = 0.1
 z_dissipation_coefficient = -1.0

diff --git a/examples/wall-bc/wall-bc_cheb.toml b/examples/wall-bc/wall-bc_cheb.toml
@@ -65,6 +65,10 @@ 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
 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 @@ 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
 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 @@ 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
 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 @@ 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
 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
@@ -78,7 +78,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
@@ -78,7 +78,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
@@ -78,7 +78,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
@@ -78,7 +78,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
@@ -85,8 +85,14 @@ 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
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split1.toml b/examples/wall-bc/wall-bc_volumerecycle_split1.toml
@@ -85,8 +85,14 @@ 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
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split2.toml b/examples/wall-bc/wall-bc_volumerecycle_split2.toml
@@ -85,8 +85,14 @@ 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
diff --git a/examples/wall-bc/wall-bc_volumerecycle_split3.toml b/examples/wall-bc/wall-bc_volumerecycle_split3.toml
@@ -85,8 +85,14 @@ 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
diff --git a/moment_kinetics/debug_test/restart_interpolation_inputs.jl b/moment_kinetics/debug_test/restart_interpolation_inputs.jl
@@ -60,7 +60,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/electron_fluid_equations.jl b/moment_kinetics/src/electron_fluid_equations.jl
@@ -138,7 +138,7 @@ function electron_energy_equation!(ppar, dens_i, fvec, moments, collisions, dt,
     # end
     # compute the contribution to the rhs of the energy equation
     # arising from artificial diffusion
-    diffusion_coefficient = num_diss_params.moment_dissipation_coefficient
+    diffusion_coefficient = num_diss_params.electron.moment_dissipation_coefficient
     if diffusion_coefficient > 0.0
         @loop_r_z ir iz begin
             ppar[iz,ir] += dt*diffusion_coefficient*moments.electron.d2ppar_dz2[iz,ir]

diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl
@@ -268,7 +268,7 @@ function update_electron_pdf_with_time_advance!(fvec, pdf, qpar, qpar_updated,
                              buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4,
                              buffer_r_5, buffer_r_6, z_spectral, z)
             # centred second derivative for dissipation
-            if num_diss_params.moment_dissipation_coefficient > 0.0
+            if num_diss_params.electron.moment_dissipation_coefficient > 0.0
                 @views derivative_z!(moments.electron.d2ppar_dz2, dppar_dz, buffer_r_1,
                                  buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z)
                 begin_serial_region()
@@ -327,7 +327,7 @@ function update_electron_pdf_with_time_advance!(fvec, pdf, qpar, qpar_updated,
         enforce_boundary_condition_on_electron_pdf!(pdf, phi, vthe, upar, vperp, vpa,
                                                     vperp_spectral, vpa_spectral,
                                                     vpa_advect, moments,
-                                                    num_diss_params.vpa_dissipation_coefficient > 0.0,
+                                                    num_diss_params.electron.vpa_dissipation_coefficient > 0.0,
                                                     composition.me_over_mi)
         #println("A pdf 1 ", pdf[:,1,1,1])
         #println("A pdf end ", pdf[:,1,end,1])
@@ -379,7 +379,7 @@ function update_electron_pdf_with_time_advance!(fvec, pdf, qpar, qpar_updated,
                              buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4,
                              buffer_r_5, buffer_r_6, z_spectral, z)
          # centred second derivative for dissipation
-        if num_diss_params.moment_dissipation_coefficient > 0.0
+        if num_diss_params.electron.moment_dissipation_coefficient > 0.0
             @views derivative_z!(moments.electron.d2ppar_dz2, dppar_dz, buffer_r_1,
                                  buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z)
         end
@@ -1442,15 +1442,15 @@ function add_dissipation_term!(residual, pdf, scratch_dummy, z_spectral, z, vpa,
     #                         buffer_r_4, z_spectral, z)
     #    @views derivative_z!(dummy_zr2, dummy_zr1, buffer_r_1, buffer_r_2, buffer_r_3,
     #                         buffer_r_4, z_spectral, z)
-    #    @. residual[ivpa,ivperp,:,:] -= num_diss_params.z_dissipation_coefficient * dummy_zr2
+    #    @. residual[ivpa,ivperp,:,:] -= num_diss_params.electron.z_dissipation_coefficient * dummy_zr2
     #end
     begin_r_z_vperp_region()
     @loop_r_z_vperp ir iz ivperp begin
         #@views derivative!(vpa.scratch, pdf[:,ivperp,iz,ir], vpa, false)
         #@views derivative!(vpa.scratch2, vpa.scratch, vpa, false)
-        #@. residual[:,ivperp,iz,ir] -= num_diss_params.vpa_dissipation_coefficient * vpa.scratch2
+        #@. residual[:,ivperp,iz,ir] -= num_diss_params.electron.vpa_dissipation_coefficient * vpa.scratch2
         @views second_derivative!(vpa.scratch, pdf[:,ivperp,iz,ir], vpa, vpa_spectral)
-        @. residual[:,ivperp,iz,ir] -= num_diss_params.vpa_dissipation_coefficient * vpa.scratch
+        @. residual[:,ivperp,iz,ir] -= num_diss_params.electron.vpa_dissipation_coefficient * vpa.scratch
     end
     #stop()
     return nothing

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]

diff --git a/moment_kinetics/src/force_balance.jl b/moment_kinetics/src/force_balance.jl
@@ -38,7 +38,7 @@ function force_balance!(pflx, density_out, fvec, moments, fields, collisions, dt
     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
             pflx[iz,ir,is] += dt*diffusion_coefficient*moments.ion.d2upar_dz2[iz,ir,is]*density[iz,ir,is]
@@ -92,7 +92,7 @@ function neutral_force_balance!(pflx, density_out, fvec, moments, fields, collis
     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
             pflx[iz,ir,isn] += dt*diffusion_coefficient*moments.neutral.d2uz_dz2[iz,ir,isn]*density[iz,ir,isn]