Skip to content

Commit

Permalink
Update Efield functionalities (#109)
Browse files Browse the repository at this point in the history
* Add control over GMRES solve absolute tolerance.

* NL solve state data need as many ghost as regular state data.

* Add a kernel for MOL-type second order electrons advection.

* Add accessor to vector of BGcharge and scale the right number of ghost
cells when dealing with the nlState

* Minor update for AMReX syntax.

* Add ParmParse of nE advection scheme order.

* Make smallvel in ABecCec laplacian consistent with advection scheme.

* Add control over the PC diffdrift operator damping and add option
for second order finite difference JtV.

* - 2nd order FD in JtV
 - 2nd order advection scheme for electrons
 - implement control of ABecCec damping
 - FillPatch nE in NLresidual, done in LinOp for phiV
 - propagate use of nGrowState throughout the NL solve.

* Update missing function and update exitNewton to exit if newton increment become too small.

* Parse new options.

* FillBoundary of Udrift.

* Fixes to ABecCec for GPUs and 3D.

* Enable user-defined Ke and prepare for tabulation.

* First pass at Ke Tabulation.

* Move nE and phiV into level data state MF

* Restart from Plt for efield.

* GPU fix.

* Couple of fixes for restarting from Plt or from non-EF cases.

* Fully enforce resetting the time if activated.

* Use AMREX OMP macro.

* Add the Kuhl bunsen case.

* Restore phiV boundaries.

* Fix restart from non-EF

* Fix using dcomp in phiV BCfill functor.

* Enlarge the electrode a bit and fix i,j,k mixup in zero_visc

* Access to more MLMG option for preconditioner.

* Split diffusion solve between neutral species and individual solve
for each ions (because BC type might change).

* Minor updates of Kuhl case.

* Fix merge blunder and update plt restart for efield.

* Remove unsused BC function.

* Fix EF derived.

* Update the IonizedWave case.

* Update Ionizedwave input files.

* Separate control of the MG depth for the preconditioner matrices.

* Add GMRES final iter count to verbose when = 2

* Update FlameSheetIons

* Expose a few more control options for efield in FlameSheetIons case input file.

* Don't restart.
  • Loading branch information
esclapez committed Jul 20, 2022
1 parent a92b9fa commit ee45a90
Show file tree
Hide file tree
Showing 46 changed files with 1,611 additions and 505 deletions.
52 changes: 28 additions & 24 deletions Exec/Efield/FlameSheetIons/input.2d-regt
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
geometry.is_periodic = 1 0 # For each dir, 0: non-perio, 1: periodic
geometry.coord_sys = 0 # 0 => cart, 1 => RZ
geometry.prob_lo = 0.0 0.0 0.0 # x_lo y_lo (z_lo)
geometry.prob_hi = 0.004 0.016 0.016 # x_hi y_hi (z_hi)
geometry.prob_hi = 0.008 0.016 0.016 # x_hi y_hi (z_hi)

# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<<
# Interior, Inflow, Outflow, Symmetry,
Expand All @@ -12,12 +12,12 @@ peleLM.hi_bc = Interior Outflow


#-------------------------AMR CONTROL----------------------------
amr.n_cell = 32 128 32 # Level 0 number of cells in each direction
amr.n_cell = 64 128 32 # Level 0 number of cells in each direction
amr.v = 1 # AMR verbose
amr.max_level = 2 # maximum level number allowed
amr.ref_ratio = 2 2 2 2 # refinement ratio
amr.regrid_int = 2 # how often to regrid
amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est
amr.n_error_buf = 2 4 2 2 # number of buffer cells in error est
amr.grid_eff = 0.7 # what constitutes an efficient grid
amr.blocking_factor = 16 # block factor in grid generation (min box size)
amr.max_grid_size = 128 # max box size
Expand All @@ -26,52 +26,56 @@ amr.max_grid_size = 128 # max box size
#--------------------------- Problem -------------------------------
prob.P_mean = 101325.0
prob.standoff = -.01
prob.pertmag = 0.0000
prob.pertmag = 0.0001
pmf.datafile = "pmf_CH4Air_1p0.dat"
pmf.do_cellAverage = 0
prob.PhiV_y_hi = 1000.0

#-------------------------PeleLM CONTROL----------------------------
peleLM.v = 1
peleLM.v = 2
peleLM.incompressible = 0
peleLM.rho = 1.17
peleLM.mu = 0.0
peleLM.use_wbar = 1
peleLM.sdc_iterMax = 2
peleLM.floor_species = 0
eleLM.floor_species = 0
peleLM.num_init_iter = 2
peleLM.advection_scheme = "Godunov_BDS"

#amr.restart = chk00005
#amr.check_int = 2000
#amr.restart = chk00070
amr.check_int = 20
amr.plot_int = 10
amr.max_step = 100
amr.dt_shrink = 0.01
amr.max_step = 20
amr.dt_shrink = 0.1
amr.stop_time = 0.001
#amr.stop_time = 1.00
#amr.fixed_dt = 1.957622853e-9
amr.cfl = 0.5
amr.derive_plot_vars = avg_pressure mag_vort mass_fractions chargedistrib efieldy
amr.cfl = 0.95
amr.derive_plot_vars = avg_pressure mag_vort mass_fractions chargedistrib efieldx efieldy LorentzFx LorentzFy

peleLM.chem_integrator = "ReactorCvode"
peleLM.use_typ_vals_chem = 1 # Use species/temp typical values in CVODE
ode.rtol = 1.0e-6 # Relative tolerance of the chemical solve
ode.atol = 1.0e-5 # Absolute tolerance factor applied on typical values
ode.rtol = 1.0e-7 # Relative tolerance of the chemical solve
ode.atol = 1.0e-6 # Absolute tolerance factor applied on typical values
cvode.solve_type = denseAJ_direct # CVODE Linear solve type (for Newton direction)
cvode.max_order = 4 # CVODE max BDF order.

godunov.use_ppm = 0

#--------------------REFINEMENT CONTROL------------------------
ef.phiV_lo_bc = Interior Dirichlet
ef.phiV_hi_bc = Interior Dirichlet
ef.phiV_polarity_lo = Neutral Cathode
ef.phiV_polarity_hi = Neutral Anode
ef.GMRES_rel_tol = 1.0e-8
ef.JFNK_lambda = 1.0e-6
ef.precond.diff_verbose = 0
ef.GMRES_rel_tol = 1.0e-6
ef.GMRES_abs_tol = 1.0e-13
ef.JFNK_lambda = 5.0e-7
ef.JFNK_diffType = 1
ef.PC_approx = 2
#ef.PC_damping = 0.95
ef.advection_scheme_order = 1
ef.precond.diff_verbose = 2
ef.precond.Stilda_verbose = 0
ef.precond.fixedIter = 10
ef.precond.fixedIter = 3
#ef.precond.max_coarsening_level_diff = 3
gmres.krylovBasis_size = 20
gmres.verbose = 2
gmres.verbose = 0
gmres.max_restart = 3

#--------------------REFINEMENT CONTROL------------------------
Expand All @@ -87,7 +91,7 @@ gmres.max_restart = 3

amr.refinement_indicators = yE
amr.yE.max_level = 3
amr.yE.value_greater = 1.0e16
amr.yE.value_greater = 1.0e17
amr.yE.field_name = nE

#amrex.fpe_trap_invalid = 1
Expand Down
8 changes: 3 additions & 5 deletions Exec/Efield/FlameSheetIons/pelelm_prob.H
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,6 @@ void pelelm_initdata(int i, int j, int k,
int is_incompressible,
amrex::Array4<amrex::Real> const& state,
amrex::Array4<amrex::Real> const& aux,
amrex::Array4<amrex::Real> const& nE,
amrex::Array4<amrex::Real> const& phiV,
amrex::GeometryData const& geomdata,
ProbParm const& prob_parm,
pele::physics::PMF::PmfData::DataContainer const * pmf_data)
Expand Down Expand Up @@ -108,7 +106,7 @@ void pelelm_initdata(int i, int j, int k,
for (int n = 0; n < NUM_SPECIES; n++) {
chargeDist += massfrac[n] * zk[n];
}
nE(i,j,k) = std::max(1.0e-24,chargeDist/elemCharge);
state(i,j,k,NE) = std::max(1.0e-24,chargeDist/elemCharge);


AMREX_D_TERM(state(i,j,k,VELX) = 0.0;,
Expand All @@ -132,8 +130,8 @@ void pelelm_initdata(int i, int j, int k,
state(i,j,k,FIRSTSPEC+n) = massfrac[n] * state(i,j,k,DENSITY);
}

nE(i,j,k) *= state(i,j,k,DENSITY);
phiV(i,j,k) = 0.0;
state(i,j,k,NE) *= state(i,j,k,DENSITY);
state(i,j,k,PHIV) = 0.0;
}

AMREX_GPU_DEVICE
Expand Down
2 changes: 2 additions & 0 deletions Exec/Efield/FlameSheetIons/pelelm_prob.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,10 @@ void PeleLM::readProbParm()
pp.query("standoff", PeleLM::prob_parm->standoff);
pp.query("pertmag", PeleLM::prob_parm->pertmag);

#ifdef PELE_USE_EFIELD
pp.query("PhiV_y_hi", PeleLM::prob_parm->phiV_hiy);
pp.query("PhiV_y_lo", PeleLM::prob_parm->phiV_loy);
#endif

PeleLM::pmf_data.initialize();
}
44 changes: 31 additions & 13 deletions Exec/Efield/IonizedAirWave/input.2d-regt
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ peleLM.hi_bc = Interior Outflow
#-------------------------AMR CONTROL----------------------------
amr.n_cell = 32 128 32 # Level 0 number of cells in each direction
amr.v = 1 # AMR verbose
amr.max_level = 1 # maximum level number allowed
amr.max_level = 2 # maximum level number allowed
amr.ref_ratio = 2 2 2 2 # refinement ratio
amr.regrid_int = 5 # how often to regrid
amr.n_error_buf = 1 1 2 2 # number of buffer cells in error est
Expand All @@ -27,8 +27,8 @@ amr.max_grid_size = 128 # max box size
prob.P_mean = 101325.0
prob.standoff = -.023
prob.pertmag = 0.0000
prob.pmf_datafile = "IonAir_pmf.dat"
prob.PhiV_y_hi = 10.0
pmf.datafile = "IonAir_pmf.dat"
prob.PhiV_y_hi = 100.0

#-------------------------PeleLM CONTROL----------------------------
peleLM.v = 1
Expand All @@ -38,27 +38,45 @@ peleLM.mu = 0.0
peleLM.use_wbar = 1
peleLM.sdc_iterMax = 2
peleLM.floor_species = 0
peleLM.num_init_iter = 3

#amr.restart = chk00005
#amr.check_int = 2000
amr.plot_int = 100
amr.max_step = 1
amr.dt_shrink = 0.0001
amr.plot_int = 2
amr.max_step = 20
amr.dt_shrink = 0.01
amr.stop_time = 0.001
#amr.stop_time = 1.00
amr.cfl = 0.5
amr.derive_plot_vars = avg_pressure mag_vort mass_fractions
amr.derive_plot_vars = avg_pressure mag_vort mass_fractions chargedistrib efieldy

cvode.solve_type = dense # CVODE Linear solve type (for Newton direction)
ode.analytical_jacobian = 0 # Provide analytical jacobian (from Fuego)

#--------------------REFINEMENT CONTROL------------------------
#--------------------EFIELD CONTROL------------------------
ef.phiV_lo_bc = Interior Dirichlet
ef.phiV_hi_bc = Interior Dirichlet
ef.phiV_polarity_lo = Neutral Cathode
ef.phiV_polarity_hi = Neutral Anode
ef.GMRES_verbose = 2
ef.GMRES_rel_tol = 1.0e-10
ef.GMRES_rel_tol = 1.0e-6
ef.GMRES_abs_tol = 1.0e-13
ef.JFNK_lambda = 1.0e-7
#ef.JFNK_diffType = 1
#ef.JFNK_maxNewton = 10
ef.PC_approx = 2
#ef.PC_damping = 0.8
ef.advection_scheme_order = 1
ef.precond.diff_verbose = 0
ef.precond.Stilda_verbose = 0
ef.precond.fixedIter = 10
#ef.precond.max_coarsening_level = 2
#ef.precond.num_pre_smooth = 4
#ef.precond.num_post_smooth = 4
ef.tabulated_Ke = 0
ef.fixed_Ke = 0.4
gmres.krylovBasis_size = 20
gmres.verbose = 0
gmres.max_restart = 3

#--------------------REFINEMENT CONTROL------------------------
#amr.refinement_indicators = temp
Expand All @@ -76,6 +94,6 @@ amr.yE.max_level = 3
amr.yE.value_greater = 1.0e12
amr.yE.field_name = nE

amrex.fpe_trap_invalid = 1
amrex.fpe_trap_zero = 1
amrex.fpe_trap_overflow = 1
#amrex.fpe_trap_invalid = 1
#amrex.fpe_trap_zero = 1
#amrex.fpe_trap_overflow = 1
99 changes: 99 additions & 0 deletions Exec/Efield/IonizedAirWave/input.2d-regt_wave
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
#----------------------DOMAIN DEFINITION------------------------
geometry.is_periodic = 1 0 # For each dir, 0: non-perio, 1: periodic
geometry.coord_sys = 0 # 0 => cart, 1 => RZ
geometry.prob_lo = 0.0 0.0 0.0 # x_lo y_lo (z_lo)
geometry.prob_hi = 0.032 0.032 0.016 # x_hi y_hi (z_hi)

# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<<
# Interior, Inflow, Outflow, Symmetry,
# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
peleLM.lo_bc = Interior Inflow
peleLM.hi_bc = Interior Outflow


#-------------------------AMR CONTROL----------------------------
amr.n_cell = 128 128 32 # Level 0 number of cells in each direction
amr.v = 1 # AMR verbose
amr.max_level = 2 # maximum level number allowed
amr.ref_ratio = 2 2 2 2 # refinement ratio
amr.regrid_int = 5 # how often to regrid
amr.n_error_buf = 1 1 2 2 # number of buffer cells in error est
amr.grid_eff = 0.7 # what constitutes an efficient grid
amr.blocking_factor = 16 # block factor in grid generation (min box size)
amr.max_grid_size = 128 # max box size


#--------------------------- Problem -------------------------------
prob.P_mean = 101325.0
prob.standoff = -.023
prob.pertmag = 0.00025
pmf.datafile = "IonAir_pmf.dat"
prob.PhiV_y_hi = 100.0

#-------------------------PeleLM CONTROL----------------------------
peleLM.v = 1
peleLM.incompressible = 0
peleLM.rho = 1.17
peleLM.mu = 0.0
peleLM.use_wbar = 1
peleLM.sdc_iterMax = 2
peleLM.floor_species = 0
peleLM.num_init_iter = 3

#amr.restart = chk00005
#amr.check_int = 2000
amr.plot_int = 2
amr.max_step = 20
amr.dt_shrink = 0.01
amr.stop_time = 0.001
#amr.stop_time = 1.00
amr.cfl = 0.5
amr.derive_plot_vars = avg_pressure mag_vort mass_fractions chargedistrib efieldx efieldy

cvode.solve_type = dense # CVODE Linear solve type (for Newton direction)
ode.analytical_jacobian = 0 # Provide analytical jacobian (from Fuego)

#--------------------EFIELD CONTROL------------------------
ef.phiV_lo_bc = Interior Dirichlet
ef.phiV_hi_bc = Interior Dirichlet
ef.phiV_polarity_lo = Neutral Cathode
ef.phiV_polarity_hi = Neutral Anode
ef.GMRES_rel_tol = 1.0e-6
ef.GMRES_abs_tol = 1.0e-13
ef.JFNK_lambda = 1.0e-7
#ef.JFNK_diffType = 1
#ef.JFNK_maxNewton = 10
ef.PC_approx = 2
#ef.PC_damping = 0.8
ef.advection_scheme_order = 1
ef.precond.diff_verbose = 0
ef.precond.Stilda_verbose = 0
ef.precond.fixedIter = 10
#ef.precond.max_coarsening_level = 2
#ef.precond.num_pre_smooth = 4
#ef.precond.num_post_smooth = 4
ef.tabulated_Ke = 0
ef.fixed_Ke = 0.4
gmres.krylovBasis_size = 20
gmres.verbose = 0
gmres.max_restart = 3

#--------------------REFINEMENT CONTROL------------------------
#amr.refinement_indicators = temp
#amr.temp.max_level = 1
#amr.temp.value_greater = 305
#amr.temp.field_name = temp

#amr.refinement_indicators = magVort
#amr.magVort.max_level = 1
#amr.magVort.value_greater = 500.0
#amr.magVort.field_name = mag_vort

amr.refinement_indicators = yE
amr.yE.max_level = 3
amr.yE.value_greater = 1.0e12
amr.yE.field_name = nE

#amrex.fpe_trap_invalid = 1
#amrex.fpe_trap_zero = 1
#amrex.fpe_trap_overflow = 1
24 changes: 11 additions & 13 deletions Exec/Efield/IonizedAirWave/pelelm_prob.H
Original file line number Diff line number Diff line change
Expand Up @@ -7,21 +7,19 @@

#include <PeleLM_Index.H>
#include <pelelm_prob_parm.H>
#include <pmf.H>
#include <pmf_data.H>
#include <PMF.H>
#include <PMFData.H>
#include <PelePhysics.H>

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void pelelm_initdata(int i, int j, int k,
int is_incompressible,
int /*is_incompressible*/,
amrex::Array4<amrex::Real> const& state,
amrex::Array4<amrex::Real> const& aux,
amrex::Array4<amrex::Real> const& nE,
amrex::Array4<amrex::Real> const& phiV,
amrex::Array4<amrex::Real> const& /*aux*/,
amrex::GeometryData const& geomdata,
ProbParm const& prob_parm,
PmfData const * pmf_data)
pele::physics::PMF::PmfData::DataContainer const * pmf_data)
{
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* prob_hi = geomdata.ProbHi();
Expand Down Expand Up @@ -58,7 +56,7 @@ void pelelm_initdata(int i, int j, int k,
y1 = (y - prob_parm.standoff - 0.5*dx[1] + pert)*100;
y2 = (y - prob_parm.standoff + 0.5*dx[1] + pert)*100;

PMF::pmf(pmf_data,y1, y2, pmf_vals);
pele::physics::PMF::pmf(pmf_data,y1, y2, pmf_vals);

state(i,j,k,TEMP) = pmf_vals[0];;

Expand Down Expand Up @@ -88,9 +86,9 @@ void pelelm_initdata(int i, int j, int k,
state(i,j,k,FIRSTSPEC+n) = massfrac[n] * state(i,j,k,DENSITY);
}

amrex::Real yc = 0.5 * (prob_hi[1] + prob_lo[1]);
nE(i,j,k) = 2.0e13 * std::exp(-(y-yc)*(y-yc)/(Ly*Ly/100.0));
phiV(i,j,k) = 0.0;
amrex::Real yc = 0.5 * (prob_hi[1] + prob_lo[1]) + pert;
state(i,j,k,NE) = 2.0e13 * std::exp(-(y-yc)*(y-yc)/(Ly*Ly/500.0));
state(i,j,k,PHIV) = 0.0;
}

AMREX_GPU_DEVICE
Expand All @@ -105,7 +103,7 @@ bcnormal(
const amrex::Real time,
amrex::GeometryData const& geomdata,
ProbParm const& prob_parm,
PmfData const *pmf_data)
pele::physics::PMF::PmfData::DataContainer const * pmf_data)
{
const amrex::Real* prob_lo = geomdata.ProbLo();
amrex::GpuArray<amrex::Real, NUM_SPECIES + 4> pmf_vals = {0.0};
Expand All @@ -114,7 +112,7 @@ bcnormal(

auto eos = pele::physics::PhysicsType::eos();
if (sgn == 1) {
PMF::pmf(pmf_data,prob_lo[idir], prob_lo[idir], pmf_vals);
pele::physics::PMF::pmf(pmf_data,prob_lo[idir], prob_lo[idir], pmf_vals);

AMREX_D_TERM(s_ext[VELX] = 0.0;,
s_ext[VELY] = pmf_vals[1]*1e-2;,
Expand Down
Loading

0 comments on commit ee45a90

Please sign in to comment.