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