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