Line data Source code
1 :
2 : #include "Hydro.H"
3 : #include "AMReX_MultiFab.H"
4 : #include "IO/ParmParse.H"
5 : #include "BC/Constant.H"
6 : #include "BC/Expression.H"
7 : #include "Numeric/Stencil.H"
8 : #include "IC/Constant.H"
9 : #include "IC/Laminate.H"
10 : #include "IC/Expression.H"
11 : #include "IC/BMP.H"
12 : #include "IC/PNG.H"
13 : #include "Solver/Local/Riemann/Roe.H"
14 : #include "Solver/Local/Riemann/HLLE.H"
15 : #include "Solver/Local/Riemann/HLLC.H"
16 : #include "AMReX_TimeIntegrator.H"
17 :
18 : #include "Model/Gas/Gas.H"
19 : #include "Model/Gas/Thermo/Thermo.H"
20 : #include "Model/Gas/Thermo/CpConstant.H"
21 : #include "Model/Gas/Transport/Transport.H"
22 : #include "Model/Gas/Transport/Mixture_Averaged.H"
23 : #include "Model/Gas/EOS/EOS.H"
24 : #include "Model/Gas/EOS/CPG.H"
25 :
26 : namespace Integrator
27 : {
28 :
29 7 : Hydro::Hydro(IO::ParmParse& pp) : Hydro()
30 : {
31 7 : pp_queryclass(*this);
32 7 : }
33 :
34 : void
35 7 : Hydro::Parse(Hydro& value, IO::ParmParse& pp)
36 : {
37 : BL_PROFILE("Integrator::Hydro::Hydro()");
38 : {
39 : // pp.query_default("r_refinement_criterion", value.r_refinement_criterion , 0.01);
40 : // energy-based refinement
41 : // pp.query_default("e_refinement_criterion", value.e_refinement_criterion , 0.01);
42 : // momentum-based refinement
43 : // pp.query_default("m_refinement_criterion", value.m_refinement_criterion , 0.01);
44 :
45 21 : pp.forbid("scheme","use integration.type instead");
46 :
47 : // eta-based refinement
48 14 : pp.query_default("eta_refinement_criterion", value.eta_refinement_criterion , 0.01);
49 : // vorticity-based refinement
50 14 : pp.query_default("omega_refinement_criterion", value.omega_refinement_criterion, 0.01);
51 : // velocity gradient-based refinement
52 14 : pp.query_default("gradu_refinement_criterion", value.gradu_refinement_criterion, 0.01);
53 : // pressure-based refinement
54 14 : pp.query_default("p_refinement_criterion", value.p_refinement_criterion, 1e100);
55 : // density-based refinement
56 21 : pp.query_default("rho_refinement_criterion", value.rho_refinement_criterion, 1e100);
57 :
58 21 : pp_forbid("gamma", "replaced by gas->gamma(...)"); // gamma for gamma law
59 14 : pp_query_required("cfl", value.cfl); // cfl condition
60 21 : pp_query_default("cfl_v", value.cfl_v,1E100); // cfl condition
61 28 : pp_forbid("mu", "replaced with gas->dynamic_viscosity(...)"); // linear viscosity coefficient
62 28 : pp_forbid("Lfactor","replaced with mu");
63 : //pp_query_default("Lfactor", value.Lfactor,1.0); // (to be removed) test factor for viscous source
64 28 : pp_forbid("Pfactor","replaced with mu");
65 : //pp_query_default("Pfactor", value.Pfactor,1.0); // (to be removed) test factor for viscous source
66 28 : pp_forbid("pref", "deprecated - use absolute pressure"); // reference pressure for Roe solver
67 :
68 28 : pp_forbid("rho.bc","--> density.bc");
69 28 : pp_forbid("p.bc","--> pressure.bc");
70 28 : pp_forbid("v.bc", "--> velocity.bc");
71 28 : pp_forbid("pressure.bc","--> energy.bc");
72 21 : pp_forbid("velocity.bc","--> momentum.bc");
73 :
74 : // Boundary condition for density
75 14 : pp.select_default<BC::Constant,BC::Expression>("density.bc",value.density_bc,1);
76 : // Boundary condition for energy
77 14 : pp.select_default<BC::Constant,BC::Expression>("energy.bc",value.energy_bc,1);
78 : // Boundary condition for momentum
79 14 : pp.select_default<BC::Constant,BC::Expression>("momentum.bc",value.momentum_bc,2);
80 :
81 7 : if (!value.managed)
82 : {
83 : // Boundary condition for phase field order parameter
84 21 : pp.select_default<BC::Constant,BC::Expression>("pf.eta.bc",value.eta_bc,1);
85 : }
86 :
87 14 : pp_query_default("small",value.small,1E-8); // small regularization value
88 21 : pp_query_default("cutoff",value.cutoff,-1E100); // cutoff value
89 21 : pp_query_default("lagrange",value.lagrange,0.0); // lagrange no-penetration factor
90 :
91 21 : pp_forbid("roefix","--> solver.roe.entropy_fix"); // Roe solver entropy fix
92 :
93 : }
94 : // Register FabFields:
95 : {
96 7 : int nghost = 1;
97 :
98 7 : if (!value.managed)
99 : {
100 7 : value.eta_mf = new Set::Field<Set::Scalar>();
101 7 : value.eta_old_mf = new Set::Field<Set::Scalar>();
102 21 : value.RegisterNewFab(*value.eta_mf, value.eta_bc, 1, nghost, "eta", true, true);
103 21 : value.RegisterNewFab(*value.eta_old_mf, value.eta_bc, 1, nghost, "eta_old", true, true);
104 : }
105 21 : value.RegisterNewFab(value.etadot_mf, value.eta_bc, 1, nghost, "etadot", true, false);
106 :
107 21 : value.RegisterNewFab(value.density_mf, value.density_bc, 1, nghost, "density", true , true);
108 21 : value.RegisterNewFab(value.density_old_mf, value.density_bc, 1, nghost, "density_old", false, true);
109 :
110 21 : value.RegisterNewFab(value.energy_mf, value.energy_bc, 1, nghost, "energy", true ,true);
111 21 : value.RegisterNewFab(value.energy_old_mf, value.energy_bc, 1, nghost, "energy_old" , false, true);
112 :
113 28 : value.RegisterNewFab(value.momentum_mf, value.momentum_bc, 2, nghost, "momentum", true ,true, {"x","y"});
114 21 : value.RegisterNewFab(value.momentum_old_mf, value.momentum_bc, 2, nghost, "momentum_old", false, true);
115 :
116 21 : value.RegisterNewFab(value.pressure_mf, &value.bc_nothing, 1, nghost, "pressure", true, false);
117 21 : value.RegisterNewFab(value.temperature_mf, &value.bc_nothing, 1, nghost, "temperature", true, false);
118 28 : value.RegisterNewFab(value.velocity_mf, &value.bc_nothing, 2, nghost, "velocity", true, false,{"x","y"});
119 21 : value.RegisterNewFab(value.vorticity_mf, &value.bc_nothing, 1, nghost, "vorticity", true, false);
120 :
121 21 : value.RegisterNewFab(value.m0_mf, &value.bc_nothing, 1, 0, "m0", true, false);
122 28 : value.RegisterNewFab(value.u0_mf, &value.bc_nothing, 2, 0, "u0", true, false, {"x","y"});
123 28 : value.RegisterNewFab(value.q_mf, &value.bc_nothing, 2, 0, "q", true, false, {"x","y"});
124 :
125 28 : value.RegisterNewFab(value.solid.momentum_mf, &value.neumann_bc_D, 2, nghost, "solid.momentum", true, false, {"x","y"});
126 21 : value.RegisterNewFab(value.solid.density_mf, &value.neumann_bc_1, 1, nghost, "solid.density", true, false);
127 21 : value.RegisterNewFab(value.solid.energy_mf, &value.neumann_bc_1, 1, nghost, "solid.energy", true, false);
128 :
129 21 : value.RegisterNewFab(value.Source_mf, &value.bc_nothing, 4, 0, "Source", true, false);
130 :
131 21 : value.RegisterNewFab(value.mass_fraction_mf, &value.bc_nothing, 1, nghost, "mass_fraction", true , true);
132 21 : value.RegisterNewFab(value.mole_fraction_mf, &value.bc_nothing, 1, nghost, "mole_fraction", true , true);
133 21 : value.RegisterNewFab(value.scratch_mf, &value.bc_nothing, 1, nghost, "scratch", false , false);
134 : }
135 :
136 28 : pp_forbid("Velocity.ic.type", "--> velocity.ic.type");
137 28 : pp_forbid("Pressure.ic", "--> pressure.ic");
138 28 : pp_forbid("SolidMomentum.ic", "--> solid.momentum.ic");
139 28 : pp_forbid("SolidDensity.ic.type", "--> solid.density.ic.type");
140 28 : pp_forbid("SolidEnergy.ic.type", "--> solid.energy.ic.type");
141 28 : pp_forbid("Density.ic.type", "--> density.ic.type");
142 28 : pp_forbid("rho_injected.ic.type","no longer using rho_injected use m0 instead");
143 21 : pp.forbid("mdot.ic.type", "replace mdot with u0");
144 :
145 :
146 : // ORDER PARAMETER
147 :
148 7 : if (!value.managed)
149 : {
150 : // eta initial condition
151 21 : pp.select_default<IC::Constant,IC::Laminate,IC::Expression,IC::BMP,IC::PNG>("eta.ic",value.eta_ic,value.geom);
152 : }
153 :
154 : // PRIMITIVE FIELD INITIAL CONDITIONS
155 :
156 : // velocity initial condition
157 14 : pp.select_default<IC::Constant,IC::Expression>("velocity.ic",value.velocity_ic,value.geom);
158 : // solid pressure initial condition
159 14 : pp.select_default<IC::Constant,IC::Expression>("pressure.ic",value.pressure_ic,value.geom);
160 : // density initial condition type
161 14 : pp.select_default<IC::Constant,IC::Expression>("density.ic",value.density_ic,value.geom);
162 :
163 :
164 : // SOLID FIELDS
165 :
166 : // solid momentum initial condition
167 14 : pp.select_default<IC::Constant,IC::Expression>("solid.momentum.ic",value.solid.momentum_ic,value.geom);
168 : // solid density initial condition
169 14 : pp.select_default<IC::Constant,IC::Expression>("solid.density.ic",value.solid.density_ic,value.geom);
170 : // solid energy initial condition
171 14 : pp.select_default<IC::Constant,IC::Expression>("solid.energy.ic",value.solid.energy_ic,value.geom);
172 :
173 :
174 : // DIFFUSE BOUNDARY SOURCES
175 :
176 : // diffuse boundary prescribed mass flux
177 14 : pp.select_default<IC::Constant,IC::Expression>("m0.ic",value.ic_m0,value.geom);
178 : // diffuse boundary prescribed velocity
179 14 : pp.select_default<IC::Constant,IC::Expression>("u0.ic",value.ic_u0,value.geom);
180 : // diffuse boundary prescribed heat flux
181 14 : pp.select_default<IC::Constant,IC::Expression>("q.ic",value.ic_q,value.geom);
182 :
183 : // Riemann solver
184 : pp.select_default< Solver::Local::Riemann::Roe,
185 : Solver::Local::Riemann::HLLE,
186 14 : Solver::Local::Riemann::HLLC>("solver",value.riemannsolver);
187 :
188 : // Gas model (Thermo, Transport, and EOS)
189 14 : pp.queryclass<Model::Gas::Gas>("gas", value.gas);
190 7 : value.nspecies = value.gas.nspecies;
191 7 : std::cout << value.nspecies << "\n";
192 :
193 7 : std::string prescribedflowmode_str;
194 : //
195 28 : pp.query_validate("prescribedflowmode",prescribedflowmode_str,{"absolute","relative"});
196 7 : if (prescribedflowmode_str == "absolute") value.prescribedflowmode = PrescribedFlowMode::Absolute;
197 0 : else if (prescribedflowmode_str == "relative") value.prescribedflowmode = PrescribedFlowMode::Relative;
198 :
199 : // Gravitational acceleration vector
200 21 : pp.queryarr_default("g",value.g,Set::Vector::Zero());
201 :
202 : bool allow_unused;
203 : // Set this to true to allow unused inputs without error.
204 : // (Not recommended.)
205 7 : pp.query_default("allow_unused",allow_unused,false);
206 7 : if (!allow_unused && pp.AnyUnusedInputs(true, false))
207 : {
208 0 : Util::Warning(INFO,"The following inputs were specified but not used:");
209 0 : pp.AllUnusedInputs();
210 0 : Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
211 : }
212 7 : }
213 :
214 :
215 7 : void Hydro::Initialize(int lev)
216 : {
217 : BL_PROFILE("Integrator::Hydro::Initialize");
218 :
219 7 : if (!managed)
220 : {
221 7 : eta_ic ->Initialize(lev, *eta_mf, 0.0);
222 7 : eta_ic ->Initialize(lev, *eta_old_mf, 0.0);
223 : }
224 7 : etadot_mf[lev] ->setVal(0.0);
225 :
226 : //flux_mf[lev] ->setVal(0.0);
227 :
228 7 : velocity_ic ->Initialize(lev, velocity_mf, 0.0);
229 7 : pressure_ic ->Initialize(lev, pressure_mf, 0.0);
230 7 : density_ic ->Initialize(lev, density_mf, 0.0);
231 :
232 7 : density_ic ->Initialize(lev, density_old_mf, 0.0);
233 :
234 7 : solid.density_ic ->Initialize(lev, solid.density_mf, 0.0);
235 7 : solid.momentum_ic->Initialize(lev, solid.momentum_mf, 0.0);
236 7 : solid.energy_ic ->Initialize(lev, solid.energy_mf, 0.0);
237 :
238 7 : ic_m0 ->Initialize(lev, m0_mf, 0.0);
239 7 : ic_u0 ->Initialize(lev, u0_mf, 0.0);
240 7 : ic_q ->Initialize(lev, q_mf, 0.0);
241 :
242 7 : Source_mf[lev] ->setVal(0.0);
243 :
244 7 : if (managed) { if (lev >= (int)mixed.size()) mixed.push_back(false);}
245 7 : else Mix(lev);
246 7 : }
247 :
248 7 : void Hydro::Mix(int lev)
249 : {
250 7 : if (managed && mixed[lev]) return;
251 :
252 14 : for (amrex::MFIter mfi(*velocity_mf[lev], true); mfi.isValid(); ++mfi)
253 : {
254 7 : const amrex::Box& bx = mfi.growntilebox();
255 :
256 7 : Set::Patch<const Set::Scalar> eta_patch = eta_old_mf->Patch(lev,mfi);
257 :
258 7 : Set::Patch<Set::Scalar> v = velocity_mf.Patch(lev,mfi);
259 7 : Set::Patch<Set::Scalar> p = pressure_mf.Patch(lev,mfi);
260 7 : Set::Patch<Set::Scalar> rho = density_mf.Patch(lev,mfi);
261 7 : Set::Patch<Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
262 7 : Set::Patch<Set::Scalar> M = momentum_mf.Patch(lev,mfi);
263 7 : Set::Patch<Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
264 7 : Set::Patch<Set::Scalar> E = energy_mf.Patch(lev,mfi);
265 7 : Set::Patch<Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
266 7 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
267 7 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
268 7 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
269 7 : Set::Patch<Set::Scalar> Y = mass_fraction_mf.Patch(lev,mfi);
270 7 : Set::Patch<Set::Scalar> X = mole_fraction_mf.Patch(lev,mfi);
271 7 : Set::Patch<Set::Scalar> T = temperature_mf.Patch(lev,mfi);
272 :
273 :
274 7 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
275 : {
276 38328 : Set::Scalar eta = invert ? 1.0-eta_patch(i,j,k)*eta_patch(i,j,k) : eta_patch(i,j,k);
277 :
278 : // Initially compute primitives (T,P,u) from given initial conditions
279 : // But from then on, compute them from mixed values to avoid zero T conditions
280 : // Except velocity - keep velocity from fluid values only
281 19164 : gas.ComputeLocalFractions(rho, Y, X, i,j,k); // Get local mole/mass fractions from fluid densities
282 19164 : Set::Scalar density = gas.ComputeD(rho, i, j, k); // If a gas mixture, this will compute the mixture density
283 38328 : T(i,j,k) = gas.ComputeT(p(i,j,k), density, X, i, j, k);
284 76656 : Set::Scalar E_fluid = gas.ComputeE(density, density*v(i,j,k,0), density*v(i,j,k,1), T(i,j,k), X, i, j, k);
285 :
286 : // Mix
287 76656 : M(i, j, k, 0) = (rho(i, j, k)*v(i, j, k, 0))*eta + M_solid(i, j, k, 0)*(1.0-eta);
288 76656 : M(i, j, k, 1) = (rho(i, j, k)*v(i, j, k, 1))*eta + M_solid(i, j, k, 1)*(1.0-eta);
289 38328 : M_old(i, j, k, 0) = M(i, j, k, 0);
290 38328 : M_old(i, j, k, 1) = M(i, j, k, 1);
291 :
292 57492 : rho(i, j, k) = eta * rho(i, j, k) + (1.0 - eta) * rho_solid(i, j, k);
293 38328 : rho_old(i, j, k) = rho(i, j, k);
294 :
295 38328 : E(i, j, k) = E_fluid*eta + E_solid(i,j,k)*(1.0-eta);
296 38328 : E_old(i, j, k) = E(i, j, k);
297 : //Util::Message(INFO,"Energy: ", E(i,j,k), " Pressure: ", p(i,j,k), " Temp: ", T(i,j,k), " Density: ",density, " R: ", gas.R(X,i,j,k), " MW: ", gas.GetMW(X,i,j,k), " Rg: ", Set::Constant::Rg);
298 :
299 : //gas.ComputeLocalFractions(rho, Y, X, i,j,k); // Get local mole/mass fractions from mixed densities
300 : //density = gas.ComputeD(rho, i, j, k);
301 : //T(i, j, k) = gas.ComputeT(density, M(i,j,k,0), M(i,j,k,1), E(i,j,k), T(i,j,k), X, i, j, k);
302 : //p(i, j, k) = gas.ComputeP(density, T(i,j,k), X, i, j, k);
303 : //v(i,j,k,0) = M(i,j,k,0)/density;
304 : //v(i,j,k,1) = M(i,j,k,1)/density;
305 19164 : });
306 : //Util::Abort(INFO);
307 7 : }
308 7 : c_max = 0.0;
309 7 : vx_max = 0.0;
310 7 : vy_max = 0.0;
311 : }
312 :
313 4650 : void Hydro::UpdateEta(int lev, Set::Scalar time)
314 : {
315 32550 : Util::Assert(INFO,TEST(!managed),"Should override this if Hydro is managed!");
316 4650 : eta_ic->Initialize(lev, *eta_mf, time);
317 4650 : }
318 :
319 0 : void Hydro::UpdateFluxes(int /*lev*/, Set::Scalar /*time*/, Set::Scalar /*dt*/)
320 : {
321 0 : Util::Assert(INFO,TEST(!managed),"Should override this if Hydro is managed!");
322 0 : }
323 :
324 4650 : void Hydro::TimeStepBegin(Set::Scalar, int /*iter*/)
325 : {
326 :
327 4650 : }
328 :
329 4650 : void Hydro::TimeStepComplete(Set::Scalar, int lev)
330 : {
331 4650 : if (dynamictimestep.on)
332 0 : Integrator::DynamicTimestep_Update();
333 4650 : return;
334 :
335 : const Set::Scalar* DX = geom[lev].CellSize();
336 :
337 : amrex::ParallelDescriptor::ReduceRealMax(c_max);
338 : amrex::ParallelDescriptor::ReduceRealMax(vx_max);
339 : amrex::ParallelDescriptor::ReduceRealMax(vy_max);
340 :
341 : Set::Scalar new_timestep = cfl / ((c_max + vx_max) / DX[0] + (c_max + vy_max) / DX[1]);
342 :
343 : Util::Assert(INFO, TEST(AMREX_SPACEDIM == 2));
344 :
345 : SetTimestep(new_timestep);
346 : }
347 :
348 4650 : void Hydro::Advance(int lev, Set::Scalar time, Set::Scalar dt)
349 : {
350 :
351 4650 : if (!managed) std::swap(*eta_old_mf, *eta_mf);
352 4650 : std::swap(density_old_mf[lev], density_mf[lev]);
353 4650 : std::swap(momentum_old_mf[lev], momentum_mf[lev]);
354 4650 : std::swap(energy_old_mf[lev], energy_mf[lev]);
355 :
356 : //
357 : // UPDATE ETA AND CALCULATE ETADOT
358 : //
359 :
360 4650 : if (!managed) UpdateEta(lev, time);
361 4650 : if (managed)
362 : {
363 0 : UpdateFluxes(lev,time,dt);
364 0 : Mix(lev);
365 : }
366 9300 : for (amrex::MFIter mfi(*(velocity_mf)[lev], true); mfi.isValid(); ++mfi)
367 : {
368 4650 : const amrex::Box& bx = mfi.growntilebox();
369 4650 : amrex::Array4<const Set::Scalar> const& eta_new = (*(*eta_mf)[lev]).array(mfi);
370 4650 : amrex::Array4<const Set::Scalar> const& eta = (*(*eta_old_mf)[lev]).array(mfi);
371 4650 : amrex::Array4<Set::Scalar> const& etadot = (*etadot_mf[lev]).array(mfi);
372 4650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
373 : {
374 :
375 41203800 : etadot(i, j, k) = (eta_new(i, j, k) - eta(i, j, k)) / dt;
376 13734600 : if (invert) etadot(i,j,k) *= 1.0;
377 :
378 13734600 : });
379 4650 : }
380 :
381 :
382 : //
383 : // DO TIME INTEGRATION (driving the RHS function)
384 : //
385 :
386 : // Organize references to the "new" solution
387 4650 : amrex::Vector<amrex::MultiFab> solution_new;
388 4650 : solution_new.emplace_back(*density_mf[lev].get(),amrex::MakeType::make_alias,0,1);
389 4650 : solution_new.emplace_back(*momentum_mf[lev].get(),amrex::MakeType::make_alias,0,2);
390 4650 : solution_new.emplace_back(*energy_mf[lev].get(),amrex::MakeType::make_alias,0,1);
391 :
392 : // Organize references to the "old" solution
393 4650 : amrex::Vector<amrex::MultiFab> solution_old;
394 4650 : solution_old.emplace_back(*density_old_mf[lev].get(),amrex::MakeType::make_alias,0,1);
395 4650 : solution_old.emplace_back(*momentum_old_mf[lev].get(),amrex::MakeType::make_alias,0,2);
396 4650 : solution_old.emplace_back(*energy_old_mf[lev].get(),amrex::MakeType::make_alias,0,1);
397 :
398 : // Create the time integrator
399 4650 : amrex::TimeIntegrator timeintegrator(solution_new, time);
400 :
401 : // Set the time integrator RHS - in this case, just relay to our current RHS function
402 4650 : timeintegrator.set_rhs([&](amrex::Vector<amrex::MultiFab> & rhs_mf, amrex::Vector<amrex::MultiFab> & solution_mf, const Set::Scalar time)
403 : {
404 5650 : RHS(lev, time,
405 : rhs_mf[0], rhs_mf[1], rhs_mf[2],
406 5650 : solution_mf[0],solution_mf[1],solution_mf[2]);
407 5650 : });
408 :
409 : // Take care of filling boundaries during stages
410 4650 : timeintegrator.set_post_stage_action([&](amrex::Vector<amrex::MultiFab> & stage_mf, Set::Scalar time)
411 : {
412 1000 : density_bc->FillBoundary(stage_mf[0],0,1,time,0);
413 1000 : stage_mf[0].FillBoundary(true);
414 1000 : momentum_bc->FillBoundary(stage_mf[1],0,2,time,0);
415 1000 : stage_mf[1].FillBoundary(true);
416 1000 : energy_bc->FillBoundary(stage_mf[2],0,1,time,0);
417 1000 : stage_mf[2].FillBoundary(true);
418 1000 : });
419 :
420 : // Do the update
421 4650 : timeintegrator.advance(solution_old, solution_new, time, dt);
422 :
423 :
424 : //
425 : // APPLY CUTOFFS AND DO DYNAMIC TIMESTEP CALCULATION
426 : //
427 :
428 4650 : Set::Scalar dt_max = std::numeric_limits<Set::Scalar>::max();
429 9300 : for (amrex::MFIter mfi(*velocity_mf[lev], false); mfi.isValid(); ++mfi)
430 : {
431 4650 : const amrex::Box& bx = mfi.validbox();
432 4650 : const Set::Scalar* DX = geom[lev].CellSize();
433 :
434 4650 : Set::Patch<const Set::Scalar> eta_patch = eta_mf->Patch(lev,mfi);
435 4650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
436 4650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
437 4650 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
438 :
439 4650 : Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
440 4650 : Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
441 4650 : Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
442 :
443 4650 : Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
444 :
445 4650 : Set::Patch<Set::Scalar> u = velocity_mf.Patch(lev,mfi);
446 4650 : Set::Patch<Set::Scalar> Source = Source_mf.Patch(lev,mfi);
447 :
448 4650 : Set::Scalar *dt_max_handle = &dt_max;
449 :
450 4650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
451 : {
452 18278400 : Set::Scalar eta = invert ? 1.0-eta_patch(i,j,k)*eta_patch(i,j,k) : eta_patch(i,j,k);
453 :
454 9139200 : if (eta < cutoff)
455 : {
456 0 : rho_new(i,j,k,0) = rho_solid(i,j,k,0);
457 0 : M_new(i,j,k,0) = M_solid(i,j,k,0);
458 0 : M_new(i,j,k,1) = M_solid(i,j,k,1);
459 0 : E_new(i,j,k,0) = E_solid(i,j,k,0);
460 : }
461 :
462 9139200 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, DX);
463 9139200 : omega(i, j, k) = eta * (gradu(1,0) - gradu(0,1));
464 :
465 9139200 : if (dynamictimestep.on)
466 : {
467 0 : *dt_max_handle = std::fabs(cfl * DX[0] / (u(i,j,k,0)*eta + small));
468 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl * DX[1] / (u(i,j,k,1)*eta + small)));
469 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[0]*DX[0] / (Source(i,j,k,1)+small)));
470 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[1]*DX[1] / (Source(i,j,k,2)+small)));
471 : }
472 9139200 : });
473 4650 : }
474 :
475 :
476 4650 : if (dynamictimestep.on)
477 : {
478 0 : this->DynamicTimestep_SyncTimeStep(lev,dt_max);
479 : }
480 :
481 4650 : }//end Advance
482 :
483 :
484 5650 : void Hydro::RHS(int lev, Set::Scalar /*time*/,
485 : amrex::MultiFab &rho_rhs_mf,
486 : amrex::MultiFab &M_rhs_mf,
487 : amrex::MultiFab &E_rhs_mf,
488 : const amrex::MultiFab &rho_mf,
489 : const amrex::MultiFab &M_mf,
490 : const amrex::MultiFab &E_mf)
491 : {
492 :
493 11300 : for (amrex::MFIter mfi(*(velocity_mf)[lev], true); mfi.isValid(); ++mfi)
494 : {
495 5650 : const amrex::Box& bx = mfi.growntilebox();
496 5650 : amrex::Array4<const Set::Scalar> const& eta_patch = (*(*eta_old_mf)[lev]).array(mfi);
497 :
498 5650 : Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi); // density
499 5650 : Set::Patch<const Set::Scalar> M = M_mf.array(mfi); // momentum
500 5650 : Set::Patch<const Set::Scalar> E = E_mf.array(mfi); // total energy (internal energy + kinetic energy) per unit volume (E/rho = e + 0.5*v^2)
501 :
502 5650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
503 5650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
504 :
505 5650 : Set::Patch<Set::Scalar> scratch = scratch_mf.Patch(lev,mfi);
506 :
507 5650 : Set::Patch<Set::Scalar> v = velocity_mf.Patch(lev,mfi);
508 5650 : Set::Patch<Set::Scalar> p = pressure_mf.Patch(lev,mfi);
509 5650 : Set::Patch<Set::Scalar> T = temperature_mf.Patch(lev,mfi);
510 5650 : Set::Patch<Set::Scalar> Y = mass_fraction_mf.Patch(lev,mfi);
511 5650 : Set::Patch<Set::Scalar> X = mole_fraction_mf.Patch(lev,mfi);
512 :
513 5650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
514 : {
515 33637200 : Set::Scalar eta = invert ? 1.0-eta_patch(i,j,k)*eta_patch(i,j,k) : eta_patch(i,j,k);
516 :
517 : // Compute T and P primitives from mixed values
518 16818600 : Set::Scalar density = gas.ComputeD(rho, i, j, k);
519 84093000 : T(i,j,k) = gas.ComputeT(density, M(i,j,k,0), M(i,j,k,1), E(i,j,k), T(i,j,k), X, i, j, k);
520 33637200 : p(i,j,k) = gas.ComputeP(density, T(i,j,k), X, i, j, k);
521 :
522 : // Compute velocity from fluid values
523 50455800 : scratch(i,j,k) = (rho(i,j,k) - rho_solid(i,j,k)*(1.0 - eta))/(eta + small);
524 16818600 : gas.ComputeLocalFractions(scratch, Y, X, i, j, k);
525 16818600 : Set::Scalar density_fluid = gas.ComputeD(scratch, i, j, k);
526 33637200 : Set::Scalar Mx_fluid = (M(i,j,k,0) - M_solid(i,j,k,0)*(1.0 - eta))/(eta + small);
527 33637200 : Set::Scalar My_fluid = (M(i,j,k,1) - M_solid(i,j,k,1)*(1.0 - eta))/(eta + small);
528 16818600 : v(i,j,k,0) = Mx_fluid/density_fluid;
529 16818600 : v(i,j,k,1) = My_fluid/density_fluid;
530 :
531 16818600 : if (eta < small)
532 : {
533 0 : v(i,j,k,0) *= eta;
534 0 : v(i,j,k,1) *= eta;
535 :
536 : #if AMREX_SPACEDIM == 3
537 0 : v(i,j,k,2) *= eta;
538 : #endif
539 : }
540 16818600 : });
541 5650 : }
542 :
543 5650 : const Set::Scalar* DX = geom[lev].CellSize();
544 5650 : amrex::Box domain = geom[lev].Domain();
545 :
546 11300 : for (amrex::MFIter mfi(*(*eta_mf)[lev], false); mfi.isValid(); ++mfi)
547 : {
548 5650 : const amrex::Box& bx = mfi.validbox();
549 :
550 : // Inputs
551 5650 : Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi);
552 5650 : Set::Patch<const Set::Scalar> E = E_mf.array(mfi);
553 5650 : Set::Patch<const Set::Scalar> M = M_mf.array(mfi);
554 :
555 : // Outputs
556 5650 : Set::Patch<Set::Scalar> rho_rhs = rho_rhs_mf.array(mfi);
557 5650 : Set::Patch<Set::Scalar> M_rhs = M_rhs_mf.array(mfi);
558 5650 : Set::Patch<Set::Scalar> E_rhs = E_rhs_mf.array(mfi);
559 :
560 :
561 : // Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
562 : // Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
563 : // Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
564 :
565 5650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
566 5650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
567 5650 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
568 :
569 5650 : Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
570 :
571 5650 : Set::Patch<const Set::Scalar> eta_patch = eta_old_mf->Patch(lev,mfi);
572 5650 : Set::Patch<const Set::Scalar> etadot = etadot_mf.Patch(lev,mfi);
573 5650 : Set::Patch<const Set::Scalar> velocity = velocity_mf.Patch(lev,mfi);
574 5650 : Set::Patch<const Set::Scalar> T = temperature_mf.Patch(lev,mfi);
575 5650 : Set::Patch<const Set::Scalar> molef = mole_fraction_mf.Patch(lev,mfi);
576 :
577 5650 : Set::Patch<const Set::Scalar> m0 = m0_mf.Patch(lev,mfi);
578 5650 : Set::Patch<const Set::Scalar> q = q_mf.Patch(lev,mfi);
579 5650 : Set::Patch<const Set::Scalar> _u0 = u0_mf.Patch(lev,mfi);
580 :
581 5650 : amrex::Array4<Set::Scalar> const& Source = (*Source_mf[lev]).array(mfi);
582 :
583 5650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
584 : {
585 11187200 : auto sten = Numeric::GetStencil(i, j, k, domain);
586 :
587 22374400 : Set::Scalar eta = invert ? 1.0-eta_patch(i,j,k)*eta_patch(i,j,k) : eta_patch(i,j,k);
588 :
589 : //Diffuse Sources
590 11187200 : Set::Vector grad_eta = Numeric::Gradient(eta_patch, i, j, k, 0, DX);
591 11187200 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
592 11187200 : Set::Matrix hess_eta = Numeric::Hessian(eta_patch, i, j, k, 0, DX);
593 11187200 : if (invert) grad_eta *= -1.0;
594 11187200 : if (invert) hess_eta *= -1.0;
595 :
596 : #if AMREX_SPACEDIM == 2
597 33561600 : Set::Vector u = Set::Vector(velocity(i, j, k, 0), velocity(i, j, k, 1)); // Velocity
598 33561600 : Set::Vector u0 = Set::Vector(_u0(i, j, k, 0), _u0(i, j, k, 1)); // Velocity
599 33561600 : Set::Vector q0 = Set::Vector(q(i,j,k,0), q(i,j,k,1));
600 : #endif
601 :
602 : #if AMREX_SPACEDIM == 3
603 0 : Set::Vector u = Set::Vector(velocity(i, j, k, 0), velocity(i, j, k, 1), velocity(i, j, k, 2)); // Velocity
604 0 : Set::Vector u0 = Set::Vector(_u0(i, j, k, 0), _u0(i, j, k, 1), _u0(i, j, k, 2)); // Velocity
605 0 : Set::Vector q0 = Set::Vector(q(i,j,k,0), q(i,j,k,1), q(i,j,k,2));
606 : #endif
607 :
608 11187200 : Set::Matrix gradM = Numeric::Gradient(M, i, j, k, DX);
609 11187200 : Set::Vector gradrho = Numeric::Gradient(rho,i,j,k,0,DX);
610 11187200 : Set::Matrix hess_rho = Numeric::Hessian(rho,i,j,k,0,DX,sten);
611 22374400 : Set::Matrix gradu = (gradM - u*gradrho.transpose()) / rho(i,j,k);
612 :
613 11187200 : if (prescribedflowmode == PrescribedFlowMode::Relative)
614 : {
615 0 : Set::Vector N = grad_eta / (grad_eta_mag + small);
616 : // Set::Vector T(N(1), -N(0));
617 : // u0 = N * u0(0) + T * u0(1);
618 :
619 : #if AMREX_SPACEDIM == 2
620 0 : Set::Vector T(N(1), -N(0));
621 0 : u0 = N * u0(0) + T * u0(1);
622 : #endif
623 :
624 : #if AMREX_SPACEDIM == 3
625 0 : Set::Vector T;
626 0 : T(0) = N(1);
627 0 : T(1) = -N(0);
628 0 : T(2) = 0;
629 0 : u0 = N*u0(0) + T * u0(1);
630 : // Might not be physcially accurate, need to find how to extend to 3 dimensions
631 : #endif
632 : }
633 :
634 :
635 11187200 : Set::Scalar mdot0 = m0(i,j,k)*grad_eta_mag;
636 11187200 : Set::Vector Pdot0 = Set::Vector::Zero(); // Linear momentum source term
637 11187200 : Set::Scalar qdot0 = q0.dot(grad_eta);
638 :
639 22374400 : Set::Scalar mu = gas.dynamic_viscosity(T(i,j,k), molef, i, j, k);
640 :
641 : // sten is necessary here because sometimes corner ghost
642 : // cells don't get filled
643 11187200 : Set::Matrix3 hess_M = Numeric::Hessian(M,i,j,k,DX);
644 11187200 : Set::Matrix3 hess_u = Set::Matrix3::Zero();
645 33561600 : for (int p = 0; p < 2; p++)
646 67123200 : for (int q = 0; q < 2; q++)
647 134246400 : for (int r = 0; r < 2; r++)
648 : {
649 89497600 : hess_u(r,p,q) =
650 89497600 : (hess_M(r,p,q) - gradu(r,q)*gradrho(p) - gradu(r,p)*gradrho(q) - u(r)*hess_rho(p,q))
651 178995200 : / rho(i,j,k);
652 : }
653 :
654 11187200 : Set::Vector Ldot0 = Set::Vector::Zero();
655 11187200 : Set::Vector div_tau = Set::Vector::Zero();
656 11187200 : Set::Scalar lambda = 0.0; //-2.0/3.0*mu_eff;
657 33561600 : for (int p = 0; p<2; p++)
658 67123200 : for (int q = 0; q<2; q++)
659 134246400 : for (int r = 0; r<2; r++)
660 268492800 : for (int s = 0; s<2; s++)
661 : {
662 178995200 : Ldot0(p) += 0.25 * (mu * ((p==r && q==s) + (p==s && q==r)) + lambda * (p==q && r==s)) * (u(r) - u0(r)) * hess_eta(q, s);
663 536985600 : div_tau(p) += 0.5 * (mu * ((p==r && q==s) + (p==s && q==r)) + lambda * (p==q && r==s)) * (hess_u(r,q,s) + hess_u(s,q,r));
664 :
665 : }
666 :
667 11187200 : Source(i,j, k, 0) = mdot0;
668 11187200 : Source(i,j, k, 1) = Pdot0(0) - Ldot0(0);
669 11187200 : Source(i,j, k, 2) = Pdot0(1) - Ldot0(1);
670 11187200 : Source(i,j, k, 3) = qdot0;// - Ldot0(0)*v(i,j,k,0) - Ldot0(1)*v(i,j,k,1);
671 :
672 : // Lagrange terms to enforce no-penetration
673 11187200 : Source(i,j,k,1) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(0);
674 11187200 : Source(i,j,k,2) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(1);
675 :
676 : //Godunov flux
677 : //states of total fields
678 11187200 : const int X = 0, Y = 1;
679 11187200 : Solver::Local::Riemann::State state_xlo(rho, M, E, i-1, j, k, X);
680 11187200 : Solver::Local::Riemann::State state_x (rho, M, E, i , j, k, X);
681 11187200 : Solver::Local::Riemann::State state_xhi(rho, M, E, i+1, j, k, X);
682 :
683 11187200 : Solver::Local::Riemann::State state_ylo(rho, M, E, i, j-1, k, Y);
684 11187200 : Solver::Local::Riemann::State state_y (rho, M, E, i, j , k, Y);
685 11187200 : Solver::Local::Riemann::State state_yhi(rho, M, E, i, j+1, k, Y);
686 :
687 : //states of solid fields
688 11187200 : Solver::Local::Riemann::State state_xlo_solid(rho_solid, M_solid, E_solid, i-1, j, k, X);
689 11187200 : Solver::Local::Riemann::State state_x_solid (rho_solid, M_solid, E_solid, i , j, k, X);
690 11187200 : Solver::Local::Riemann::State state_xhi_solid(rho_solid, M_solid, E_solid, i+1, j, k, X);
691 :
692 11187200 : Solver::Local::Riemann::State state_ylo_solid(rho_solid, M_solid, E_solid, i, j-1, k, Y);
693 11187200 : Solver::Local::Riemann::State state_y_solid (rho_solid, M_solid, E_solid, i, j , k, Y);
694 11187200 : Solver::Local::Riemann::State state_yhi_solid(rho_solid, M_solid, E_solid, i, j+1, k, Y);
695 :
696 11187200 : Solver::Local::Riemann::State state_xlo_fluid = invert ?
697 0 : (state_xlo - (eta_patch(i-1,j,k))*state_xlo_solid) / (1.0 - eta_patch(i-1,j,k) + small) :
698 33561600 : (state_xlo - (1.0 - eta_patch(i-1,j,k))*state_xlo_solid) / (eta_patch(i-1,j,k) + small);
699 11187200 : Solver::Local::Riemann::State state_x_fluid = invert ?
700 0 : (state_x - (eta_patch(i,j,k) )*state_x_solid ) / (1.0 - eta_patch(i,j,k) + small):
701 33561600 : (state_x - (1.0 - eta_patch(i,j,k) )*state_x_solid ) / (eta_patch(i,j,k) + small);
702 11187200 : Solver::Local::Riemann::State state_xhi_fluid = invert ?
703 0 : (state_xhi - (eta_patch(i+1,j,k))*state_xhi_solid) / (1.0 - eta_patch(i+1,j,k) + small) :
704 33561600 : (state_xhi - (1.0 - eta_patch(i+1,j,k))*state_xhi_solid) / (eta_patch(i+1,j,k) + small);
705 11187200 : Solver::Local::Riemann::State state_ylo_fluid = invert ?
706 0 : (state_ylo - (eta_patch(i,j-1,k))*state_ylo_solid) / (1.0 - eta_patch(i,j-1,k) + small):
707 33561600 : (state_ylo - (1.0 - eta_patch(i,j-1,k))*state_ylo_solid) / (eta_patch(i,j-1,k) + small);
708 11187200 : Solver::Local::Riemann::State state_y_fluid = invert ?
709 0 : (state_y - (eta_patch(i,j,k) )*state_y_solid ) / (1.0 - eta_patch(i,j,k) + small):
710 33561600 : (state_y - (1.0 - eta_patch(i,j,k) )*state_y_solid ) / (eta_patch(i,j,k) + small);
711 11187200 : Solver::Local::Riemann::State state_yhi_fluid = invert ?
712 0 : (state_yhi - (eta_patch(i,j+1,k))*state_yhi_solid) / (1.0 - eta_patch(i,j+1,k) + small):
713 33561600 : (state_yhi - (1.0 - eta_patch(i,j+1,k))*state_yhi_solid) / (eta_patch(i,j+1,k) + small);
714 :
715 11187200 : Solver::Local::Riemann::Flux flux_xlo, flux_ylo, flux_xhi, flux_yhi;
716 :
717 : try
718 : {
719 : //lo interface fluxes
720 11187200 : flux_xlo = riemannsolver->Solve(state_xlo_fluid, state_x_fluid, gas, molef, i, j, k, 0, small) * eta;
721 11187200 : flux_ylo = riemannsolver->Solve(state_ylo_fluid, state_y_fluid, gas, molef, i, j, k, 2, small) * eta;
722 :
723 : //hi interface fluxes
724 11187200 : flux_xhi = riemannsolver->Solve(state_x_fluid, state_xhi_fluid, gas, molef, i, j, k, 1, small) * eta;
725 11187200 : flux_yhi = riemannsolver->Solve(state_y_fluid, state_yhi_fluid, gas, molef, i, j, k, 3, small) * eta;
726 : }
727 0 : catch(...)
728 : {
729 0 : Util::ParallelMessage(INFO,"lev=",lev);
730 0 : Util::ParallelMessage(INFO,"i=",i,"j=",j);
731 0 : Util::Abort(INFO);
732 0 : }
733 :
734 :
735 : Set::Scalar drhof_dt =
736 11187200 : (flux_xlo.mass - flux_xhi.mass) / DX[0] +
737 11187200 : (flux_ylo.mass - flux_yhi.mass) / DX[1] +
738 11187200 : Source(i, j, k, 0);
739 :
740 22374400 : rho_rhs(i,j,k) =
741 : // rho_new(i, j, k) = rho(i, j, k) +
742 : //(
743 11187200 : drhof_dt +
744 : // todo add drhos_dt term if want time-evolving rhos
745 44748800 : etadot(i,j,k) * (rho(i,j,k) - rho_solid(i,j,k)) / (eta + small)
746 : // ) * dt;
747 : ;
748 :
749 :
750 :
751 : Set::Scalar dMxf_dt =
752 11187200 : (flux_xlo.momentum_normal - flux_xhi.momentum_normal ) / DX[0] +
753 22374400 : (flux_ylo.momentum_tangent - flux_yhi.momentum_tangent) / DX[1] +
754 11187200 : div_tau(0) * eta +
755 11187200 : g(0)*rho(i,j,k) +
756 11187200 : Source(i, j, k, 1);
757 :
758 22374400 : M_rhs(i,j,k,0) =
759 : //M_new(i, j, k, 0) = M(i, j, k, 0) +
760 : // (
761 11187200 : dMxf_dt +
762 : // todo add dMs_dt term if want time-evolving Ms
763 44748800 : etadot(i,j,k)*(M(i,j,k,0) - M_solid(i,j,k,0)) / (eta + small)
764 : // ) * dt;
765 : ;
766 :
767 : Set::Scalar dMyf_dt =
768 11187200 : (flux_xlo.momentum_tangent - flux_xhi.momentum_tangent) / DX[0] +
769 22374400 : (flux_ylo.momentum_normal - flux_yhi.momentum_normal ) / DX[1] +
770 11187200 : div_tau(1) * eta +
771 11187200 : g(1)*rho(i,j,k) +
772 11187200 : Source(i, j, k, 2);
773 :
774 22374400 : M_rhs(i,j,k,1) =
775 : //M_new(i, j, k, 1) = M(i, j, k, 1) +
776 : //(
777 11187200 : dMyf_dt +
778 : // todo add dMs_dt term if want time-evolving Ms
779 44748800 : etadot(i,j,k)*(M(i,j,k,1) - M_solid(i,j,k,1)) / (eta+small)
780 : // )*dt;
781 : ;
782 :
783 : Set::Scalar dEf_dt =
784 11187200 : (flux_xlo.energy - flux_xhi.energy) / DX[0] +
785 11187200 : (flux_ylo.energy - flux_yhi.energy) / DX[1] +
786 11187200 : Source(i, j, k, 3);
787 :
788 22374400 : E_rhs(i,j,k) =
789 : // E_new(i, j, k) = E(i, j, k) +
790 : // (
791 11187200 : dEf_dt +
792 : // todo add dEs_dt term if want time-evolving Es
793 44748800 : etadot(i,j,k)*(E(i,j,k) - E_solid(i,j,k)) / (eta+small)
794 : // ) * dt;
795 : ;
796 :
797 : #ifdef AMREX_DEBUG
798 : if ((rho_rhs(i,j,k) != rho_rhs(i,j,k)) ||
799 : (M_rhs(i,j,k,0) != M_rhs(i,j,k,0)) ||
800 : (M_rhs(i,j,k,1) != M_rhs(i,j,k,1)) ||
801 : (E_rhs(i,j,k) != E_rhs(i,j,k)))
802 : {
803 : Util::ParallelMessage(INFO,"rho_rhs=",rho_rhs(i,j,k));
804 : Util::ParallelMessage(INFO,"Mx_rhs=",M_rhs(i,j,k,0));
805 : Util::ParallelMessage(INFO,"Mx_rhs=",M_rhs(i,j,k,1));
806 : Util::ParallelMessage(INFO,"E_rhs=",E_rhs(i,j,k));
807 :
808 : Util::ParallelMessage(INFO,"lev=",lev);
809 : Util::ParallelMessage(INFO,"i=",i," j=",j);
810 : Util::ParallelMessage(INFO,"drhof_dt ",drhof_dt); // dies
811 : Util::ParallelMessage(INFO,"flux_xlo.mass ",flux_xlo.mass);
812 : Util::ParallelMessage(INFO,"flux_xhi.mass ",flux_xhi.mass); // dies, depends on state_xx, state_xhi, state_x_solid, state_xhi_solid, eta, small
813 : Util::ParallelMessage(INFO,"flux_ylo.mass ",flux_ylo.mass);
814 : Util::ParallelMessage(INFO,"flux_xhi.mass ",flux_yhi.mass);
815 : Util::ParallelMessage(INFO,"eta ",eta);
816 : Util::ParallelMessage(INFO,"etadot ",etadot(i,j,k));
817 : Util::ParallelMessage(INFO,"Source ",Source(i,j,k,0));
818 : Util::ParallelMessage(INFO,"state_x ",state_x); // <<<<
819 : Util::ParallelMessage(INFO,"state_y ",state_y);
820 : Util::ParallelMessage(INFO,"state_x_solid ",state_x_solid); // <<<<
821 : Util::ParallelMessage(INFO,"state_y_solid ",state_y_solid);
822 : Util::ParallelMessage(INFO,"state_xhi ",state_xhi); // <<<<
823 : Util::ParallelMessage(INFO,"state_yhi ",state_yhi);
824 : Util::ParallelMessage(INFO,"state_xhi_solid ",state_xhi_solid);
825 : Util::ParallelMessage(INFO,"state_yhi_solids ",state_yhi_solid);
826 : Util::ParallelMessage(INFO,"state_xlo ",state_xlo);
827 : Util::ParallelMessage(INFO,"state_ylo ",state_ylo);
828 : Util::ParallelMessage(INFO,"state_xlo_solid ",state_xlo_solid);
829 : Util::ParallelMessage(INFO,"state_ylo_solid ",state_ylo_solid);
830 :
831 : Util::ParallelMessage(INFO,"Mx_solid ",M_solid(i,j,k,0));
832 : Util::ParallelMessage(INFO,"My_solid ",M_solid(i,j,k,1));
833 : Util::ParallelMessage(INFO,"small ",small);
834 : Util::ParallelMessage(INFO,"Mx ",M(i,j,k,0));
835 : Util::ParallelMessage(INFO,"My ",M(i,j,k,1));
836 : Util::ParallelMessage(INFO,"dMx/dt ",dMxf_dt);
837 : Util::ParallelMessage(INFO,"dMy/dt ",dMyf_dt);
838 :
839 :
840 : Util::Message(INFO,flux_xlo.momentum_tangent);
841 : Util::Message(INFO,flux_xhi.momentum_tangent);
842 : Util::Message(INFO,DX[0]);
843 : Util::Message(INFO,flux_ylo.momentum_normal);
844 : Util::Message(INFO,flux_yhi.momentum_normal);
845 : Util::Message(INFO,DX[1]);
846 : Util::Message(INFO,div_tau);
847 : Util::Message(INFO,Source(i, j, k, 2));
848 :
849 : Util::Message(INFO,hess_eta);
850 : Util::Message(INFO,velocity(i,j,k,0));
851 : Util::Message(INFO,velocity(i,j,k,1));
852 :
853 : Util::Exception(INFO);
854 : }
855 : #endif
856 :
857 :
858 :
859 : // todo - may need to move this for higher order schemes...
860 11187200 : omega(i, j, k) = eta * (gradu(1,0) - gradu(0,1));
861 11187200 : });
862 5650 : }
863 5650 : }
864 :
865 0 : void Hydro::Regrid(int lev, Set::Scalar /* time */)
866 : {
867 : BL_PROFILE("Integrator::Hydro::Regrid");
868 0 : Source_mf[lev]->setVal(0.0);
869 0 : if (lev < finest_level) return;
870 :
871 0 : Util::Message(INFO, "Regridding on level", lev);
872 : }//end regrid
873 :
874 : //void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, Set::Scalar time, int ngrow)
875 0 : void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar, int)
876 : {
877 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
878 :
879 0 : const Set::Scalar* DX = geom[lev].CellSize();
880 0 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
881 :
882 : // Eta criterion for refinement
883 0 : for (amrex::MFIter mfi(*(*eta_mf)[lev], true); mfi.isValid(); ++mfi) {
884 0 : const amrex::Box& bx = mfi.tilebox();
885 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
886 0 : amrex::Array4<const Set::Scalar> const& eta = (*(*eta_mf)[lev]).array(mfi);
887 :
888 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
889 0 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
890 0 : if (grad_eta.lpNorm<2>() * dr * 2 > eta_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
891 0 : });
892 0 : }
893 :
894 : // Vorticity criterion for refinement
895 0 : for (amrex::MFIter mfi(*vorticity_mf[lev], true); mfi.isValid(); ++mfi) {
896 0 : const amrex::Box& bx = mfi.tilebox();
897 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
898 0 : amrex::Array4<const Set::Scalar> const& omega = (*vorticity_mf[lev]).array(mfi);
899 :
900 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
901 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
902 0 : Set::Vector grad_omega = Numeric::Gradient(omega, i, j, k, 0, DX, sten);
903 0 : if (grad_omega.lpNorm<2>() * dr * 2 > omega_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
904 0 : });
905 0 : }
906 :
907 : // Gradu criterion for refinement
908 0 : for (amrex::MFIter mfi(*velocity_mf[lev], true); mfi.isValid(); ++mfi) {
909 0 : const amrex::Box& bx = mfi.tilebox();
910 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
911 0 : amrex::Array4<const Set::Scalar> const& v = (*velocity_mf[lev]).array(mfi);
912 :
913 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
914 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
915 0 : Set::Matrix grad_u = Numeric::Gradient(v, i, j, k, DX, sten);
916 0 : if (grad_u.lpNorm<2>() * dr * 2 > gradu_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
917 0 : });
918 0 : }
919 :
920 : // Pressure criterion for refinement
921 0 : for (amrex::MFIter mfi(*pressure_mf[lev], true); mfi.isValid(); ++mfi) {
922 0 : const amrex::Box& bx = mfi.tilebox();
923 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
924 0 : amrex::Array4<const Set::Scalar> const& p = (*pressure_mf[lev]).array(mfi);
925 :
926 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
927 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
928 0 : Set::Vector grad_p = Numeric::Gradient(p, i, j, k, 0, DX, sten);
929 0 : if (grad_p.lpNorm<2>() * dr * 2 > p_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
930 0 : });
931 0 : }
932 :
933 : // Density criterion for refinement
934 0 : for (amrex::MFIter mfi(*density_mf[lev], true); mfi.isValid(); ++mfi) {
935 0 : const amrex::Box& bx = mfi.tilebox();
936 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
937 0 : amrex::Array4<const Set::Scalar> const& rho = (*density_mf[lev]).array(mfi);
938 :
939 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
940 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
941 0 : Set::Vector grad_rho = Numeric::Gradient(rho, i, j, k, 0, DX, sten);
942 0 : if (grad_rho.lpNorm<2>() * dr * 2 > rho_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
943 0 : });
944 0 : }
945 :
946 0 : }
947 :
948 : }
|
