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