From cc8bfea3ca14c598ef08e8877a793b566f2f7054 Mon Sep 17 00:00:00 2001 From: Lucas Esclapez <13371051+esclapez@users.noreply.github.com> Date: Fri, 1 Jul 2022 08:14:07 -0700 Subject: [PATCH] Minor updates (#105) * Add an enstrophy derived (not EB aware yet) and get enstrophy integral in temporals. * Pass reaction data to derived to accomodate more cases. Really need to update the entire process. * Update EB_Pipe to accomodate turbulence in between plates (SINTEF). * Add channel flow input file. * Add dt to temporal data, always useful ... * Add new derived to doc. * Add min/max of state components to the temporals. * Add a description of the temporal diagnostics in the doc. * Update key for mass_balance to make it consistent with others. * Update all the input files to reflect changes in temporals. * Minor summer cleaning of trailing whitespaces. --- Docs/source/LMeXControls.rst | 29 ++- Exec/Cases/ChallengeProblem/input.3d-regt | 3 +- Exec/Cases/ChallengeProblem/input.3d-regt_AMR | 3 +- .../ChallengeProblem/input.3d-regt_Quarter | 3 +- Exec/Cases/CounterFlow/input.2d-regt | 3 +- Exec/Cases/CounterFlowSpray/input.2d-regt | 3 +- Exec/RegTests/EB_EnclosedFlame/input.2d-regt | 3 +- .../EB_EnclosedFlame/input.2d-regt_Hypre | 3 +- Exec/RegTests/EB_EnclosedFlame/input.3d-regt | 3 +- Exec/RegTests/EB_EnclosedVortex/input.2d-regt | 3 +- .../EB_EnclosedVortex/input.2d-regt_Hypre | 3 +- Exec/RegTests/EB_EnclosedVortex/input.3d-regt | 3 +- .../EB_PipeFlow/input.3d-channel_noEB | 68 +++++++ Exec/RegTests/EB_PipeFlow/pelelm_prob.H | 178 +++++++++++------- Exec/RegTests/EB_PipeFlow/pelelm_prob.cpp | 5 + Exec/RegTests/EB_PipeFlow/pelelm_prob_parm.H | 1 + Exec/RegTests/EnclosedFlame/input.2d-regt | 3 +- Exec/RegTests/EnclosedFlame/input.3d-regt | 3 +- Exec/RegTests/EnclosedInjection/input.2d-regt | 3 +- Exec/RegTests/FlameSheet/input.2d-regt | 5 +- Exec/RegTests/FlameSheet/input.3d-regt | 3 +- Exec/RegTests/FlameSheet/inputs.3d_Dodecane | 4 +- .../RegTests/FlameSheet/inputs.3d_DodecaneQSS | 4 +- Exec/RegTests/SprayTest/input.2d | 3 +- Exec/RegTests/TurbInflow/input.3d | 3 +- Exec/RegTests/TurbInflow/input.3d_BoxLoZ | 3 +- Exec/RegTests/TurbInflow/input.3d_negY | 3 +- Exec/RegTests/TurbInflow/input.3d_posX | 3 +- Exec/RegTests/TurbInflow/input.3d_twoInjs | 3 +- Source/DiffusionOp.cpp | 6 +- Source/PeleLM.H | 2 + Source/PeleLM.cpp | 6 +- Source/PeleLMAdvance.cpp | 4 +- Source/PeleLMAdvection.cpp | 22 +-- Source/PeleLMBC.cpp | 16 +- Source/PeleLMData.cpp | 10 +- Source/PeleLMDerive.H | 7 +- Source/PeleLMDerive.cpp | 4 +- Source/PeleLMDeriveFunc.H | 32 ++-- Source/PeleLMDeriveFunc.cpp | 134 +++++++++++-- Source/PeleLMDiffusion.cpp | 12 +- Source/PeleLMEB.cpp | 12 +- Source/PeleLMEos.cpp | 12 +- Source/PeleLMEvaluate.cpp | 8 +- Source/PeleLMEvolve.cpp | 6 +- Source/PeleLMForces.cpp | 16 +- Source/PeleLMInit.cpp | 4 +- Source/PeleLMPlot.cpp | 36 ++-- Source/PeleLMReactions.cpp | 14 +- Source/PeleLMSetup.cpp | 31 +-- Source/PeleLMTagging.cpp | 6 +- Source/PeleLMTemporals.cpp | 65 +++++-- Source/PeleLMTimestep.cpp | 2 +- Source/PeleLMTransportProp.cpp | 6 +- Source/PeleLMUMac.cpp | 16 +- Source/PeleLMUtils.H | 2 +- Source/PeleLMUtils.cpp | 78 ++++---- Source/PeleLM_Index.H | 6 +- Source/main.cpp | 4 +- 59 files changed, 603 insertions(+), 333 deletions(-) create mode 100644 Exec/RegTests/EB_PipeFlow/input.3d-channel_noEB diff --git a/Docs/source/LMeXControls.rst b/Docs/source/LMeXControls.rst index a7456d5a..f07deaad 100644 --- a/Docs/source/LMeXControls.rst +++ b/Docs/source/LMeXControls.rst @@ -127,7 +127,13 @@ The following list of derived variables are available in PeleLMeX: - Vorticity (2D) or vorticity magnitude (3D) * - `kinetic_energy` - 1 - - Kinetic energy + - Kinetic energy: 0.5 * rho * (u^2+v^2+w^2) + * - `enstrophy` + - 1 + - enstrophy: 0.5 * rho * (\omega_x^2+\omega_y^2+\omega_z^2) + * - `HeatRelease` + - 1 + - Heat release rate from chem. reactions Note that `mixture_fraction` and `progress_variable` requires additional inputs from the users as described below. @@ -203,6 +209,27 @@ Linear solvers are a key component of PeleLMeX algorithm, separate controls are Run-time diagnostics -------------------- +PeleLMeX provides a few diagnostics to check you simulations while it is running as well as adding basic analysis ingredients. + +It is often usefull to have an estimate of integrated quantities (kinetic energy, heat release rate, ,..), state extremas +or other overall balance information to get a sense of the status and sanity of the simulation. To this end, it is possible +to activate `temporal` diagnostics performing these reductions at given intervals: + +:: + + #-------------------------TEMPORALS--------------------------- + peleLM.do_temporals = 1 # [OPT, DEF=0] Activate temporal diagnostics + peleLM.temporal_int = 10 # [OPT, DEF=5] Temporal freq. + peleLM.do_extremas = 1 # [OPT, DEF=0] Trigger extremas, if temporals activated + peleLM.do_mass_balance = 1 # [OPT, DEF=0] Compute mass balance, if temporals activated + +The `do_temporal` flag will trigger the creation of a `temporals` folder in your run directory and the following entries +will be appended to an ASCII `temporals/tempState` file: step, time, dt, kin. energy integral, enstrophy integral, mean pressure +, fuel consumption rate integral, heat release rate integral. Additionnally, if the `do_temporal` flag is activated, one can +turn on state extremas (stored in `temporals/tempExtremas` as min/max for each state entry) and mass balance (stored in +`temporals/tempMass`) computing the total mass, dMdt and mass fluxes across the domain boundary as well as the error in +the balance (dMdt - sum of fluxes). + Combustion diagnostics often involve the use of a mixture fraction and/or a progress variable, both of which can be defined at run time and added to the derived variables included in the plotfile. If `mixture_fraction` or `progress_variable` is added to the `amr.derive_plot_vars` list, one need to provide input for defining those. The mixture fraction is based on diff --git a/Exec/Cases/ChallengeProblem/input.3d-regt b/Exec/Cases/ChallengeProblem/input.3d-regt index ac59debd..1c9c638c 100644 --- a/Exec/Cases/ChallengeProblem/input.3d-regt +++ b/Exec/Cases/ChallengeProblem/input.3d-regt @@ -106,8 +106,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 peleLM.num_init_iter = 1 peleLM.num_divu_iter = 1 diff --git a/Exec/Cases/ChallengeProblem/input.3d-regt_AMR b/Exec/Cases/ChallengeProblem/input.3d-regt_AMR index 6b5bedf1..b5e90fb3 100644 --- a/Exec/Cases/ChallengeProblem/input.3d-regt_AMR +++ b/Exec/Cases/ChallengeProblem/input.3d-regt_AMR @@ -106,8 +106,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 peleLM.num_init_iter = 1 peleLM.num_divu_iter = 1 diff --git a/Exec/Cases/ChallengeProblem/input.3d-regt_Quarter b/Exec/Cases/ChallengeProblem/input.3d-regt_Quarter index 3d72c52a..28e97f7c 100644 --- a/Exec/Cases/ChallengeProblem/input.3d-regt_Quarter +++ b/Exec/Cases/ChallengeProblem/input.3d-regt_Quarter @@ -65,8 +65,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 peleLM.num_init_iter = 1 peleLM.num_divu_iter = 1 diff --git a/Exec/Cases/CounterFlow/input.2d-regt b/Exec/Cases/CounterFlow/input.2d-regt index 3b092a6c..ab3ac493 100644 --- a/Exec/Cases/CounterFlow/input.2d-regt +++ b/Exec/Cases/CounterFlow/input.2d-regt @@ -72,8 +72,7 @@ mac_proj.verbose = 0 #diffusion.verbose = 2 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #--------------------REFINEMENT CONTROL------------------------ amr.refinement_indicators = temp diff --git a/Exec/Cases/CounterFlowSpray/input.2d-regt b/Exec/Cases/CounterFlowSpray/input.2d-regt index a5e8a925..477f5386 100644 --- a/Exec/Cases/CounterFlowSpray/input.2d-regt +++ b/Exec/Cases/CounterFlowSpray/input.2d-regt @@ -110,8 +110,7 @@ nodal_proj.rtol = 1.0e-10 #diffusion.verbose = 2 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #--------------------REFINEMENT CONTROL------------------------ amr.refinement_indicators = temp diff --git a/Exec/RegTests/EB_EnclosedFlame/input.2d-regt b/Exec/RegTests/EB_EnclosedFlame/input.2d-regt index 9d03f8f7..d76d5577 100644 --- a/Exec/RegTests/EB_EnclosedFlame/input.2d-regt +++ b/Exec/RegTests/EB_EnclosedFlame/input.2d-regt @@ -41,8 +41,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 0 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00700 amr.check_int = 20 diff --git a/Exec/RegTests/EB_EnclosedFlame/input.2d-regt_Hypre b/Exec/RegTests/EB_EnclosedFlame/input.2d-regt_Hypre index 252d912d..29ab051a 100644 --- a/Exec/RegTests/EB_EnclosedFlame/input.2d-regt_Hypre +++ b/Exec/RegTests/EB_EnclosedFlame/input.2d-regt_Hypre @@ -41,8 +41,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 0 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00030 amr.check_int = 20 diff --git a/Exec/RegTests/EB_EnclosedFlame/input.3d-regt b/Exec/RegTests/EB_EnclosedFlame/input.3d-regt index 130194dd..f660252a 100644 --- a/Exec/RegTests/EB_EnclosedFlame/input.3d-regt +++ b/Exec/RegTests/EB_EnclosedFlame/input.3d-regt @@ -41,8 +41,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 1 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00700 amr.check_int = 20 diff --git a/Exec/RegTests/EB_EnclosedVortex/input.2d-regt b/Exec/RegTests/EB_EnclosedVortex/input.2d-regt index 77473af5..def0a630 100644 --- a/Exec/RegTests/EB_EnclosedVortex/input.2d-regt +++ b/Exec/RegTests/EB_EnclosedVortex/input.2d-regt @@ -39,8 +39,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.check_int = 100 amr.plot_int = 50 diff --git a/Exec/RegTests/EB_EnclosedVortex/input.2d-regt_Hypre b/Exec/RegTests/EB_EnclosedVortex/input.2d-regt_Hypre index 212d7c9a..06ae1f05 100644 --- a/Exec/RegTests/EB_EnclosedVortex/input.2d-regt_Hypre +++ b/Exec/RegTests/EB_EnclosedVortex/input.2d-regt_Hypre @@ -40,8 +40,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00005 amr.check_int = 5 diff --git a/Exec/RegTests/EB_EnclosedVortex/input.3d-regt b/Exec/RegTests/EB_EnclosedVortex/input.3d-regt index 5683d3be..f17f4f3f 100644 --- a/Exec/RegTests/EB_EnclosedVortex/input.3d-regt +++ b/Exec/RegTests/EB_EnclosedVortex/input.3d-regt @@ -40,8 +40,7 @@ peleLM.num_divu_iter = 0 peleLM.num_init_iter = 1 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk01000 amr.check_int = 100 diff --git a/Exec/RegTests/EB_PipeFlow/input.3d-channel_noEB b/Exec/RegTests/EB_PipeFlow/input.3d-channel_noEB new file mode 100644 index 00000000..9e345cce --- /dev/null +++ b/Exec/RegTests/EB_PipeFlow/input.3d-channel_noEB @@ -0,0 +1,68 @@ +#----------------------DOMAIN DEFINITION------------------------ +# https://amrex-combustion.github.io/PeleLMeX/LMeXControls.html +geometry.is_periodic = 1 0 1 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = 0.0 -0.005 -0.0075 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.03 0.005 0.0075 # x_hi y_hi (z_hi) + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +# STM: NoSlipWallIsotherm crashes? +peleLM.lo_bc = Interior NoSlipWallAdiab Interior +peleLM.hi_bc = Interior NoSlipWallAdiab Interior + + +#-------------------------AMR CONTROL---------------------------- +#amr.n_cell = 96 32 48 # Level 0 number of cells in each direction +#amr.n_cell = 1152 384 576 # Level 0 number of cells in each direction +amr.n_cell = 384 128 192 +amr.n_cell = 576 192 256 +amr.v = 1 # AMR verbose +amr.max_level = 1 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 1000 # how often to regrid +amr.n_error_buf = 2 2 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.T_mean = 750.0 +prob.P_mean = 101325.0 +prob.problem_type = 2 +prob.meanFlowMag = 5.0 + +#-------------------------PeleLM CONTROL---------------------------- +peleLM.v = 1 +peleLM.incompressible = 0 +peleLM.gradP0 = -709.79 0.0 0.0 + +peleLM.do_temporals = 1 +peleLM.temporal_int = 10 + +#amr.restart = chk30000 +amr.max_step = 1000 +amr.stop_time = 0.2 +amr.plot_int = 100 +amr.check_int = 100 +amr.cfl = 0.5 +amr.derive_plot_vars = avg_pressure mag_vort + +peleLM.chem_integrator = "ReactorNull" # Chemistry integrator, from PelePhysics available list + +#--------------------REFINEMENT CONTROL------------------------ +#amr.refinement_indicators = temp +#amr.temp.max_level = 2 +#amr.temp.value_greater = 305 +#amr.temp.field_name = temp + +amr.refinement_indicators = yLow yHigh +amr.yLow.in_box_lo = -0.001 -0.0052 -0.0085 +amr.yLow.in_box_hi = 0.031 -0.0045 0.0085 +amr.yHigh.in_box_lo = -0.001 0.0045 -0.0085 +amr.yHigh.in_box_hi = 0.031 0.0052 0.0085 + +#amrex.fpe_trap_invalid = 1 +#amrex.fpe_trap_zero = 1 +#amrex.fpe_trap_overflow = 1 diff --git a/Exec/RegTests/EB_PipeFlow/pelelm_prob.H b/Exec/RegTests/EB_PipeFlow/pelelm_prob.H index 85564838..8689cce6 100644 --- a/Exec/RegTests/EB_PipeFlow/pelelm_prob.H +++ b/Exec/RegTests/EB_PipeFlow/pelelm_prob.H @@ -43,37 +43,75 @@ void pelelm_initdata(int i, int j, int k, massfrac[O2_ID] = 0.233; massfrac[N2_ID] = 0.767; - switch(prob_parm.meanFlowDir) { - case 1 : - AMREX_D_TERM(state(i,j,k,VELX) = prob_parm.meanFlowMag;, - state(i,j,k,VELY) = 5.0 * std::sin(2 * Pi * 4 * y / Ly) * std::sin(2 * Pi * 11 * x / Lx);, - state(i,j,k,VELZ) = 5.0 * std::sin(2 * Pi * 5 * z / Lz) * std::sin(2 * Pi * 4 * x / Lx)); - break; - case -1 : - AMREX_D_TERM(state(i,j,k,VELX) = -prob_parm.meanFlowMag;, - state(i,j,k,VELY) = 0.0;, - state(i,j,k,VELZ) = 0.0); - break; - case 2 : - AMREX_D_TERM(state(i,j,k,VELX) = 0.0;, - state(i,j,k,VELY) = prob_parm.meanFlowMag;, - state(i,j,k,VELZ) = 0.0); - break; - case -2 : - AMREX_D_TERM(state(i,j,k,VELX) = 0.0;, - state(i,j,k,VELY) = -prob_parm.meanFlowMag;, - state(i,j,k,VELZ) = 0.0); - break; - case 3 : - AMREX_D_TERM(state(i,j,k,VELX) = prob_parm.meanFlowMag;, - state(i,j,k,VELY) = prob_parm.meanFlowMag;, - state(i,j,k,VELZ) = 0.0); - break; - case -3 : - AMREX_D_TERM(state(i,j,k,VELX) = -prob_parm.meanFlowMag;, - state(i,j,k,VELY) = -prob_parm.meanFlowMag;, - state(i,j,k,VELZ) = 0.0); - break; + if ( prob_parm.flowType == 1 ) { + switch(prob_parm.meanFlowDir) { + case 1 : + AMREX_D_TERM(state(i,j,k,VELX) = prob_parm.meanFlowMag;, + state(i,j,k,VELY) = 5.0 * std::sin(2 * Pi * 4 * y / Ly) * std::sin(2 * Pi * 11 * x / Lx);, + state(i,j,k,VELZ) = 5.0 * std::sin(2 * Pi * 5 * z / Lz) * std::sin(2 * Pi * 4 * x / Lx)); + break; + case -1 : + AMREX_D_TERM(state(i,j,k,VELX) = -prob_parm.meanFlowMag;, + state(i,j,k,VELY) = 0.0;, + state(i,j,k,VELZ) = 0.0); + break; + case 2 : + AMREX_D_TERM(state(i,j,k,VELX) = 0.0;, + state(i,j,k,VELY) = prob_parm.meanFlowMag;, + state(i,j,k,VELZ) = 0.0); + break; + case -2 : + AMREX_D_TERM(state(i,j,k,VELX) = 0.0;, + state(i,j,k,VELY) = -prob_parm.meanFlowMag;, + state(i,j,k,VELZ) = 0.0); + break; + case 3 : + AMREX_D_TERM(state(i,j,k,VELX) = prob_parm.meanFlowMag;, + state(i,j,k,VELY) = prob_parm.meanFlowMag;, + state(i,j,k,VELZ) = 0.0); + break; + case -3 : + AMREX_D_TERM(state(i,j,k,VELX) = -prob_parm.meanFlowMag;, + state(i,j,k,VELY) = -prob_parm.meanFlowMag;, + state(i,j,k,VELZ) = 0.0); + break; + } + } else if ( prob_parm.flowType == 2 ) { + const amrex::Real dpdx = 709.79; // Pa/m + const amrex::Real c_f = 0.0062418; + amrex::GpuArray u; + u[0] = sqrt(dpdx*Ly/(0.4688*c_f))*(1.0-pow(fabs(2.0*y/Ly),10.0)); + + // Can't get AMReX random engine in here without changing + // the initdata function signature -> just use not random IC + + //amrex::Real randp = 0.0; + //amrex::Real randr = 0.0; + //amrex::Real noisefac = 0.0; + // srand(time(NULL) + getpid() + getppid()); + //randp = amrex::Random()/(amrex::Real)(RAND_MAX/1.0); + //randr = amrex::Random()/(amrex::Real)(RAND_MAX/1.0); + //// x direction + //noisefac = sqrt(-2.0*log(randr))*sin(2.0*Pi*randp); + //u[0] += prob_parm.meanFlowMag * noisefac; + //// y direction. Note the cos(). + //noisefac = sqrt(-2.0*log(randr))*cos(2.0*Pi*randp); + //u[1] += prob_parm.meanFlowMag * noisefac; + //// z direction + //randp = (amrex::Real)rand()/(amrex::Real)(RAND_MAX/1.0); + //randr = (amrex::Real)rand()/(amrex::Real)(RAND_MAX/1.0); + //noisefac = sqrt(-2.0*log(randr))*sin(2.0*Pi*randp); + //u[2] += prob_parm.meanFlowMag * noisefac; + + //AMREX_D_TERM(state(i,j,k,VELX) = u[0];, + // state(i,j,k,VELY) = u[1];, + // state(i,j,k,VELZ) = u[2]); + AMREX_D_TERM(state(i,j,k,VELX) = u[0] * (1.0 + std::sin(2 * Pi * 5 * y / Ly) * std::sin(2 * Pi * 14 * z / Lz) + + std::sin(2 * Pi * 7 * y / Ly) * std::sin(2 * Pi * 6 * (z-0.00541) / Lz));, + state(i,j,k,VELY) = 0.2*u[0] * ( std::sin(2 * Pi * 5 * y / Ly) * std::sin(2 * Pi * 7 * (x-0.0013) / Lx) + + std::sin(2 * Pi * 6 * y / Ly) * std::sin(2 * Pi * 13 * x / Lx));, + state(i,j,k,VELZ) = 0.2*u[0] * ( std::sin(2 * Pi * 11 * z / Lz) * std::sin(2 * Pi * 4 * x / Lx) + + std::sin(2 * Pi * 3 * z / Lz) * std::sin(2 * Pi * 9 * (x-0.00779) / Lx))); } if (is_incompressible) return; @@ -113,46 +151,48 @@ bcnormal( ProbParm const& prob_parm, pele::physics::PMF::PmfData::DataContainer const * /*pmf_data*/) { - const amrex::Real* prob_lo = geomdata.ProbLo(); - //amrex::GpuArray pmf_vals = {0.0}; - amrex::Real massfrac[NUM_SPECIES] = {0.0}; - - auto eos = pele::physics::PhysicsType::eos(); - - switch(prob_parm.meanFlowDir) { - case 1 : - AMREX_D_TERM(s_ext[VELX] = prob_parm.meanFlowMag;, - s_ext[VELY] = 0.0;, - s_ext[VELZ] = 0.0); - break; - case 2 : - AMREX_D_TERM(s_ext[VELY] = prob_parm.meanFlowMag;, - s_ext[VELX] = 0.0;, - s_ext[VELZ] = 0.0); - break; - case 3 : - AMREX_D_TERM(s_ext[VELZ] = prob_parm.meanFlowMag;, - s_ext[VELX] = 0.0;, - s_ext[VELY] = 0.0); - break; - } - - massfrac[O2_ID] = 0.333; - massfrac[N2_ID] = 0.667; - - s_ext[TEMP] = prob_parm.T_mean; - - amrex::Real rho_cgs, P_cgs, RhoH_temp; - P_cgs = prob_parm.P_mean * 10.0; - - eos.PYT2R(P_cgs, massfrac, s_ext[TEMP], rho_cgs); - s_ext[DENSITY] = rho_cgs * 1.0e3; + const amrex::Real* prob_lo = geomdata.ProbLo(); + //amrex::GpuArray pmf_vals = {0.0}; + amrex::Real massfrac[NUM_SPECIES] = {0.0}; - eos.TY2H(s_ext[TEMP], massfrac, RhoH_temp); - s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY]; // CGS -> MKS conversion + auto eos = pele::physics::PhysicsType::eos(); - for (int n = 0; n < NUM_SPECIES; n++) { - s_ext[FIRSTSPEC+n] = massfrac[n] * s_ext[DENSITY]; + if ( prob_parm.flowType == 1 ) { + switch(prob_parm.meanFlowDir) { + case 1 : + AMREX_D_TERM(s_ext[VELX] = prob_parm.meanFlowMag;, + s_ext[VELY] = 0.0;, + s_ext[VELZ] = 0.0); + break; + case 2 : + AMREX_D_TERM(s_ext[VELY] = prob_parm.meanFlowMag;, + s_ext[VELX] = 0.0;, + s_ext[VELZ] = 0.0); + break; + case 3 : + AMREX_D_TERM(s_ext[VELZ] = prob_parm.meanFlowMag;, + s_ext[VELX] = 0.0;, + s_ext[VELY] = 0.0); + break; + } + + massfrac[O2_ID] = 0.333; + massfrac[N2_ID] = 0.667; + + s_ext[TEMP] = prob_parm.T_mean; + + amrex::Real rho_cgs, P_cgs, RhoH_temp; + P_cgs = prob_parm.P_mean * 10.0; + + eos.PYT2R(P_cgs, massfrac, s_ext[TEMP], rho_cgs); + s_ext[DENSITY] = rho_cgs * 1.0e3; + + eos.TY2H(s_ext[TEMP], massfrac, RhoH_temp); + s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY]; // CGS -> MKS conversion + + for (int n = 0; n < NUM_SPECIES; n++) { + s_ext[FIRSTSPEC+n] = massfrac[n] * s_ext[DENSITY]; + } } } diff --git a/Exec/RegTests/EB_PipeFlow/pelelm_prob.cpp b/Exec/RegTests/EB_PipeFlow/pelelm_prob.cpp index 869c6b7b..27c32549 100644 --- a/Exec/RegTests/EB_PipeFlow/pelelm_prob.cpp +++ b/Exec/RegTests/EB_PipeFlow/pelelm_prob.cpp @@ -9,7 +9,11 @@ void PeleLM::readProbParm() pp.query("P_mean", prob_parm->P_mean); pp.query("meanFlowDir", prob_parm->meanFlowDir); pp.query("meanFlowMag", prob_parm->meanFlowMag); + pp.query("problem_type",prob_parm->flowType); + AMREX_ALWAYS_ASSERT(prob_parm->flowType == 1 || + prob_parm->flowType == 2); + /* if (!m_incompressible) { auto& trans_parm = PeleLM::trans_parms.host_trans_parm(); amrex::ParmParse pptr("transport"); @@ -19,4 +23,5 @@ void PeleLM::readProbParm() pp.query("const_diffusivity", trans_parm.const_diffusivity); PeleLM::trans_parms.sync_to_device(); } + */ } diff --git a/Exec/RegTests/EB_PipeFlow/pelelm_prob_parm.H b/Exec/RegTests/EB_PipeFlow/pelelm_prob_parm.H index 614c3130..54c664a5 100644 --- a/Exec/RegTests/EB_PipeFlow/pelelm_prob_parm.H +++ b/Exec/RegTests/EB_PipeFlow/pelelm_prob_parm.H @@ -13,5 +13,6 @@ struct ProbParm amrex::Real P_mean = 101325.0_rt; amrex::Real meanFlowMag = 5.0; int meanFlowDir = 1; + int flowType = 1; }; #endif diff --git a/Exec/RegTests/EnclosedFlame/input.2d-regt b/Exec/RegTests/EnclosedFlame/input.2d-regt index 1d6c4b0b..e42073c4 100644 --- a/Exec/RegTests/EnclosedFlame/input.2d-regt +++ b/Exec/RegTests/EnclosedFlame/input.2d-regt @@ -42,8 +42,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00005 amr.check_int = 100 diff --git a/Exec/RegTests/EnclosedFlame/input.3d-regt b/Exec/RegTests/EnclosedFlame/input.3d-regt index 7f569c9d..dce43e6e 100644 --- a/Exec/RegTests/EnclosedFlame/input.3d-regt +++ b/Exec/RegTests/EnclosedFlame/input.3d-regt @@ -41,8 +41,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00005 amr.check_int = 5 diff --git a/Exec/RegTests/EnclosedInjection/input.2d-regt b/Exec/RegTests/EnclosedInjection/input.2d-regt index d95d599a..aa6c6c3c 100644 --- a/Exec/RegTests/EnclosedInjection/input.2d-regt +++ b/Exec/RegTests/EnclosedInjection/input.2d-regt @@ -51,8 +51,7 @@ peleLM.typical_values_reset_int = 10 # Advection # Temporals peleLM.do_temporals = 1 -peleLM.temporal_int = 1 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 # Chamber peleLM.closed_chamber = 1 diff --git a/Exec/RegTests/FlameSheet/input.2d-regt b/Exec/RegTests/FlameSheet/input.2d-regt index 729b85ba..52a11f3f 100644 --- a/Exec/RegTests/FlameSheet/input.2d-regt +++ b/Exec/RegTests/FlameSheet/input.2d-regt @@ -39,8 +39,9 @@ peleLM.sdc_iterMax = 2 peleLM.floor_species = 0 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.temporal_int = 10 +peleLM.do_extremas = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00005 #amr.check_int = 2000 diff --git a/Exec/RegTests/FlameSheet/input.3d-regt b/Exec/RegTests/FlameSheet/input.3d-regt index fae57de1..a57d5011 100644 --- a/Exec/RegTests/FlameSheet/input.3d-regt +++ b/Exec/RegTests/FlameSheet/input.3d-regt @@ -39,8 +39,7 @@ peleLM.sdc_iterMax = 2 peleLM.floor_species = 0 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 #amr.restart = chk00005 #amr.check_int = 2000 diff --git a/Exec/RegTests/FlameSheet/inputs.3d_Dodecane b/Exec/RegTests/FlameSheet/inputs.3d_Dodecane index 06dcb3b9..282553fa 100644 --- a/Exec/RegTests/FlameSheet/inputs.3d_Dodecane +++ b/Exec/RegTests/FlameSheet/inputs.3d_Dodecane @@ -38,9 +38,7 @@ peleLM.use_wbar = 1 peleLM.sdc_iterMax = 2 peleLM.floor_species = 0 -peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_temporals = 0 #amr.restart = chk00005 #amr.check_int = 2000 diff --git a/Exec/RegTests/FlameSheet/inputs.3d_DodecaneQSS b/Exec/RegTests/FlameSheet/inputs.3d_DodecaneQSS index bff1a766..446a1921 100644 --- a/Exec/RegTests/FlameSheet/inputs.3d_DodecaneQSS +++ b/Exec/RegTests/FlameSheet/inputs.3d_DodecaneQSS @@ -38,9 +38,7 @@ peleLM.use_wbar = 1 peleLM.sdc_iterMax = 2 peleLM.floor_species = 0 -peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_temporals = 0 #amr.restart = chk00005 #amr.check_int = 2000 diff --git a/Exec/RegTests/SprayTest/input.2d b/Exec/RegTests/SprayTest/input.2d index ca28a993..7177fcba 100644 --- a/Exec/RegTests/SprayTest/input.2d +++ b/Exec/RegTests/SprayTest/input.2d @@ -48,8 +48,7 @@ peleLM.sdc_iterMax = 1 peleLM.floor_species = 0 peleLM.do_temporals = 0 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 0 peleLM.num_init_iter = 0 amr.plot_int = 100 diff --git a/Exec/RegTests/TurbInflow/input.3d b/Exec/RegTests/TurbInflow/input.3d index 962ef7ba..fb491395 100644 --- a/Exec/RegTests/TurbInflow/input.3d +++ b/Exec/RegTests/TurbInflow/input.3d @@ -57,8 +57,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 1 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.plot_int = 10 amr.max_step = 300 diff --git a/Exec/RegTests/TurbInflow/input.3d_BoxLoZ b/Exec/RegTests/TurbInflow/input.3d_BoxLoZ index 7de53b43..ba355bf2 100644 --- a/Exec/RegTests/TurbInflow/input.3d_BoxLoZ +++ b/Exec/RegTests/TurbInflow/input.3d_BoxLoZ @@ -56,8 +56,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.plot_int = 10 amr.max_step = 300 diff --git a/Exec/RegTests/TurbInflow/input.3d_negY b/Exec/RegTests/TurbInflow/input.3d_negY index 5b1a5eba..c2dea9d6 100644 --- a/Exec/RegTests/TurbInflow/input.3d_negY +++ b/Exec/RegTests/TurbInflow/input.3d_negY @@ -56,8 +56,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.plot_int = 10 amr.max_step = 100 diff --git a/Exec/RegTests/TurbInflow/input.3d_posX b/Exec/RegTests/TurbInflow/input.3d_posX index 58d66469..e3cda3ab 100644 --- a/Exec/RegTests/TurbInflow/input.3d_posX +++ b/Exec/RegTests/TurbInflow/input.3d_posX @@ -56,8 +56,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.plot_int = 10 amr.max_step = 100 diff --git a/Exec/RegTests/TurbInflow/input.3d_twoInjs b/Exec/RegTests/TurbInflow/input.3d_twoInjs index dcac28f7..5ff28f7e 100644 --- a/Exec/RegTests/TurbInflow/input.3d_twoInjs +++ b/Exec/RegTests/TurbInflow/input.3d_twoInjs @@ -69,8 +69,7 @@ peleLM.num_divu_iter = 1 peleLM.num_init_iter = 3 peleLM.do_temporals = 1 -peleLM.temporal_int = 2 -peleLM.mass_balance = 1 +peleLM.do_mass_balance = 1 amr.plot_int = 10 amr.max_step = 100 diff --git a/Source/DiffusionOp.cpp b/Source/DiffusionOp.cpp index e8ce571d..a1e1f11d 100644 --- a/Source/DiffusionOp.cpp +++ b/Source/DiffusionOp.cpp @@ -193,7 +193,7 @@ void DiffusionOp::diffuse_scalar(Vector const& a_phi, int phi_comp, if (have_bcoeff) { int doZeroVisc = 1; Vector subBCRec = {a_bcrec.begin()+comp,a_bcrec.begin()+comp+m_ncomp}; - Array bcoeff_ec = m_pelelm->getDiffusivity(lev, bcoeff_comp+comp, m_ncomp, + Array bcoeff_ec = m_pelelm->getDiffusivity(lev, bcoeff_comp+comp, m_ncomp, doZeroVisc, subBCRec, *a_bcoeff[lev]); #ifdef AMREX_USE_EB m_scal_solve_op->setBCoeffs(lev, GetArrOfConstPtrs(bcoeff_ec), MLMG::Location::FaceCentroid); @@ -325,11 +325,11 @@ void DiffusionOp::computeDiffLap(Vector const& a_laps, int lap_comp, m_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(bcoeff_ec)); #endif m_scal_apply_op->setLevelBC(lev, &component[lev]); - } + } MLMG mlmg(*m_scal_apply_op); mlmg.apply(GetVecOfPtrs(laps), GetVecOfPtrs(component)); - } + } } void DiffusionOp::computeDiffFluxes(Vector> const& a_flux, int flux_comp, diff --git a/Source/PeleLM.H b/Source/PeleLM.H index 69462b4d..32c21a69 100644 --- a/Source/PeleLM.H +++ b/Source/PeleLM.H @@ -1181,6 +1181,7 @@ class PeleLM : public amrex::AmrCore { // Temporals int m_do_temporals = 0; int m_temp_int = 5; + int m_do_extremas = 0; int m_do_massBalance = 0; int m_do_energyBalance = 0; int m_do_speciesBalance = 0; @@ -1196,6 +1197,7 @@ class PeleLM : public amrex::AmrCore { amrex::Array m_domainUmacFlux; std::ofstream tmpStateFile; + std::ofstream tmpExtremasFile; std::ofstream tmpMassFile; std::ofstream tmpSpecFile; diff --git a/Source/PeleLM.cpp b/Source/PeleLM.cpp index 0577aab5..189d7a9f 100644 --- a/Source/PeleLM.cpp +++ b/Source/PeleLM.cpp @@ -13,11 +13,11 @@ PeleLM::~PeleLM() for (int lev = 0; lev <= finest_level; ++lev) { ClearLevel(lev); } - + if (!m_incompressible) { trans_parms.deallocate(); m_reactor->close(); - } + } closeTempFile(); typical_values.clear(); @@ -35,7 +35,7 @@ PeleLM::LevelData* PeleLM::getLevelDataPtr(int lev, const PeleLM::TimeStamp &a_time, int /*useUMac*/) { AMREX_ASSERT(a_time==AmrOldTime || a_time==AmrNewTime || a_time==AmrHalfTime); - if ( a_time == AmrOldTime ) { + if ( a_time == AmrOldTime ) { return m_leveldata_old[lev].get(); } else if ( a_time == AmrNewTime ) { return m_leveldata_new[lev].get(); diff --git a/Source/PeleLMAdvance.cpp b/Source/PeleLMAdvance.cpp index cd4c0e6b..9d93b7a6 100644 --- a/Source/PeleLMAdvance.cpp +++ b/Source/PeleLMAdvance.cpp @@ -271,11 +271,11 @@ void PeleLM::oneSDC(int sdcIter, averageDownnE(AmrNewTime); #endif fillPatchState(AmrNewTime); - + calcDiffusivity(AmrNewTime); computeDifferentialDiffusionTerms(AmrNewTime,diffData); if (m_has_divu) { - int is_initialization = 0; // Not here + int is_initialization = 0; // Not here int computeDiffusionTerm = 0; // Nope, we just did that int do_avgDown = 1; // Always calcDivU(is_initialization,computeDiffusionTerm,do_avgDown,AmrNewTime,diffData); diff --git a/Source/PeleLMAdvection.cpp b/Source/PeleLMAdvection.cpp index 0d611a75..ed1f9015 100644 --- a/Source/PeleLMAdvection.cpp +++ b/Source/PeleLMAdvection.cpp @@ -37,7 +37,7 @@ void PeleLM::computeVelocityAdvTerm(std::unique_ptr &advData) int add_gradP = 1; getVelForces(AmrOldTime,GetVecOfPtrs(divtau),GetVecOfPtrs(velForces),nGrow_force,add_gradP); - auto bcRecVel = fetchBCRecArray(VELX,AMREX_SPACEDIM); + auto bcRecVel = fetchBCRecArray(VELX,AMREX_SPACEDIM); auto bcRecVel_d = convertToDeviceVector(bcRecVel); auto AdvTypeVel = fetchAdvTypeArray(VELX,AMREX_SPACEDIM); auto AdvTypeVel_d = convertToDeviceVector(AdvTypeVel); @@ -102,11 +102,11 @@ void PeleLM::computeVelocityAdvTerm(std::unique_ptr &advData) m_advection_type); } #ifdef AMREX_USE_EB - EB_set_covered_faces(GetArrOfPtrs(fluxes[lev]),0.); - EB_set_covered_faces(GetArrOfPtrs(faces[lev]),0.); + EB_set_covered_faces(GetArrOfPtrs(fluxes[lev]),0.); + EB_set_covered_faces(GetArrOfPtrs(faces[lev]),0.); #endif } - + //---------------------------------------------------------------- // Average down fluxes to ensure C/F consistency for (int lev = finest_level; lev > 0; --lev) { @@ -283,11 +283,11 @@ void PeleLM::computeScalarAdvTerms(std::unique_ptr &advData) auto bcRecSpec_d = convertToDeviceVector(bcRecSpec); auto AdvTypeSpec = fetchAdvTypeArray(FIRSTSPEC,NUM_SPECIES); auto AdvTypeSpec_d = convertToDeviceVector(AdvTypeSpec); - auto bcRecTemp = fetchBCRecArray(TEMP,1); + auto bcRecTemp = fetchBCRecArray(TEMP,1); auto bcRecTemp_d = convertToDeviceVector(bcRecTemp); auto AdvTypeTemp = fetchAdvTypeArray(TEMP,1); auto AdvTypeTemp_d = convertToDeviceVector(AdvTypeTemp); - auto bcRecRhoH = fetchBCRecArray(RHOH,1); + auto bcRecRhoH = fetchBCRecArray(RHOH,1); auto bcRecRhoH_d = convertToDeviceVector(bcRecRhoH); auto AdvTypeRhoH = fetchAdvTypeArray(RHOH,1); auto AdvTypeRhoH_d = convertToDeviceVector(AdvTypeRhoH); @@ -304,7 +304,7 @@ void PeleLM::computeScalarAdvTerms(std::unique_ptr &advData) //---------------------------------------------------------------- // Loop over levels and get the fluxes - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) { // Get level data ptr Old auto ldata_p = getLevelDataPtr(lev,AmrOldTime); @@ -459,7 +459,7 @@ void PeleLM::computeScalarAdvTerms(std::unique_ptr &advData) { rho_ed(i,j,k) = 0.0; }); - } else if (flagfab.getType(ebx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(ebx) != FabType::regular ) { // EB containing boxes const auto& afrac = areafrac[idim]->array(mfi); amrex::ParallelFor(ebx, [rho_ed, rhoY_ed, afrac] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -548,7 +548,7 @@ void PeleLM::computeScalarAdvTerms(std::unique_ptr &advData) { rhoHm(i,j,k) = 0.0; }); - } else if (flagfab.getType(ebx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(ebx) != FabType::regular ) { // EB containing boxes const auto& afrac = areafrac[idim]->array(mfi); amrex::ParallelFor(ebx, [rho, rhoY, T, rhoHm, afrac] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -610,7 +610,7 @@ void PeleLM::computeScalarAdvTerms(std::unique_ptr &advData) m_advection_type); } #ifdef AMREX_USE_EB - EB_set_covered_faces(GetArrOfPtrs(fluxes[lev]),0.); + EB_set_covered_faces(GetArrOfPtrs(fluxes[lev]),0.); #endif } @@ -781,7 +781,7 @@ void PeleLM::computePassiveAdvTerms(std::unique_ptr &advData, //---------------------------------------------------------------- // Loop over levels and get the fluxes - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) { // Get level data ptr Old auto ldata_p = getLevelDataPtr(lev,AmrOldTime); diff --git a/Source/PeleLMBC.cpp b/Source/PeleLMBC.cpp index 1034337d..d0922296 100644 --- a/Source/PeleLMBC.cpp +++ b/Source/PeleLMBC.cpp @@ -85,14 +85,14 @@ InterpBase* PeleLM::getInterpolator() { // // Get EB-aware interpolater when needed // -#ifdef AMREX_USE_EB +#ifdef AMREX_USE_EB return (EBFactory(0).isAllRegular()) ? &mf_cell_cons_interp : &eb_mf_cell_cons_interp; #else return &mf_cell_cons_interp; #endif } - + void PeleLM::setBoundaryConditions() { // Initialize the BCRecs @@ -116,7 +116,7 @@ void PeleLM::setBoundaryConditions() { m_bcrec_state[VELX+idim].setHi(idim2,tang_vel_bc[hi_bc[idim2]]); } } - } + } // General forces: use int_dir in interior and foextrap otherwise for (int i = 0; i < sizeForceBC; i++) { @@ -387,7 +387,7 @@ void PeleLM::fillpatch_density(int lev, auto* mapper = getInterpolator(); // Density - PhysBCFunct> crse_bndry_func_rho(geom[lev-1], fetchBCRecArray(DENSITY,1), + PhysBCFunct> crse_bndry_func_rho(geom[lev-1], fetchBCRecArray(DENSITY,1), PeleLMCCFillExtDirDens{lprobparm, lpmfdata, m_nAux}); PhysBCFunct> fine_bndry_func_rho(geom[lev], fetchBCRecArray(DENSITY,1), PeleLMCCFillExtDirDens{lprobparm, lpmfdata, m_nAux}); @@ -559,7 +559,7 @@ void PeleLM::fillpatch_nE(int lev, auto* mapper = getInterpolator(); // nE - PhysBCFunct> crse_bndry_func(geom[lev-1], fetchBCRecArray(NE,1), + PhysBCFunct> crse_bndry_func(geom[lev-1], fetchBCRecArray(NE,1), PeleLMCCFillExtDirnE{lprobparm, lpmfdata, m_nAux}); PhysBCFunct> fine_bndry_func(geom[lev], fetchBCRecArray(NE,1), PeleLMCCFillExtDirnE{lprobparm, lpmfdata, m_nAux}); @@ -597,7 +597,7 @@ void PeleLM::fillpatch_phiV(int lev, auto* mapper = getInterpolator(); // Density - PhysBCFunct> crse_bndry_func(geom[lev-1], fetchBCRecArray(PHIV,1), + PhysBCFunct> crse_bndry_func(geom[lev-1], fetchBCRecArray(PHIV,1), PeleLMCCFillExtDirPhiV{lprobparm, lpmfdata, m_nAux}); PhysBCFunct> fine_bndry_func(geom[lev], fetchBCRecArray(PHIV,1), PeleLMCCFillExtDirPhiV{lprobparm, lpmfdata, m_nAux}); @@ -819,7 +819,7 @@ void PeleLM::setInflowBoundaryVel(MultiFab &a_vel, int lev, TimeStamp a_time) { BL_PROFILE_VAR("PeleLM::setInflowBoundaryVel()", setInflowBoundaryVel); - + Real time = getTime(lev, a_time); // Create a dummy BCRec from Velocity BCRec keeping only Inflow and set the other to bogus @@ -838,7 +838,7 @@ void PeleLM::setInflowBoundaryVel(MultiFab &a_vel, dummyVelBCRec[idim].setHi(idim2,BCType::bogus); } } - } + } fillTurbInflow(a_vel, 0, lev, time); diff --git a/Source/PeleLMData.cpp b/Source/PeleLMData.cpp index dd911385..34e6ed51 100644 --- a/Source/PeleLMData.cpp +++ b/Source/PeleLMData.cpp @@ -5,7 +5,7 @@ using namespace amrex; PeleLM::LevelData::LevelData(amrex::BoxArray const& ba, amrex::DistributionMapping const& dm, amrex::FabFactory const& factory, - int a_incompressible, int a_has_divu, + int a_incompressible, int a_has_divu, int a_nAux, int a_nGrowState) { if (a_incompressible ) { @@ -18,7 +18,7 @@ PeleLM::LevelData::LevelData(amrex::BoxArray const& ba, dm, 1 , 1 , MFInfo(), factory); visc_cc.define( ba, dm, 1 , 1 , MFInfo(), factory); if (! a_incompressible ) { - if (a_has_divu) { + if (a_has_divu) { divu.define (ba, dm, 1 , 1 , MFInfo(), factory); } diff_cc.define (ba, dm, NUM_SPECIES+2 , 1 , MFInfo(), factory); @@ -36,7 +36,7 @@ PeleLM::LevelData::LevelData(amrex::BoxArray const& ba, } PeleLM::LevelDataReact::LevelDataReact(const amrex::BoxArray &ba, - const amrex::DistributionMapping &dm, + const amrex::DistributionMapping &dm, const amrex::FabFactory &factory) { int IRsize = NUM_SPECIES; @@ -67,7 +67,7 @@ PeleLM::LevelDataNLSolve::LevelDataNLSolve(amrex::BoxArray const& ba, PeleLM::AdvanceDiffData::AdvanceDiffData(int a_finestLevel, const amrex::Vector &ba, - const amrex::Vector &dm, + const amrex::Vector &dm, const amrex::Vector>> &factory, int nGrowAdv, int a_use_wbar, @@ -109,7 +109,7 @@ PeleLM::AdvanceDiffData::AdvanceDiffData(int a_finestLevel, PeleLM::AdvanceAdvData::AdvanceAdvData(int a_finestLevel, const amrex::Vector &ba, - const amrex::Vector &dm, + const amrex::Vector &dm, const amrex::Vector>> &factory, int a_incompressible, int nGrowAdv, diff --git a/Source/PeleLMDerive.H b/Source/PeleLMDerive.H index 91d6398b..e3a77e75 100644 --- a/Source/PeleLMDerive.H +++ b/Source/PeleLMDerive.H @@ -12,19 +12,20 @@ class PeleLM; typedef void (*PeleLMDeriveFunc) (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); typedef void (*PeleLMMFVecDeriveFunc) (amrex::Vector derfab, int dcomp, int ncomp, amrex::Vector statefab, + amrex::Vector reactfab, amrex::Vector pressfab, const amrex::Geometry &geomdata, amrex::Real time, const amrex::Vector &bcrec); class PeleLMDeriveRec -{ +{ friend class PeleLMDeriveList; - + public: typedef amrex::Box (*DeriveBoxMap)(const amrex::Box&); diff --git a/Source/PeleLMDerive.cpp b/Source/PeleLMDerive.cpp index d3a2a99f..9499fa19 100644 --- a/Source/PeleLMDerive.cpp +++ b/Source/PeleLMDerive.cpp @@ -66,7 +66,7 @@ PeleLMDeriveRec::PeleLMDeriveRec (const std::string& a_name, bx_map(box_map) {} -PeleLMDeriveRec::~PeleLMDeriveRec () +PeleLMDeriveRec::~PeleLMDeriveRec () { func = nullptr; mapper = 0; @@ -187,7 +187,7 @@ PeleLMDeriveList::dlist () } bool -PeleLMDeriveList::canDerive (const std::string& name) const +PeleLMDeriveList::canDerive (const std::string& name) const { for (std::list::const_iterator li = lst.begin(), End = lst.end(); li != End; diff --git a/Source/PeleLMDeriveFunc.H b/Source/PeleLMDeriveFunc.H index 55d21270..e9476edd 100644 --- a/Source/PeleLMDeriveFunc.H +++ b/Source/PeleLMDeriveFunc.H @@ -12,57 +12,67 @@ class PeleLM; void pelelm_dermassfrac (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_dermolefrac (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_dertemp (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); +void pelelm_derheatrelease (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, + const amrex::Geometry& geomdata, + amrex::Real time, const amrex::Vector &bcrec, int level); + void pelelm_derkineticenergy (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); +void pelelm_derenstrophy (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, + const amrex::Geometry& geomdata, + amrex::Real time, const amrex::Vector &bcrec, int level); + void pelelm_deravgpress (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_dermgvort (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_dermixfrac (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_derprogvar (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_dervisc (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_derdiffc (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); void pelelm_derlambda (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab, + const amrex::FArrayBox& statefab, const amrex::FArrayBox& reactfab, const amrex::FArrayBox& pressfab, const amrex::Geometry& geomdata, amrex::Real time, const amrex::Vector &bcrec, int level); #endif diff --git a/Source/PeleLMDeriveFunc.cpp b/Source/PeleLMDeriveFunc.cpp index c829af4e..83042279 100644 --- a/Source/PeleLMDeriveFunc.cpp +++ b/Source/PeleLMDeriveFunc.cpp @@ -11,7 +11,7 @@ using namespace amrex; // Extract temp // void pelelm_dertemp (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -29,11 +29,44 @@ void pelelm_dertemp (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dco }); } +// +// Compute heat release +// +void pelelm_derheatrelease (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, + const FArrayBox& statefab, const FArrayBox& reactfab, const FArrayBox& /*pressfab*/, + const Geometry& /*geomdata*/, + Real /*time*/, const Vector& /*bcrec*/, int level) + +{ + AMREX_ASSERT(derfab.box().contains(bx)); + AMREX_ASSERT(statefab.box().contains(bx)); + AMREX_ASSERT(derfab.nComp() >= dcomp + ncomp); + AMREX_ASSERT(!a_pelelm->m_incompressible); + + FArrayBox EnthFab; + EnthFab.resize(bx,NUM_SPECIES); + Elixir Enthi = EnthFab.elixir(); + + auto const temp = statefab.const_array(TEMP); + auto const react = reactfab.const_array(0); + auto const& Hi = EnthFab.array(); + auto HRR = derfab.array(dcomp); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + getHGivenT( i, j, k, temp, Hi ); + HRR(i,j,k) = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + HRR(i,j,k) -= Hi(i,j,k,n) * react(i,j,k,n); + } + }); +} + // // Extract species mass fractions Y_n // void pelelm_dermassfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -58,7 +91,7 @@ void pelelm_dermassfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int // Extract species mole fractions X_n // void pelelm_dermolefrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) { @@ -91,7 +124,7 @@ void pelelm_dermolefrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int // Compute cell averaged pressure from nodes // void pelelm_deravgpress (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, - const FArrayBox& /*statefab*/, const FArrayBox& pressfab, + const FArrayBox& /*statefab*/, const FArrayBox& /*reactfab*/, const FArrayBox& pressfab, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -119,7 +152,7 @@ void pelelm_deravgpress (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int // Compute vorticity magnitude // void pelelm_dermgvort (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& geomdata, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -164,7 +197,7 @@ void pelelm_dermgvort (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int d // Compute the kinetic energy // void pelelm_derkineticenergy (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -197,11 +230,80 @@ void pelelm_derkineticenergy (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab } } +// +// Compute enstrophy +// +void pelelm_derenstrophy (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, + const Geometry& geomdata, + Real /*time*/, const Vector& /*bcrec*/, int /*level*/) + +{ + + AMREX_D_TERM(const amrex::Real idx = geomdata.InvCellSize(0);, + const amrex::Real idy = geomdata.InvCellSize(1);, + const amrex::Real idz = geomdata.InvCellSize(2);); + + // TODO : EB + // TODO : BCs + + if (a_pelelm->m_incompressible) { + auto const& dat_arr = statefab.const_array(VELX); + auto const& ens_arr = derfab.array(dcomp); + amrex::ParallelFor(bx, [=,rho=a_pelelm->m_rho] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { +#if ( AMREX_SPACEDIM == 2 ) + amrex::Real vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + amrex::Real uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; + ens_arr(i,j,k) = 0.5 * rho * (vx-uy)*(vx-uy); + +#elif ( AMREX_SPACEDIM == 3 ) + amrex::Real vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + amrex::Real wx = 0.5 * (dat_arr(i+1,j,k,2) - dat_arr(i-1,j,k,2)) * idx; + + amrex::Real uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; + amrex::Real wy = 0.5 * (dat_arr(i,j+1,k,2) - dat_arr(i,j-1,k,2)) * idy; + + amrex::Real uz = 0.5 * (dat_arr(i,j,k+1,0) - dat_arr(i,j,k-1,0)) * idz; + amrex::Real vz = 0.5 * (dat_arr(i,j,k+1,1) - dat_arr(i,j,k-1,1)) * idz; + + ens_arr(i,j,k) = 0.5 * rho * ((wy-vz)*(wy-vz) + (uz-wx)*(uz-wx) + (vx-uy)*(vx-uy)); +#endif + }); + } else { + auto const& dat_arr = statefab.const_array(VELX); + auto const& rho_arr = statefab.const_array(DENSITY); + auto const& ens_arr = derfab.array(dcomp); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { +#if ( AMREX_SPACEDIM == 2 ) + amrex::Real vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + amrex::Real uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; + ens_arr(i,j,k) = 0.5 * rho_arr(i,j,k) * (vx-uy)*(vx-uy); + +#elif ( AMREX_SPACEDIM == 3 ) + amrex::Real vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + amrex::Real wx = 0.5 * (dat_arr(i+1,j,k,2) - dat_arr(i-1,j,k,2)) * idx; + + amrex::Real uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; + amrex::Real wy = 0.5 * (dat_arr(i,j+1,k,2) - dat_arr(i,j-1,k,2)) * idy; + + amrex::Real uz = 0.5 * (dat_arr(i,j,k+1,0) - dat_arr(i,j,k-1,0)) * idz; + amrex::Real vz = 0.5 * (dat_arr(i,j,k+1,1) - dat_arr(i,j,k-1,1)) * idz; + + ens_arr(i,j,k) = 0.5 * rho_arr(i,j,k) * ((wy-vz)*(wy-vz) + (uz-wx)*(uz-wx) + (vx-uy)*(vx-uy)); +#endif + }); + } + +} + + // // Compute mixture fraction // void pelelm_dermixfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) @@ -240,14 +342,14 @@ void pelelm_dermixfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int // Compute progress variable // void pelelm_derprogvar (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) { AMREX_ASSERT(derfab.box().contains(bx)); AMREX_ASSERT(statefab.box().contains(bx)); - AMREX_ASSERT(ncomp == 1); + AMREX_ASSERT(ncomp == 1); if (a_pelelm->m_C0 < 0.0) amrex::Abort("Progress variable not initialized"); @@ -262,22 +364,22 @@ void pelelm_derprogvar (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int amrex::GpuArray Cweights; for (int n=0; nm_Cweights[n]; - } + } amrex::ParallelFor(bx, [=,revert=a_pelelm->m_Crevert] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { + { amrex::Real rho_inv = 1.0_rt / density(i,j,k); prog_var(i,j,k) = 0.0_rt; for (int n = 0; n& /*bcrec*/, int /*level*/) { @@ -312,7 +414,7 @@ void pelelm_dervisc (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dco // Extract mixture averaged species diffusion coefficients // void pelelm_derdiffc (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) { @@ -339,7 +441,7 @@ void pelelm_derdiffc (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dc // Extract thermal diffusivity // void pelelm_derlambda (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp, - const FArrayBox& statefab, const FArrayBox& /*pressfab*/, + const FArrayBox& statefab, const FArrayBox& /*reactfab*/, const FArrayBox& /*pressfab*/, const Geometry& /*geomdata*/, Real /*time*/, const Vector& /*bcrec*/, int /*level*/) { diff --git a/Source/PeleLMDiffusion.cpp b/Source/PeleLMDiffusion.cpp index 7d689d84..95c17eab 100644 --- a/Source/PeleLMDiffusion.cpp +++ b/Source/PeleLMDiffusion.cpp @@ -88,13 +88,13 @@ void PeleLM::computeDifferentialDiffusionTerms(const TimeStamp &a_time, GetVecOfPtrs(diffData->Dn), 0, GetVecOfArrOfPtrs(fluxes), 0, NUM_SPECIES, 1, bcRecSpec_d.dataPtr(), -1.0, m_dt); - auto bcRecTemp = fetchBCRecArray(TEMP,1); + auto bcRecTemp = fetchBCRecArray(TEMP,1); auto bcRecTemp_d = convertToDeviceVector(bcRecTemp); fluxDivergenceRD(GetVecOfConstPtrs(getTempVect(AmrOldTime)), 0, GetVecOfPtrs(diffData->Dn), NUM_SPECIES, GetVecOfArrOfPtrs(fluxes), NUM_SPECIES, 1, 1, bcRecTemp_d.dataPtr(), -1.0, m_dt); - auto bcRecRhoH = fetchBCRecArray(RHOH,1); + auto bcRecRhoH = fetchBCRecArray(RHOH,1); auto bcRecRhoH_d = convertToDeviceVector(bcRecRhoH); fluxDivergenceRD(GetVecOfConstPtrs(getRhoHVect(AmrOldTime)), 0, GetVecOfPtrs(diffData->Dn), NUM_SPECIES+1, @@ -111,13 +111,13 @@ void PeleLM::computeDifferentialDiffusionTerms(const TimeStamp &a_time, GetVecOfPtrs(diffData->Dnp1), 0, GetVecOfArrOfPtrs(fluxes), 0, NUM_SPECIES, 1, bcRecSpec_d.dataPtr(), -1.0, m_dt); - auto bcRecTemp = fetchBCRecArray(TEMP,1); + auto bcRecTemp = fetchBCRecArray(TEMP,1); auto bcRecTemp_d = convertToDeviceVector(bcRecTemp); fluxDivergenceRD(GetVecOfConstPtrs(getTempVect(AmrNewTime)), 0, GetVecOfPtrs(diffData->Dnp1), NUM_SPECIES, GetVecOfArrOfPtrs(fluxes), NUM_SPECIES, 1, 1, bcRecTemp_d.dataPtr(), -1.0, m_dt); - auto bcRecRhoH = fetchBCRecArray(RHOH,1); + auto bcRecRhoH = fetchBCRecArray(RHOH,1); auto bcRecRhoH_d = convertToDeviceVector(bcRecRhoH); fluxDivergenceRD(GetVecOfConstPtrs(getRhoHVect(AmrNewTime)), 0, GetVecOfPtrs(diffData->Dnp1), NUM_SPECIES+1, @@ -505,7 +505,7 @@ void PeleLM::computeSpeciesEnthalpyFlux(const Vector scratch = tmpfab.array(0); if (m_adv_redist_type == "FluxRedist") - { + { amrex::ParallelFor(Box(scratch), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { scratch(i,j,k) = 1.;}); // TODO might want to test volfrac @@ -153,7 +153,7 @@ void PeleLM::redistributeDiff(int a_lev, Elixir eli = tmpfab.elixir(); Array4 scratch = tmpfab.array(0); if (m_diff_redist_type == "FluxRedist") - { + { amrex::ParallelFor(Box(scratch), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { scratch(i,j,k) = 1.;}); // TODO might want to test volfrac @@ -311,7 +311,7 @@ void PeleLM::initialRedistribution() void PeleLM::getEBDistance(int a_lev, MultiFab &a_signDistLev) { - + if (a_lev == 0) { MultiFab::Copy(a_signDistLev,*m_signedDist0,0,0,1,0); @@ -329,8 +329,8 @@ void PeleLM::getEBDistance(int a_lev, // Get signDist on coarsen fineBA BoxArray coarsenBA(grids[lev].size()); - for (int j = 0, N = coarsenBA.size(); j < N; ++j) - { + for (int j = 0, N = coarsenBA.size(); j < N; ++j) + { coarsenBA.set(j,interpolater.CoarseBox(grids[lev][j], refRatio(lev-1))); } MultiFab coarsenSignDist(coarsenBA,dmap[lev],1,0); @@ -350,7 +350,7 @@ void PeleLM::getEBDistance(int a_lev, interpolater.interp(coarsenSignDist, 0, *currentSignDist, 0, 1, IntVect(0), - Geom(lev-1), Geom(lev), + Geom(lev-1), Geom(lev), Geom(lev).Domain(), refRatio(lev-1), {m_bcrec_force},0); diff --git a/Source/PeleLMEos.cpp b/Source/PeleLMEos.cpp index 3c52d4ae..af56b539 100644 --- a/Source/PeleLMEos.cpp +++ b/Source/PeleLMEos.cpp @@ -125,7 +125,7 @@ void PeleLM::calcDivU(int is_init, { divu(i,j,k) = 0.0; }); - } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes amrex::ParallelFor(bx, [ rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, use_react, flag] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -285,8 +285,8 @@ PeleLM::adjustPandDivU(std::unique_ptr &advData) amrex::ParallelFor(bx, [=,pOld=m_pOld,pNew=m_pNew] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real gammaInv_o = getGammaInv(i,j,k,rhoYo,T_o); - Real gammaInv_n = getGammaInv(i,j,k,rhoYn,T_n); + Real gammaInv_o = getGammaInv(i,j,k,rhoYo,T_o); + Real gammaInv_n = getGammaInv(i,j,k,rhoYn,T_n); theta(i,j,k) = 0.5 * (gammaInv_o/pOld + gammaInv_n/pNew); }); } @@ -312,11 +312,11 @@ PeleLM::adjustPandDivU(std::unique_ptr &advData) // mac_divu is now delta_S advData->mac_divu[lev].plus(-Sbar,0,1); } - + // Compute 1/Volume * int(U_inflow)dA across all boundary faces amrex::Real umacFluxBalance = AMREX_D_TERM( m_domainUmacFlux[0] + m_domainUmacFlux[1], - + m_domainUmacFlux[2] + m_domainUmacFlux[3], - + m_domainUmacFlux[4] + m_domainUmacFlux[5]); + + m_domainUmacFlux[2] + m_domainUmacFlux[3], + + m_domainUmacFlux[4] + m_domainUmacFlux[5]); Real divu_vol = umacFluxBalance/m_uncoveredVol; // Advance the ambient pressure diff --git a/Source/PeleLMEvaluate.cpp b/Source/PeleLMEvaluate.cpp index 5a01a9f7..a8623d2e 100644 --- a/Source/PeleLMEvaluate.cpp +++ b/Source/PeleLMEvaluate.cpp @@ -57,7 +57,7 @@ void PeleLM::Evaluate() { // Regular derived functions and State entries are called on a per level basis // derive function can handle both derived and state entries - } else if ( derive_lst.canDerive(m_evaluatePlotVars[ivar]) + } else if ( derive_lst.canDerive(m_evaluatePlotVars[ivar]) || isStateVariable(m_evaluatePlotVars[ivar]) ) { for (int lev = 0; lev <= finest_level; ++lev) { std::unique_ptr mf; @@ -67,7 +67,7 @@ void PeleLM::Evaluate() { } } cnt += cntIncr; - } + } //---------------------------------------------------------------- // Write the evaluated variables to disc @@ -119,12 +119,12 @@ PeleLM::MLevaluate(const Vector &a_MFVec, Abort("advTerm not available yet, soon hopefully"); } else if ( a_var == "instRR" ) { for (int lev = 0; lev <= finest_level; ++lev) { - std::unique_ptr I_RR = std::make_unique (*a_MFVec[lev],amrex::make_alias,a_comp,NUM_SPECIES); + std::unique_ptr I_RR = std::make_unique (*a_MFVec[lev],amrex::make_alias,a_comp,NUM_SPECIES); computeInstantaneousReactionRate(lev, AmrNewTime, I_RR.get()); } nComp = NUM_SPECIES; } else if ( a_var == "transportCC" ) { - // Cell-centered transport coefficients functions go through the level + // Cell-centered transport coefficients functions go through the level // data container. Simply copy once the later has been filled. calcViscosity(AmrNewTime); calcDiffusivity(AmrNewTime); diff --git a/Source/PeleLMEvolve.cpp b/Source/PeleLMEvolve.cpp index e52c4559..4d2a87c8 100644 --- a/Source/PeleLMEvolve.cpp +++ b/Source/PeleLMEvolve.cpp @@ -10,7 +10,7 @@ void PeleLM::Evolve() { int plt_justDidIt = 0; int chk_justDidIt = 0; - + while(!do_not_evolve) { plt_justDidIt = 0; @@ -93,10 +93,10 @@ PeleLM::writePlotNow() if ( m_plot_int > 0 && (m_nstep % m_plot_int == 0) ) { write_now = true; - + } else if ( m_plot_per_exact > 0.0 && (std::abs(std::remainder(m_cur_time, m_plot_per_exact)) < 1.e-12) ) { write_now = true; - + } else if (m_plot_per_approx > 0.0) { // Check to see if we've crossed a plot_per interval by comparing // the number of intervals that have elapsed for both the current diff --git a/Source/PeleLMForces.cpp b/Source/PeleLMForces.cpp index 52dc2cab..beebc79a 100644 --- a/Source/PeleLMForces.cpp +++ b/Source/PeleLMForces.cpp @@ -51,7 +51,7 @@ void PeleLM::getVelForces(const TimeStamp &a_time, #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(*a_velForce,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + { const auto& bx = mfi.tilebox(); FArrayBox DummyFab(bx,1); const auto& vel_arr = ldata_p->state.const_array(mfi,VELX); @@ -78,20 +78,20 @@ void PeleLM::getVelForces(const TimeStamp &a_time, if ( add_gradP || has_divTau ) { const auto& gp_arr = (add_gradP) ? ldataGP_p->gp.const_array(mfi) : DummyFab.array(); const auto& divTau_arr = (has_divTau) ? a_divTau->const_array(mfi) : DummyFab.array(); - amrex::ParallelFor(bx, [incomp_rho_inv, is_incomp, add_gradP, has_divTau, rho_arr, + amrex::ParallelFor(bx, [incomp_rho_inv, is_incomp, add_gradP, has_divTau, rho_arr, gp_arr, divTau_arr, force_arr] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { if ( is_incomp ) { for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - if (add_gradP) force_arr(i,j,k,idim) -= gp_arr(i,j,k,idim); - if (has_divTau) force_arr(i,j,k,idim) += divTau_arr(i,j,k,idim); + if (add_gradP) force_arr(i,j,k,idim) -= gp_arr(i,j,k,idim); + if (has_divTau) force_arr(i,j,k,idim) += divTau_arr(i,j,k,idim); force_arr(i,j,k,idim) *= incomp_rho_inv; } } else { for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { - if (add_gradP) force_arr(i,j,k,idim) -= gp_arr(i,j,k,idim); - if (has_divTau) force_arr(i,j,k,idim) += divTau_arr(i,j,k,idim); + if (add_gradP) force_arr(i,j,k,idim) -= gp_arr(i,j,k,idim); + if (has_divTau) force_arr(i,j,k,idim) += divTau_arr(i,j,k,idim); force_arr(i,j,k,idim) /= rho_arr(i,j,k); } } @@ -138,8 +138,8 @@ void PeleLM::getVelForces(int lev, amrex::ParallelFor(bx, [=,grav=m_gravity, gp0=m_background_gp] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - makeVelForce(i,j,k, is_incomp, rho_incomp, + makeVelForce(i,j,k, is_incomp, rho_incomp, pseudo_gravity, a_time, grav, gp0, dV_control, dx, vel, rho, rhoY, rhoh, temp, extMom, extRho, force); - }); + }); } diff --git a/Source/PeleLMInit.cpp b/Source/PeleLMInit.cpp index 5c8f7830..9dd05474 100644 --- a/Source/PeleLMInit.cpp +++ b/Source/PeleLMInit.cpp @@ -104,7 +104,7 @@ void PeleLM::MakeNewLevelFromScratch( int lev, if ( lev == 0 && m_signDistNeeded) { // Set up CC signed distance container to control EB refinement m_signedDist0.reset(new MultiFab(grids[lev], dmap[lev], 1, 1, MFInfo(), *m_factory[lev])); - + // Estimate the maximum distance we need in terms of level 0 dx: Real extentFactor = static_cast(nErrorBuf(0)); for (int ilev = 1; ilev <= max_level; ++ilev) { @@ -384,7 +384,7 @@ void PeleLM::initLevelData(int lev) { { pelelm_initdata(i, j, k, m_incompressible, state_arr, aux_arr, #ifdef PELE_USE_EFIELD - ne_arr, phiV_arr, + ne_arr, phiV_arr, #endif geomdata, *lprobparm, lpmfdata); }); diff --git a/Source/PeleLMPlot.cpp b/Source/PeleLMPlot.cpp index 2c29df8b..cb2f7109 100644 --- a/Source/PeleLMPlot.cpp +++ b/Source/PeleLMPlot.cpp @@ -353,36 +353,36 @@ void PeleLM::WriteHeader(const std::string& name, bool is_checkpoint) const void PeleLM::WriteCheckPointFile() { BL_PROFILE("PeleLM::WriteCheckPointFile()"); - + const std::string& checkpointname = amrex::Concatenate(m_check_file, m_nstep, m_ioDigits); - + if (m_verbose) { amrex::Print() << "\n Writting checkpoint file: " << checkpointname << "\n"; } - + amrex::PreBuildDirectorHierarchy(checkpointname, level_prefix, finest_level + 1, true); bool is_checkpoint = true; WriteHeader(checkpointname, is_checkpoint); WriteJobInfo(checkpointname); - + for(int lev = 0; lev <= finest_level; ++lev) - { + { VisMF::Write(m_leveldata_new[lev]->state, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "state")); - + VisMF::Write(m_leveldata_new[lev]->gp, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "gradp")); - + VisMF::Write(m_leveldata_new[lev]->press, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "p")); - + if (!m_incompressible) { if (m_has_divu) { VisMF::Write(m_leveldata_new[lev]->divu, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "divU")); } - + if (m_do_react) { VisMF::Write(m_leveldatareact[lev]->I_R, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "I_R")); @@ -395,7 +395,7 @@ void PeleLM::WriteCheckPointFile() VisMF::Write(m_leveldata_new[lev]->nE, amrex::MultiFabFileFullPrefix(lev, checkpointname, level_prefix, "nE")); #endif - } + } } #ifdef PELELM_USE_SPRAY if (theSprayPC() != nullptr && do_spray_particles) { @@ -431,12 +431,12 @@ void PeleLM::ReadCheckPointFile() Vector fileCharPtr; ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr); std::string fileCharPtrString(fileCharPtr.dataPtr()); - std::istringstream is(fileCharPtrString, std::istringstream::in); + std::istringstream is(fileCharPtrString, std::istringstream::in); std::string line, word; - // Start reading from checkpoint file - + // Start reading from checkpoint file + // Title line std::getline(is, line); @@ -590,8 +590,8 @@ void PeleLM::initLevelDataFromPlt(int a_lev, int idT = -1, idV = -1, idY = -1, nSpecPlt = 0; for (int i = 0; i < plt_vars.size(); ++i) { std::string firstChars = plt_vars[i].substr(0, 2); - if (plt_vars[i] == "temp") idT = i; - if (plt_vars[i] == "x_velocity") idV = i; + if (plt_vars[i] == "temp") idT = i; + if (plt_vars[i] == "x_velocity") idV = i; if (firstChars == "Y(" && idY < 0 ) { // species might not be ordered in the order of the current mech. idY = i; } @@ -662,7 +662,7 @@ void PeleLM::initLevelDataFromPlt(int a_lev, sumYs += massfrac[n]; } } - massfrac[N2_ID] = 1.0 - sumYs; + massfrac[N2_ID] = 1.0 - sumYs; // Get density Real P_cgs = lprobparm->P_mean * 10.0; @@ -674,7 +674,7 @@ void PeleLM::initLevelDataFromPlt(int a_lev, Real h_cgs = 0.0; eos.TY2H(temp_arr(i,j,k), massfrac, h_cgs); rhoH_arr(i,j,k) = h_cgs * 1.0e-4 * rho_arr(i,j,k); - + // Fill rhoYs for (int n = 0; n < NUM_SPECIES; n++){ rhoY_arr(i,j,k,n) = massfrac[n] * rho_arr(i,j,k); @@ -766,7 +766,7 @@ void PeleLM::WriteJobInfo(const std::string& path) const jobInfoFile << " maximum zones = "; for(int idim = 0; idim < AMREX_SPACEDIM; idim++) { - jobInfoFile << geom[lev].Domain().length(idim) << " "; + jobInfoFile << geom[lev].Domain().length(idim) << " "; } jobInfoFile << "\n\n"; } diff --git a/Source/PeleLMReactions.cpp b/Source/PeleLMReactions.cpp index ca533318..7034e2ee 100644 --- a/Source/PeleLMReactions.cpp +++ b/Source/PeleLMReactions.cpp @@ -177,7 +177,7 @@ void PeleLM::advanceChemistry(int lev, // This advanceChemistry is called on all but the finest level // It works with BoxArrays built such that each box is either covered -// or uncovered and chem. integrator is called only on uncovered boxes +// or uncovered and chem. integrator is called only on uncovered boxes // the averaged down version of I_R is linearly added to the AD forcing // to build the t^{np1} solution on covered boxes. void PeleLM::advanceChemistry(int lev, @@ -278,7 +278,7 @@ void PeleLM::advanceChemistry(int lev, dt_incr, time_chem #ifdef AMREX_USE_GPU , amrex::Gpu::gpuStream() -#endif +#endif ); } else { // Use forcing and averaged down IR to advance species/rhoH/temp @@ -414,7 +414,7 @@ void PeleLM::computeInstantaneousReactionRate(int lev, { rhoYdot(i,j,k,n) = 0.0; }); - } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes amrex::ParallelFor(bx, [rhoY, rhoH, T, rhoYdot, flag] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -497,13 +497,13 @@ void PeleLM::getHeatRelease(int a_lev, auto const& Hi = EnthFab.array(); auto const& HRR = a_HR->array(mfi); amrex::ParallelFor(bx, [T, Hi, HRR, react] - AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { + AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { getHGivenT( i, j, k, T, Hi ); - HRR(i,j,k) = 0.0; + HRR(i,j,k) = 0.0; for (int n = 0; n < NUM_SPECIES; n++) { HRR(i,j,k) -= Hi(i,j,k,n) * react(i,j,k,n); - } + } }); } } diff --git a/Source/PeleLMSetup.cpp b/Source/PeleLMSetup.cpp index a2e0156d..2c735557 100644 --- a/Source/PeleLMSetup.cpp +++ b/Source/PeleLMSetup.cpp @@ -134,7 +134,7 @@ void PeleLM::Setup() { m_pNew = prob_parm->P_mean; // Copy problem parameters into device copy - Gpu::copy(Gpu::hostToDevice, prob_parm, prob_parm+1,prob_parm_d); + Gpu::copy(Gpu::hostToDevice, prob_parm, prob_parm+1,prob_parm_d); } void PeleLM::readParameters() { @@ -227,7 +227,7 @@ void PeleLM::readParameters() { // Get the phiV bc ppef.getarr("phiV_lo_bc",lo_bc_char,0,AMREX_SPACEDIM); ppef.getarr("phiV_hi_bc",hi_bc_char,0,AMREX_SPACEDIM); - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { if (lo_bc_char[idim] == "Interior"){ m_phiV_bc.setLo(idim,0); @@ -237,7 +237,7 @@ void PeleLM::readParameters() { m_phiV_bc.setLo(idim,2); } else { amrex::Abort("Wrong PhiV bc. Should be : Interior, Dirichlet or Neumann"); - } + } if (hi_bc_char[idim] == "Interior"){ m_phiV_bc.setHi(idim,0); } else if (hi_bc_char[idim] == "Dirichlet") { @@ -298,12 +298,12 @@ void PeleLM::readParameters() { pp.queryarr("gravity", grav, 0, AMREX_SPACEDIM); Vector gp0(AMREX_SPACEDIM,0); pp.queryarr("gradP0", gp0, 0, AMREX_SPACEDIM); - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { m_background_gp[idim] = gp0[idim]; m_gravity[idim] = grav[idim]; } - + // ----------------------------------------- // diffusion @@ -312,7 +312,7 @@ void PeleLM::readParameters() { pp.query("deltaT_iterMax",m_deltaTIterMax); pp.query("deltaT_tol",m_deltaT_norm_max); pp.query("deltaT_crashIfFailing",m_crashOnDeltaTFail); - + // ----------------------------------------- // initialization @@ -359,7 +359,7 @@ void PeleLM::readParameters() { m_Godunov_ppm = 0; } else { Abort("Unknown 'advection_scheme'. Recognized options are: Godunov_PLM, Godunov_PPM or Godunov_BDS"); - } + } m_predict_advection_type = "Godunov"; // Only option at this point. This will disapear when predict_velocity support BDS. // ----------------------------------------- @@ -381,8 +381,11 @@ void PeleLM::readParameters() { // Temporals // ----------------------------------------- pp.query("do_temporals",m_do_temporals); - pp.query("temporal_int",m_temp_int); - pp.query("mass_balance",m_do_massBalance); + if (m_do_temporals) { + pp.query("temporal_int",m_temp_int); + pp.query("do_extremas",m_do_extremas); + pp.query("do_mass_balance",m_do_massBalance); + } // ----------------------------------------- // Time stepping control @@ -410,7 +413,7 @@ void PeleLM::readParameters() { m_EB_refine_type != "Adaptive" ) { Abort("refine_EB_type can only be 'Static' or 'Adaptive'"); } - // Default EB refinement level is max_level + // Default EB refinement level is max_level m_EB_refine_LevMax = max_level; pp.query("refine_EB_max_level",m_EB_refine_LevMax); pp.query("refine_EB_buffer",m_derefineEBBuffer); @@ -508,7 +511,7 @@ void PeleLM::variablesSetup() { std::string PrettyLine = std::string(78, '=') + "\n"; //---------------------------------------------------------------- - // Variables ordering is defined through compiler macro in PeleLM_Index.H + // Variables ordering is defined through compiler macro in PeleLM_Index.H // Simply print on screen the state layout and append to the stateComponents list Print() << "\n"; Print() << PrettyLine; @@ -722,6 +725,9 @@ void PeleLM::derivedSetup() derive_lst.add("diffcoeff",IndexType::TheCellType(),NUM_SPECIES, var_names_massfrac,pelelm_derdiffc,the_same_box); + // Heat Release + derive_lst.add("HeatRelease",IndexType::TheCellType(),1,pelelm_derheatrelease,the_same_box); + // Thermal diffusivity derive_lst.add("lambda",IndexType::TheCellType(),1,pelelm_derlambda,the_same_box); @@ -745,6 +751,9 @@ void PeleLM::derivedSetup() // Kinetic energy derive_lst.add("kinetic_energy",IndexType::TheCellType(),1,pelelm_derkineticenergy,the_same_box); + // Enstrophy + derive_lst.add("enstrophy",IndexType::TheCellType(),1,pelelm_derenstrophy,grow_box_by_two); + #ifdef PELE_USE_EFIELD // Charge distribution derive_lst.add("chargedistrib",IndexType::TheCellType(),1,pelelm_derchargedist,the_same_box); diff --git a/Source/PeleLMTagging.cpp b/Source/PeleLMTagging.cpp index bb6b310f..7bb152d0 100644 --- a/Source/PeleLMTagging.cpp +++ b/Source/PeleLMTagging.cpp @@ -39,7 +39,7 @@ PeleLM::ErrorEst( int lev, MultiFab signDist(grids[lev],dmap[lev],1,0,MFInfo(),EBFactory(lev)); getEBDistance(lev, signDist); //VisMF::Write(signDist,"signDistLev"+std::to_string(lev)); - + // Estimate how far I need to derefine Real diagFac = std::sqrt(2.0) * m_derefineEBBuffer; Real clearTagDist = Geom(m_EB_refine_LevMax).CellSize(0) * static_cast(nErrorBuf(m_EB_refine_LevMax)) * diagFac; @@ -53,9 +53,9 @@ PeleLM::ErrorEst( int lev, #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(tags,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + { const auto& bx = mfi.tilebox(); - const auto& dist = signDist.const_array(mfi); + const auto& dist = signDist.const_array(mfi); auto tag = tags.array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) diff --git a/Source/PeleLMTemporals.cpp b/Source/PeleLMTemporals.cpp index 766efba7..38cf3254 100644 --- a/Source/PeleLMTemporals.cpp +++ b/Source/PeleLMTemporals.cpp @@ -41,8 +41,8 @@ void PeleLM::massBalance() m_massNew = MFSum(GetVecOfConstPtrs(getDensityVect(AmrNewTime)),0); Real dmdt = (m_massNew - m_massOld) / m_dt; Real massFluxBalance = AMREX_D_TERM( m_domainMassFlux[0] + m_domainMassFlux[1], - + m_domainMassFlux[2] + m_domainMassFlux[3], - + m_domainMassFlux[4] + m_domainMassFlux[5]); + + m_domainMassFlux[2] + m_domainMassFlux[3], + + m_domainMassFlux[4] + m_domainMassFlux[5]); tmpMassFile << m_nstep << " " << m_cur_time // Time info << " " << m_massNew // mass @@ -80,7 +80,7 @@ void PeleLM::addMassFluxes(const Array &a_fluxes Real sumLo = 0.0; Real sumHi = 0.0; - sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -96,7 +96,7 @@ void PeleLM::addMassFluxes(const Array &a_fluxes return t; }); - sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -149,7 +149,7 @@ void PeleLM::addUmacFluxes(std::unique_ptr &advData, const Geome Real sumLo = 0.0; Real sumHi = 0.0; - sumLo = amrex::ReduceSum(advData->umac[lev][idim], 0, [=] + sumLo = amrex::ReduceSum(advData->umac[lev][idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -163,7 +163,7 @@ void PeleLM::addUmacFluxes(std::unique_ptr &advData, const Geome return t; }); - sumHi = amrex::ReduceSum(advData->umac[lev][idim], 0, [=] + sumHi = amrex::ReduceSum(advData->umac[lev][idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -190,8 +190,8 @@ void PeleLM::rhoHBalance() m_RhoHNew = MFSum(GetVecOfConstPtrs(getRhoHVect(AmrNewTime)),0); Real dRhoHdt = (m_RhoHNew - m_RhoHOld) / m_dt; Real rhoHFluxBalance = AMREX_D_TERM( m_domainRhoHFlux[0] + m_domainRhoHFlux[1], - + m_domainRhoHFlux[2] + m_domainRhoHFlux[3], - + m_domainRhoHFlux[4] + m_domainRhoHFlux[5]); + + m_domainRhoHFlux[2] + m_domainRhoHFlux[3], + + m_domainRhoHFlux[4] + m_domainRhoHFlux[5]); tmpMassFile << m_nstep << " " << m_cur_time // Time info << " " << m_RhoHNew // RhoH @@ -229,7 +229,7 @@ void PeleLM::addRhoHFluxes(const Array &a_fluxes Real sumLo = 0.0; Real sumHi = 0.0; - sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -244,7 +244,7 @@ void PeleLM::addRhoHFluxes(const Array &a_fluxes return t; }); - sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -299,7 +299,7 @@ void PeleLM::addRhoYFluxes(const Array &a_fluxes Real sumLo = 0.0; Real sumHi = 0.0; - sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumLo = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -314,7 +314,7 @@ void PeleLM::addRhoYFluxes(const Array &a_fluxes return t; }); - sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] + sumHi = amrex::ReduceSum(*a_fluxes[idim], 0, [=] AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& flux) -> Real { Real t = 0.0; @@ -348,33 +348,49 @@ void PeleLM::writeTemporals() //---------------------------------------------------------------- // State - // Get kinetic energy + // Get kinetic energy and enstrophy Vector> kinEnergy(finest_level + 1); + Vector> enstrophy(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { kinEnergy[lev] = derive("kinetic_energy", m_cur_time, lev, 0); + enstrophy[lev] = derive("enstrophy", m_cur_time, lev, 0); } - Real kinetic_energy = MFSum(GetVecOfConstPtrs(kinEnergy),0); + Real kinenergy_int = MFSum(GetVecOfConstPtrs(kinEnergy),0); + Real enstrophy_int = MFSum(GetVecOfConstPtrs(enstrophy),0); // Combustion Real fuelConsumptionInt = 0.0; Real heatReleaseRateInt = 0.0; if (fuelID > 0 && !(m_chem_integrator == "ReactorNull")) { fuelConsumptionInt = MFSum(GetVecOfConstPtrs(getIRVect()),fuelID); - for (int lev = 0; lev <= finest_level; ++lev) { + for (int lev = 0; lev <= finest_level; ++lev) { getHeatRelease(lev, kinEnergy[lev].get()); // Re-use kinEnergy container } heatReleaseRateInt = MFSum(GetVecOfConstPtrs(kinEnergy),0); } - // Get min/max/mean for non-species state components - - tmpStateFile << m_nstep << " " << m_cur_time // Time - << " " << kinetic_energy // Kinetic energy + tmpStateFile << m_nstep << " " << m_cur_time << " " << m_dt // Time + << " " << kinenergy_int // Kinetic energy + << " " << enstrophy_int // Enstrophy << " " << m_pNew // Thermo. pressure << " " << fuelConsumptionInt // Integ fuel burning rate << " " << heatReleaseRateInt // Integ heat release rate << " \n"; - tmpStateFile.flush(); + tmpStateFile.flush(); + + // Get min/max for state components + auto stateMax = ( m_incompressible ) ? MLmax(GetVecOfConstPtrs(getStateVect(AmrNewTime)),0,AMREX_SPACEDIM) + : MLmax(GetVecOfConstPtrs(getStateVect(AmrNewTime)),0,NVAR); + auto stateMin = ( m_incompressible ) ? MLmin(GetVecOfConstPtrs(getStateVect(AmrNewTime)),0,AMREX_SPACEDIM) + : MLmin(GetVecOfConstPtrs(getStateVect(AmrNewTime)),0,NVAR); + + tmpExtremasFile << m_nstep << " " << m_cur_time; // Time + for (int n = 0; n < stateMax.size(); ++n) { // Min & max of each state variable + tmpExtremasFile << " " << stateMin[n] << " " << stateMax[n]; + } + tmpExtremasFile << " \n"; + tmpExtremasFile.flush(); + } void PeleLM::openTempFile() @@ -393,6 +409,11 @@ void PeleLM::openTempFile() tmpMassFile.open(tempFileName.c_str(),std::ios::out | std::ios::app | std::ios_base::binary); tmpMassFile.precision(12); } + if (m_do_extremas) { + std::string tempFileName = "temporals/tempExtremas"; + tmpExtremasFile.open(tempFileName.c_str(),std::ios::out | std::ios::app | std::ios_base::binary); + tmpExtremasFile.precision(12); + } } } @@ -407,5 +428,9 @@ void PeleLM::closeTempFile() tmpMassFile.flush(); tmpMassFile.close(); } + if (m_do_extremas) { + tmpExtremasFile.flush(); + tmpExtremasFile.close(); + } } } diff --git a/Source/PeleLMTimestep.cpp b/Source/PeleLMTimestep.cpp index 5ee16293..a90cd02b 100644 --- a/Source/PeleLMTimestep.cpp +++ b/Source/PeleLMTimestep.cpp @@ -144,7 +144,7 @@ PeleLM::estDivUDt(const TimeStamp &a_time) { // Note: only method 1 of PeleLM is available here AMREX_ASSERT(m_divu_checkFlag==1); - + for (int lev = 0; lev <= finest_level; ++lev) { auto ldata_p = getLevelDataPtr(lev, a_time); diff --git a/Source/PeleLMTransportProp.cpp b/Source/PeleLMTransportProp.cpp index 79847668..a91b8641 100644 --- a/Source/PeleLMTransportProp.cpp +++ b/Source/PeleLMTransportProp.cpp @@ -122,7 +122,7 @@ PeleLM::getDiffusivity(int lev, int beta_comp, int ncomp, int doZeroVisc, for (MFIter mfi(beta_cc,TilingIfNotGPU()); mfi.isValid();++mfi) { for (int idim = 0; idim < AMREX_SPACEDIM; idim++) - { + { const Box ebx = mfi.nodaltilebox(idim); const Box& edomain = amrex::surroundingNodes(domain,idim); const auto& diff_c = beta_cc.const_array(mfi,beta_comp); @@ -140,11 +140,11 @@ PeleLM::getDiffusivity(int lev, int beta_comp, int ncomp, int doZeroVisc, } } #endif - + // Enable zeroing diffusivity on faces to produce walls if (doZeroVisc) { const auto geomdata = geom[lev].data(); - for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { + for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { const Box& edomain = amrex::surroundingNodes(domain,idim); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Source/PeleLMUMac.cpp b/Source/PeleLMUMac.cpp index 90ac47d4..0714fd24 100644 --- a/Source/PeleLMUMac.cpp +++ b/Source/PeleLMUMac.cpp @@ -10,7 +10,7 @@ void PeleLM::predictVelocity(std::unique_ptr &advData) { BL_PROFILE("PeleLM::predictVelocity()"); - // set umac boundaries to zero + // set umac boundaries to zero if ( advData->umac[0][0].nGrow() > 0 ) { for (int lev=0; lev <= finest_level; ++lev) { @@ -40,7 +40,7 @@ void PeleLM::predictVelocity(std::unique_ptr &advData) //---------------------------------------------------------------- // Predict face velocities at t^{n+1/2} with Godunov - auto bcRecVel = fetchBCRecArray(VELX,AMREX_SPACEDIM); + auto bcRecVel = fetchBCRecArray(VELX,AMREX_SPACEDIM); auto bcRecVel_d = convertToDeviceVector(bcRecVel); for (int lev = 0; lev <= finest_level; ++lev) { @@ -216,7 +216,7 @@ PeleLM::create_constrained_umac_grown(int a_lev, int a_nGrow, // Divergence preserving interp Interpolater* mapper = &face_divfree_interp; - + // Set BCRec for Umac Vector bcrec(1); for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { @@ -226,16 +226,16 @@ PeleLM::create_constrained_umac_grown(int a_lev, int a_nGrow, } else { bcrec[0].setLo(idim,BCType::foextrap); bcrec[0].setHi(idim,BCType::foextrap); - } - } + } + } Array,AMREX_SPACEDIM> bcrecArr = {AMREX_D_DECL(bcrec,bcrec,bcrec)}; - + PhysBCFunct> crse_bndry_func(*crse_geom, bcrec, umacFill{}); Array>,AMREX_SPACEDIM> cbndyFuncArr = {AMREX_D_DECL(crse_bndry_func,crse_bndry_func,crse_bndry_func)}; - + PhysBCFunct> fine_bndry_func(*fine_geom, bcrec, umacFill{}); Array>,AMREX_SPACEDIM> fbndyFuncArr = {AMREX_D_DECL(fine_bndry_func,fine_bndry_func,fine_bndry_func)}; - + // Use piecewise constant interpolation in time, so create dummy variable for time Real dummy = 0.; FillPatchTwoLevels(u_mac_fine, IntVect(a_nGrow), dummy, diff --git a/Source/PeleLMUtils.H b/Source/PeleLMUtils.H index a0246a28..1ccd1a2b 100644 --- a/Source/PeleLMUtils.H +++ b/Source/PeleLMUtils.H @@ -3,7 +3,7 @@ template amrex::Gpu::DeviceVector convertToDeviceVector ( amrex::Vector v) -{ +{ int ncomp = v.size(); amrex::Gpu::DeviceVector v_d(ncomp); #ifdef AMREX_USE_GPU diff --git a/Source/PeleLMUtils.cpp b/Source/PeleLMUtils.cpp index 5136c42d..d91d9b23 100644 --- a/Source/PeleLMUtils.cpp +++ b/Source/PeleLMUtils.cpp @@ -51,10 +51,10 @@ void PeleLM::fluxDivergenceRD(const Vector &a_state, } else { // Fluxes are extensive extFluxDivergenceLevel(lev, divTmp, 0, a_fluxes[lev], flux_comp, ncomp, scale); } - + // Need FillBoundary before redistribution divTmp.FillBoundary(geom[lev].periodicity()); - + // Redistribute diffusion term redistributeDiff(lev, a_dt, divTmp, 0, @@ -115,7 +115,7 @@ void PeleLM::extFluxDivergenceLevel(int lev, { divergence(i,j,k,n) = 0.0; }); - } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes auto vfrac = ebfact.getVolFrac().const_array(mfi); amrex::ParallelFor(bx, [ncomp, flag, vfrac, divergence, AMREX_D_DECL(fluxX, fluxY, fluxZ), vol, scale] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -213,7 +213,7 @@ void PeleLM::intFluxDivergenceLevel(int lev, { divergence(i,j,k,n) = 0.0; }); - } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes + } else if (flagfab.getType(bx) != FabType::regular ) { // EB containing boxes auto vfrac = ebfact.getVolFrac().const_array(mfi); AMREX_D_TERM( const auto& afrac_x = areafrac[0]->array(mfi);, const auto& afrac_y = areafrac[1]->array(mfi);, @@ -313,14 +313,14 @@ void PeleLM::advFluxDivergence(int a_lev, #ifdef AMREX_USE_EB auto const& ebfact = EBFactory(a_lev); #else - amrex::ignore_unused(a_lev); + amrex::ignore_unused(a_lev); #endif #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(a_divergence,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - + Box const& bx = mfi.tilebox(); // Get the divergence @@ -536,6 +536,7 @@ PeleLM::derive(const std::string &a_name, std::unique_ptr statemf = fillPatchState(lev, a_time, m_nGrowState); // Get pressure: TODO no fillpatch for pressure just yet, simply get new state auto ldata_p = getLevelDataPtr(lev,AmrNewTime); + auto ldataR_p = getLevelDataReactPtr(lev); auto stateBCs = fetchBCRecArray(VELX,NVAR); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -545,8 +546,9 @@ PeleLM::derive(const std::string &a_name, const Box& bx = mfi.growntilebox(nGrow); FArrayBox& derfab = (*mf)[mfi]; FArrayBox const& statefab = (*statemf)[mfi]; + FArrayBox const& reactfab = ldataR_p->I_R[mfi]; FArrayBox const& pressfab = ldata_p->press[mfi]; - rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, pressfab, geom[lev], a_time, stateBCs, lev); + rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, reactfab, pressfab, geom[lev], a_time, stateBCs, lev); } } else { // This is a state variable mf.reset(new MultiFab(grids[lev], dmap[lev], 1, nGrow)); @@ -582,6 +584,7 @@ PeleLM::deriveComp(const std::string &a_name, std::unique_ptr statemf = fillPatchState(lev, a_time, m_nGrowState); // Get pressure: TODO no fillpatch for pressure just yet, simply get new state auto ldata_p = getLevelDataPtr(lev,AmrNewTime); + auto ldataR_p = getLevelDataReactPtr(lev); auto stateBCs = fetchBCRecArray(VELX,NVAR); // Temp MF for all the derive components @@ -594,8 +597,9 @@ PeleLM::deriveComp(const std::string &a_name, const Box& bx = mfi.growntilebox(nGrow); FArrayBox& derfab = derTemp[mfi]; FArrayBox const& statefab = (*statemf)[mfi]; + FArrayBox const& reactfab = ldataR_p->I_R[mfi]; FArrayBox const& pressfab = ldata_p->press[mfi]; - rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, pressfab, geom[lev], a_time, stateBCs, lev); + rec->derFunc()(this, bx, derfab, 0, rec->numDerive(), statefab, reactfab, pressfab, geom[lev], a_time, stateBCs, lev); } // Copy into outgoing unique_ptr int derComp = rec->variableComp(a_name); @@ -634,19 +638,19 @@ PeleLM::initProgressVariable() int entryCount = pp.countval("progressVariable.weights"); stringIn.resize(entryCount); pp.getarr("progressVariable.weights",stringIn,0,entryCount); - parseVars(varNames, stringIn, weightsIn); + parseVars(varNames, stringIn, weightsIn); // Cold side/Hot side Vector coldState(NUM_SPECIES+1,0.0); entryCount = pp.countval("progressVariable.coldState"); stringIn.resize(entryCount); pp.getarr("progressVariable.coldState",stringIn,0,entryCount); - parseVars(varNames, stringIn, coldState); + parseVars(varNames, stringIn, coldState); Vector hotState(NUM_SPECIES+1,0.0); entryCount = pp.countval("progressVariable.hotState"); stringIn.resize(entryCount); pp.getarr("progressVariable.hotState",stringIn,0,entryCount); - parseVars(varNames, stringIn, hotState); + parseVars(varNames, stringIn, hotState); m_C0 = 0.0; m_C1 = 0.0; for (int i = 0; i < NUM_SPECIES+1; ++i) { @@ -804,7 +808,7 @@ PeleLM::fetchDiffTypeArray(int scomp, int ncomp) return types; } -Real +Real PeleLM::MFSum (const Vector &a_mf, int comp) { // Get the integral of the MF, not including the fine-covered and @@ -822,11 +826,11 @@ PeleLM::MFSum (const Vector &a_mf, int comp) // Use amrex::ReduceSum auto const& ebfact = dynamic_cast(Factory(lev)); auto const& vfrac = ebfact.getVolFrac(); - + Real sm = 0.0; if ( lev != finest_level ) { sm = amrex::ReduceSum(*a_mf[lev], vfrac, *m_coveredMask[lev], 0, [vol, comp] - AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& mf_arr, + AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& mf_arr, Array4 const& vf_arr, Array4 const& covered_arr) -> Real { @@ -839,7 +843,7 @@ PeleLM::MFSum (const Vector &a_mf, int comp) }); } else { sm = amrex::ReduceSum(*a_mf[lev], vfrac, 0, [vol, comp] - AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& mf_arr, + AMREX_GPU_HOST_DEVICE (Box const& bx, Array4 const& mf_arr, Array4 const& vf_arr) -> Real { Real sum = 0.0; @@ -912,7 +916,7 @@ void PeleLM::setTypicalValues(const TimeStamp &a_time, int is_init) for (int idim = 0; idim < AMREX_SPACEDIM; idim++) { typical_values[idim] = std::max(stateMax[VELX+idim],std::abs(stateMin[VELX+idim])); } - + if (!m_incompressible) { // Average between max/min typical_values[DENSITY] = 0.5 * (stateMax[DENSITY] + stateMin[DENSITY]); @@ -981,7 +985,7 @@ PeleLM::MFmax(const MultiFab *a_MF, #ifdef AMREX_USE_EB if ( a_MF->hasEBFabFactory() ) - { + { const auto& ebfactory = dynamic_cast(a_MF->Factory()); auto const& flags = ebfactory.getMultiEBCellFlagFab(); #ifdef AMREX_USE_GPU @@ -991,13 +995,13 @@ PeleLM::MFmax(const MultiFab *a_MF, auto const& mask = a_mask.const_arrays(); mx = ParReduce(TypeList{}, TypeList{}, *a_MF, IntVect(0), [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple - { + { if (flagsma[box_no](i,j,k).isCovered() || !mask[box_no](i,j,k)) { return AMREX_REAL_LOWEST; } else { return ma[box_no](i,j,k,comp); - } + } }); } else #endif @@ -1035,7 +1039,7 @@ PeleLM::MFmax(const MultiFab *a_MF, return AMREX_REAL_LOWEST; } else { return ma[box_no](i,j,k,comp); - } + } }); } else #endif @@ -1070,7 +1074,7 @@ PeleLM::MFmin(const MultiFab *a_MF, #ifdef AMREX_USE_EB if ( a_MF->hasEBFabFactory() ) - { + { const auto& ebfactory = dynamic_cast(a_MF->Factory()); auto const& flags = ebfactory.getMultiEBCellFlagFab(); #ifdef AMREX_USE_GPU @@ -1080,13 +1084,13 @@ PeleLM::MFmin(const MultiFab *a_MF, auto const& mask = a_mask.const_arrays(); mn = ParReduce(TypeList{}, TypeList{}, *a_MF, IntVect(0), [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple - { + { if (flagsma[box_no](i,j,k).isCovered() || !mask[box_no](i,j,k)) { return AMREX_REAL_MAX; } else { return ma[box_no](i,j,k,comp); - } + } }); } else #endif @@ -1124,7 +1128,7 @@ PeleLM::MFmin(const MultiFab *a_MF, return AMREX_REAL_MAX; } else { return ma[box_no](i,j,k,comp); - } + } }); } else #endif @@ -1246,11 +1250,11 @@ PeleLM::initMixtureFraction() int entryCount = pp.countval("mixtureFraction.oxidTank"); compositionIn.resize(entryCount); pp.getarr("mixtureFraction.oxidTank",compositionIn,0,entryCount); - parseComposition(compositionIn, MFCompoType, YO); + parseComposition(compositionIn, MFCompoType, YO); entryCount = pp.countval("mixtureFraction.fuelTank"); compositionIn.resize(entryCount); pp.getarr("mixtureFraction.fuelTank",compositionIn,0,entryCount); - parseComposition(compositionIn, MFCompoType, YF); + parseComposition(compositionIn, MFCompoType, YF); } else if ( !MFformat.compare("RealList")) { // use a list of Real. MUST contains an entry for each species in the mixture std::string MFCompoType; pp.query("mixtureFraction.type", MFCompoType); @@ -1335,29 +1339,29 @@ PeleLM::parseComposition(Vector compositionIn, pele::physics::eos::speciesNames(specNames); // For each entry in the user-provided composition, parse name and value - std::string delimiter = ":"; + std::string delimiter = ":"; int specCountIn = compositionIn.size(); for (int i = 0; i < specCountIn; i++ ) { long unsigned sep = compositionIn[i].find(delimiter); if ( sep == std::string::npos ) { Abort("Error parsing '"+compositionIn[i]+"' --> unable to find delimiter :"); - } + } std::string specNameIn = compositionIn[i].substr(0, sep); Real value = std::stod(compositionIn[i].substr(sep+1,compositionIn[i].length())); - int foundIt = 0; + int foundIt = 0; for (int k = 0; k < NUM_SPECIES; k++ ) { if ( specNameIn == specNames[k] ) { compoIn[k] = value; - foundIt = 1; - } - } + foundIt = 1; + } + } if ( !foundIt ) { Abort("Error parsing '"+compositionIn[i]+"' --> unable to match to any species name"); - } + } } // Ensure that it sums to 1.0: - Real sum = 0.0; + Real sum = 0.0; for (int k = 0; k < NUM_SPECIES; k++ ) { sum += compoIn[k]; } @@ -1391,7 +1395,7 @@ void PeleLM::extendSignedDistance(MultiFab *a_signDist, const auto& ebfactory = dynamic_cast(a_signDist->Factory()); const auto& flags = ebfactory.getMultiEBCellFlagFab(); int nGrowFac = flags.nGrow()+1; - + // First set the region far away at the max value we need #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -1404,7 +1408,7 @@ void PeleLM::extendSignedDistance(MultiFab *a_signDist, AMREX_GPU_DEVICE (int i, int j, int k) noexcept { if ( sd_cc(i,j,k) >= maxSignedDist - 1e-12 ) { - const Real* dx = geomdata.CellSize(); + const Real* dx = geomdata.CellSize(); sd_cc(i,j,k) = nGrowFac*dx[0]*a_extendFactor; } }); @@ -1457,7 +1461,7 @@ void PeleLM::extendSignedDistance(MultiFab *a_signDist, } }); } - a_signDist->FillBoundary(geom[0].periodicity()); + a_signDist->FillBoundary(geom[0].periodicity()); } } #endif diff --git a/Source/PeleLM_Index.H b/Source/PeleLM_Index.H index f13219e0..8290abbb 100644 --- a/Source/PeleLM_Index.H +++ b/Source/PeleLM_Index.H @@ -1,9 +1,9 @@ #ifndef PELELM_INDEX_H_ #define PELELM_INDEX_H_ -#define VELX 0 -#define VELY 1 -#define VELZ 2 +#define VELX 0 +#define VELY 1 +#define VELZ 2 #define DENSITY AMREX_SPACEDIM #define FIRSTSPEC DENSITY + 1 #define RHOH FIRSTSPEC + NUM_SPECIES diff --git a/Source/main.cpp b/Source/main.cpp index c2954c87..e62d8a97 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -39,14 +39,14 @@ int main(int argc, char* argv[]) { pelelm.Init(); // Switch between Evolve and UnitTest mode - if ( pelelm.runMode() == "normal" ) { + if ( pelelm.runMode() == "normal" ) { // Advance solution to final time pelelm.Evolve(); } else if ( pelelm.runMode() == "evaluate" ) { - // + // pelelm.Evaluate(); } else {