Line data Source code
1 : #include "Flame.H"
2 : #include "IO/ParmParse.H"
3 : #include "BC/Constant.H"
4 : #include "Numeric/Stencil.H"
5 : #include "IC/Laminate.H"
6 : #include "IC/Constant.H"
7 : #include "IC/PSRead.H"
8 : #include "Numeric/Function.H"
9 : #include "IC/Expression.H"
10 : #include "IC/BMP.H"
11 : #include "IC/PNG.H"
12 : #include "Base/Mechanics.H"
13 : #include "Util/Util.H"
14 : #include "Model/Propellant/Propellant.H"
15 : #include "Model/Propellant/FullFeedback.H"
16 : #include "Model/Propellant/Homogenize.H"
17 :
18 : #include <cmath>
19 :
20 : namespace Integrator
21 : {
22 :
23 3 : Flame::Flame() : Base::Mechanics<model_type>() {}
24 :
25 3 : Flame::Flame(IO::ParmParse& pp) : Flame()
26 : {
27 9 : pp_queryclass(*this);
28 3 : }
29 :
30 :
31 : void
32 3 : Flame::Forbids(IO::ParmParse& pp)
33 : {
34 24 : pp.forbid("pressure.P","use chamber.pressure instead");
35 :
36 24 : pp.forbid("geometry.x_len","This is specified by geometry.prob_lo/hi");
37 24 : pp.forbid("geometry.y_len","This is specified by geometry.prob_lo/hi");
38 24 : pp.forbid("amr.ghost_cells", "This should not be adjustable ");
39 :
40 24 : pp.forbid("pf.gamma","implemented in regression rate law objects now");
41 :
42 24 : pp.forbid("pressure.r_ap", "use regression.powerlaw.r_ap");
43 24 : pp.forbid("pressure.r_htpb", "use regression.powerlaw.r_htpb");
44 24 : pp.forbid("pressure.r_comb", "use regression.powerlaw.r_comb");
45 24 : pp.forbid("pressure.n_ap", "use regression.powerlaw.n_ap");
46 24 : pp.forbid("pressure.n_htpb", "use regression.powerlaw.n_htpb");
47 24 : pp.forbid("pressure.n_comb", "use regression.powerlaw.n_comb");
48 :
49 24 : pp.forbid("thermal.bound", "use thermal.Tref");
50 24 : pp.forbid("thermal.T_fluid", "use thermal.Tfluid (or nothing)");
51 24 : pp.forbid("thermal.m_ap", "use regression.arrhenius.m_ap");
52 24 : pp.forbid("thermal.m_htpb", "use regression.arrhenius.m_htpb");
53 24 : pp.forbid("thermal.E_ap", "use regression.arrhenius.E_ap");
54 24 : pp.forbid("thermal.E_htpb", "use regression.arrhenius.E_htpb");
55 24 : pp.forbid("thermal.modeling_ap", "Old debug variable. Should equal 1 ");
56 24 : pp.forbid("thermal.modeling_htpb", "Old debug variable. Should equal 1");
57 :
58 24 : pp.forbid("pressure.a1", "use propellant.mesoscale.a1 instead");
59 24 : pp.forbid("pressure.a2", "use propellant.mesoscale.a2 instead");
60 24 : pp.forbid("pressure.a3", "use propellant.mesoscale.a3 instead");
61 24 : pp.forbid("pressure.b1", "use propellant.mesoscale.b1 instead");
62 24 : pp.forbid("pressure.b2", "use propellant.mesoscale.b2 instead");
63 24 : pp.forbid("pressure.b3", "use propellant.mesoscale.b3 instead");
64 24 : pp.forbid("pressure.c1", "use propellant.mesoscale.c1 instead");
65 24 : pp.forbid("pressure.mob_ap", "no longer used");
66 24 : pp.forbid("pressure.dependency", "use propellant.mesoscale.arrhenius_dependency");
67 24 : pp.forbid("pressure.h1", "use propellant.mesoscale.h1 instead");
68 24 : pp.forbid("pressure.h2", "use propellant.mesoscale.h2 instead");
69 24 : pp.forbid("thermal.mlocal_ap", "use propellant.mesoscale.mlocal_ap");
70 24 : pp.forbid("thermal.mlocal_comb", "use propellant.mesoscale.mlocal_comb");
71 24 : pp.forbid("thermal.mlocal_htpb", "this actually did **nothing** - it was overridden by a hard code using massfraction.");
72 :
73 24 : pp.forbid("thermal.disperssion1", "use propellant.mesoscale.dispersion1");
74 24 : pp.forbid("thermal.disperssion2", "use propellant.mesoscale.dispersion2");
75 24 : pp.forbid("thermal.disperssion3", "use propellant.mesoscale.dispersion3");
76 :
77 24 : pp.forbid("thermal.rho_ap", "use propellant.XX.rho_ap ");
78 24 : pp.forbid("thermal.rho_htpb","use propellant.XX.rho_htpb ");
79 24 : pp.forbid("thermal.k_ap", "use propellant.XX.k_ap ");
80 24 : pp.forbid("thermal.k_htpb", "use propellant.XX.k_htpb ");
81 24 : pp.forbid("thermal.cp_ap", "use propellant.XX.cp_ap ");
82 21 : pp.forbid("thermal.cp_htpb","use propellant.XX.cp_htpb ");
83 3 : }
84 :
85 :
86 : // [parser]
87 : void
88 3 : Flame::Parse(Flame& value, IO::ParmParse& pp)
89 : {
90 : BL_PROFILE("Integrator::Flame::Flame()");
91 :
92 3 : Forbids(pp);
93 :
94 : // Whether to include extra fields (such as mdot, etc) in the plot output
95 18 : pp_query_default("plot_field",value.plot_field,true);
96 :
97 : //
98 : // PHASE FIELD VARIABLES
99 : //
100 :
101 18 : pp_query_default("pf.eps", value.pf.eps, 0.0); // Burn width thickness
102 18 : pp_query_default("pf.kappa", value.pf.kappa, 0.0); // Interface energy param
103 18 : pp_query_default("pf.lambda", value.pf.lambda, 0.0); // Chemical potential multiplier
104 18 : pp_query_default("pf.w1", value.pf.w1, 0.0); // Unburned rest energy
105 18 : pp_query_default("pf.w12", value.pf.w12, 0.0); // Barrier energy
106 15 : pp_query_default("pf.w0", value.pf.w0, 0.0); // Burned rest energy
107 :
108 : // Boundary conditions for phase field order params
109 6 : pp.select<BC::Constant>("pf.eta.bc", value.bc_eta, 1 );
110 9 : value.RegisterNewFab(value.eta_mf, value.bc_eta, 1, 2, "eta", true);
111 9 : value.RegisterNewFab(value.eta_old_mf, value.bc_eta, 1, 2, "eta_old", false);
112 :
113 : // phase field initial condition
114 6 : pp.select<IC::Laminate,IC::Constant,IC::Expression,IC::BMP,IC::PNG>("pf.eta.ic",value.ic_eta,value.geom);
115 :
116 :
117 :
118 : // Select reduced order model to capture heat feedback
119 : pp.select< Model::Propellant::PowerLaw,
120 : Model::Propellant::FullFeedback,
121 : Model::Propellant::Homogenize>
122 9 : ("propellant",value.propellant);
123 :
124 :
125 : // Whether to use the Thermal Transport Model
126 18 : pp_query_default("thermal.on", value.thermal.on, false);
127 :
128 : // Reference temperature
129 : // Used to set all other reference temperatures by default.
130 15 : pp_query_default("thermal.Tref", value.thermal.Tref, 300.0);
131 :
132 3 : if (value.thermal.on) {
133 :
134 : // Used to change heat flux units
135 6 : pp_query_default("thermal.hc", value.thermal.hc, 1.0);
136 : // System AP mass fraction
137 : //pp_query_default("thermal.massfraction", value.thermal.massfraction, 0.8);
138 :
139 : // Effective fluid temperature
140 5 : pp_query_default("thermal.Tfluid", value.thermal.Tfluid, value.thermal.Tref);
141 :
142 : //Temperature boundary condition
143 2 : pp.select_default<BC::Constant>("thermal.temp.bc", value.bc_temp, 1);
144 :
145 3 : value.RegisterNewFab(value.temp_mf, value.bc_temp, 1, 3, "temp", true);
146 3 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, 1, 3, "temp_old", false);
147 3 : value.RegisterNewFab(value.temps_mf, value.bc_temp, 1, 0, "temps", false);
148 :
149 3 : value.RegisterNewFab(value.mdot_mf, value.bc_temp, 1, 0, "mdot", value.plot_field);
150 3 : value.RegisterNewFab(value.alpha_mf, value.bc_temp, 1, 0, "alpha", value.plot_field);
151 3 : value.RegisterNewFab(value.heatflux_mf, value.bc_temp, 1, 0, "heatflux", value.plot_field);
152 3 : value.RegisterNewFab(value.laser_mf, value.bc_temp, 1, 0, "laser", value.plot_field);
153 :
154 2 : value.RegisterIntegratedVariable(&value.chamber.volume, "volume");
155 2 : value.RegisterIntegratedVariable(&value.chamber.area, "area");
156 2 : value.RegisterIntegratedVariable(&value.chamber.massflux, "mass_flux");
157 :
158 : // laser initial condition
159 2 : pp.select_default<IC::Constant,IC::Expression>("laser.ic",value.ic_laser, value.geom);
160 :
161 : // thermal initial condition
162 3 : pp.select_default<IC::Constant,IC::Expression,IC::BMP,IC::PNG>("temp.ic",value.thermal.ic_temp,value.geom);
163 : }
164 :
165 :
166 : // Constant pressure value
167 18 : pp_query_default("chamber.pressure", value.chamber.pressure, 1.0);
168 :
169 : // Whether to compute the pressure evolution
170 18 : pp_query_default("variable_pressure", value.variable_pressure, 0);
171 :
172 : // Refinement criterion for eta field
173 18 : pp_query_default("amr.refinement_criterion", value.m_refinement_criterion, 0.001);
174 :
175 : // Refinement criterion for temperature field
176 18 : pp_query_default("amr.refinement_criterion_temp", value.t_refinement_criterion, 0.001);
177 :
178 : // Eta value to restrict the refinament for the temperature field
179 18 : pp_query_default("amr.refinament_restriction", value.t_refinement_restriction, 0.1);
180 :
181 : // Refinement criterion for phi field [infinity]
182 18 : pp_query_default("amr.phi_refinement_criterion", value.phi_refinement_criterion, 1.0e100);
183 :
184 : // Minimum allowable threshold for $\eta$
185 15 : pp_query_default("small", value.small, 1.0e-8);
186 :
187 : // Initial condition for $\phi$ field.
188 : pp.select_default<IC::PSRead,IC::Laminate,IC::Expression,IC::Constant,IC::BMP,IC::PNG>
189 6 : ("phi.ic",value.ic_phi,value.geom);
190 :
191 9 : value.RegisterNodalFab(value.phi_mf, 1, 2, "phi", true);
192 :
193 : // Whether to use Neo-hookean Elastic model
194 18 : pp_query_default("elastic.on", value.elastic.on, 0);
195 :
196 : // Body force
197 18 : pp_query_default("elastic.traction", value.elastic.traction, 0.0);
198 :
199 : // Phi refinement criteria
200 18 : pp_query_default("elastic.phirefinement", value.elastic.phirefinement, 1);
201 :
202 15 : pp.queryclass<Base::Mechanics<model_type>>("elastic",value);
203 :
204 3 : if (value.m_type != Type::Disable)
205 : {
206 : // Reference temperature for thermal expansion
207 : // (temperature at which the material is strain-free)
208 0 : pp_query_default("Telastic", value.elastic.Telastic, value.thermal.Tref);
209 : // elastic model of AP
210 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_ap", value.elastic.model_ap);
211 : // elastic model of HTPB
212 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_htpb", value.elastic.model_htpb);
213 :
214 : // Use our current eta field as the psi field for the solver
215 0 : value.psi_on = false;
216 0 : value.solver.setPsi(value.eta_mf);
217 : }
218 :
219 : bool allow_unused;
220 : // Set this to true to allow unused inputs without error.
221 : // (Not recommended.)
222 15 : pp.query_default("allow_unused",allow_unused,false);
223 3 : if (!allow_unused && pp.AnyUnusedInputs())
224 : {
225 0 : Util::Warning(INFO,"The following inputs were specified but not used:");
226 0 : pp.AllUnusedInputs();
227 0 : Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
228 : }
229 3 : }
230 :
231 14 : void Flame::Initialize(int lev)
232 : {
233 : BL_PROFILE("Integrator::Flame::Initialize");
234 14 : Base::Mechanics<model_type>::Initialize(lev);
235 :
236 14 : ic_eta->Initialize(lev, eta_mf);
237 14 : ic_eta->Initialize(lev, eta_old_mf);
238 14 : ic_phi->Initialize(lev, phi_mf);
239 : //ic_phicell->Initialize(lev, phicell_mf);
240 :
241 14 : if (elastic.on) {
242 0 : rhs_mf[lev]->setVal(Set::Vector::Zero());
243 : }
244 14 : if (thermal.on) {
245 8 : if (thermal.ic_temp)
246 : {
247 8 : thermal.ic_temp->Initialize(lev,temp_mf);
248 8 : thermal.ic_temp->Initialize(lev,temp_old_mf);
249 8 : thermal.ic_temp->Initialize(lev,temps_mf);
250 : }
251 : else
252 : {
253 0 : temp_mf[lev]->setVal(thermal.Tref);
254 0 : temp_old_mf[lev]->setVal(thermal.Tref);
255 0 : temps_mf[lev]->setVal(thermal.Tref);
256 : }
257 8 : alpha_mf[lev]->setVal(0.0);
258 8 : mdot_mf[lev]->setVal(0.0);
259 8 : heatflux_mf[lev]->setVal(0.0);
260 8 : thermal.w1 = 0.2 * chamber.pressure + 0.9;
261 8 : ic_laser->Initialize(lev, laser_mf);
262 : }
263 14 : if (variable_pressure) chamber.pressure = 1.0;
264 14 : }
265 :
266 0 : void Flame::UpdateModel(int /*a_step*/, Set::Scalar /*a_time*/)
267 : {
268 0 : if (m_type == Base::Mechanics<model_type>::Type::Disable) return;
269 :
270 0 : for (int lev = 0; lev <= finest_level; ++lev)
271 : {
272 0 : amrex::Box domain = this->geom[lev].Domain();
273 0 : domain.convert(amrex::IntVect::TheNodeVector());
274 0 : const Set::Scalar* DX = geom[lev].CellSize();
275 :
276 0 : phi_mf[lev]->FillBoundary();
277 0 : eta_mf[lev]->FillBoundary();
278 0 : temp_mf[lev]->FillBoundary();
279 :
280 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
281 : {
282 0 : amrex::Box smallbox = mfi.nodaltilebox();
283 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
284 0 : Set::Patch<model_type> model = model_mf.Patch(lev,mfi);
285 0 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
286 0 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
287 0 : Set::Patch<Set::Vector> rhs = rhs_mf.Patch(lev,mfi);
288 :
289 0 : if (elastic.on)
290 : {
291 0 : Set::Patch <const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
292 0 : amrex::ParallelFor(smallbox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
293 : {
294 0 : Set::Vector grad_eta = Numeric::CellGradientOnNode(eta, i, j, k, 0, DX);
295 0 : rhs(i, j, k) = elastic.traction * grad_eta;
296 0 : });
297 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
298 : {
299 0 : Set::Scalar phi_avg = phi(i, j, k, 0);
300 0 : Set::Scalar temp_avg = Numeric::Interpolate::CellToNodeAverage(temp, i, j, k, 0);
301 0 : model_type model_ap = elastic.model_ap;
302 0 : model_ap.F0 -= Set::Matrix::Identity();
303 0 : model_ap.F0 *= (temp_avg - elastic.Telastic);
304 0 : model_ap.F0 += Set::Matrix::Identity();
305 0 : model_type model_htpb = elastic.model_htpb;
306 0 : model_htpb.F0 -= Set::Matrix::Identity();
307 0 : model_htpb.F0 *= (temp_avg - elastic.Telastic);
308 0 : model_htpb.F0 += Set::Matrix::Identity();
309 :
310 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
311 0 : });
312 : }
313 : else
314 : {
315 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
316 : {
317 0 : Set::Scalar phi_avg = Numeric::Interpolate::CellToNodeAverage(phi, i, j, k, 0);
318 : //phi_avg = phi(i,j,k,0);
319 0 : model_type model_ap = elastic.model_ap;
320 0 : model_ap.F0 *= Set::Matrix::Zero();
321 0 : model_type model_htpb = elastic.model_htpb;
322 0 : model_htpb.F0 *= Set::Matrix::Zero();
323 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
324 0 : });
325 : }
326 0 : }
327 0 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
328 :
329 : }
330 : }
331 :
332 818 : void Flame::TimeStepBegin(Set::Scalar a_time, int a_iter)
333 : {
334 : BL_PROFILE("Integrator::Flame::TimeStepBegin");
335 818 : Base::Mechanics<model_type>::TimeStepBegin(a_time, a_iter);
336 818 : if (thermal.on) {
337 90 : for (int lev = 0; lev <= finest_level; ++lev)
338 80 : ic_laser->Initialize(lev, laser_mf, a_time);
339 : }
340 818 : }
341 :
342 818 : void Flame::TimeStepComplete(Set::Scalar /*a_time*/, int /*a_iter*/)
343 : {
344 : BL_PROFILE("Integrator::Flame::TimeStepComplete");
345 818 : if (variable_pressure) {
346 : //Set::Scalar x_len = geom[0].ProbDomain().length(0);
347 : //Set::Scalar y_len = geom[0].ProbDomain().length(1);
348 : // Set::Scalar domain_area = x_len * y_len;
349 0 : Util::Message(INFO, "Mass = ", chamber.massflux);
350 0 : Util::Message(INFO, "Pressure = ", chamber.pressure);
351 : }
352 818 : }
353 :
354 8206 : void Flame::Advance(int lev, Set::Scalar time, Set::Scalar dt)
355 : {
356 : BL_PROFILE("Integrador::Flame::Advance");
357 8206 : Base::Mechanics<model_type>::Advance(lev, time, dt);
358 8206 : const Set::Scalar* DX = geom[lev].CellSize();
359 :
360 8206 : std::swap(eta_old_mf[lev], eta_mf[lev]);
361 :
362 : //
363 : // Chamber pressure update
364 : //
365 8206 : if (variable_pressure) {
366 0 : chamber.pressure = exp(0.00075 * chamber.massflux);
367 0 : if (chamber.pressure > 10.0) {
368 0 : chamber.pressure = 10.0;
369 : }
370 0 : else if (chamber.pressure <= 0.99) {
371 0 : chamber.pressure = 0.99;
372 : }
373 0 : elastic.traction = chamber.pressure;
374 : }
375 :
376 :
377 : //
378 : // Multi-well chemical potential
379 : //
380 : Numeric::Function::Polynomial<4> w( pf.w0,
381 : 0.0,
382 8206 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
383 8206 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
384 8206 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0);
385 8206 : Numeric::Function::Polynomial<3> dw = w.D();
386 :
387 8206 : propellant.set_pressure(chamber.pressure);
388 :
389 30068 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
390 : {
391 21862 : const amrex::Box& bx = mfi.tilebox();
392 : // Phase fields
393 21862 : Set::Patch<Set::Scalar> etanew = eta_mf.Patch(lev,mfi);
394 21862 : Set::Patch<const Set::Scalar> eta = eta_old_mf.Patch(lev,mfi);
395 21862 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
396 : // Heat transfer fields
397 21862 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
398 21862 : Set::Patch<Set::Scalar> alpha = alpha_mf.Patch(lev,mfi);
399 21862 : Set::Patch<Set::Scalar> laser = laser_mf.Patch(lev,mfi);
400 : // Diagnostic fields
401 21862 : Set::Patch<Set::Scalar> mdot = mdot_mf.Patch(lev,mfi);
402 21862 : Set::Patch<Set::Scalar> heatflux = heatflux_mf.Patch(lev,mfi);
403 :
404 :
405 21862 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
406 : {
407 : //
408 : // CALCULATE PHI-AVERAGED QUANTITIES
409 : //
410 59108224 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
411 60856064 : Set::Scalar T = thermal.on ? temp(i,j,k) : NAN;
412 :
413 59108224 : Set::Scalar K = propellant.get_K(phi_avg);
414 :
415 59108224 : Set::Scalar rho = propellant.get_rho(phi_avg);
416 :
417 59108224 : Set::Scalar cp = propellant.get_cp(phi_avg);
418 :
419 : //
420 : // CALCULATE MOBILITY
421 : //
422 59108224 : Set::Scalar L = propellant.get_L( phi_avg, T);
423 :
424 : //
425 : // EVOLVE PHASE FIELD (ETA)
426 : //
427 :
428 59108224 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
429 118216448 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
430 118216448 : etanew(i, j, k) = eta(i, j, k) - L * dt * df_deta;
431 118216448 : if (etanew(i, j, k) <= small) etanew(i, j, k) = small;
432 :
433 59108224 : if (thermal.on)
434 : {
435 : //
436 : // Calculate thermal diffisivity and store for later gradient
437 : //
438 :
439 1747840 : alpha(i, j, k) = K / rho / cp;
440 :
441 : //
442 : // CALCULATE MASS FLUX BASED ON EVOLVING ETA
443 : //
444 :
445 5243520 : mdot(i, j, k) = rho * fabs(eta(i, j, k) - etanew(i, j, k)) / dt;
446 :
447 : //
448 : // CALCULATE HEAT FLUX BASED ON THE CALCULATED MASS FLUX
449 : //
450 :
451 3495680 : Set::Scalar q0 = propellant.get_qdot(mdot(i,j,k), phi_avg);
452 5243520 : heatflux(i,j,k) = ( thermal.hc*q0 + laser(i,j,k) ) / K;
453 :
454 : }
455 :
456 59108224 : });
457 :
458 8206 : } // MFi For loop
459 :
460 :
461 : //
462 : // THERMAL TRANSPORT
463 : //
464 8206 : if (thermal.on)
465 : {
466 2550 : std::swap(temp_old_mf[lev], temp_mf[lev]);
467 :
468 13100 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
469 : {
470 10550 : const amrex::Box& bx = mfi.tilebox();
471 :
472 10550 : Set::Patch<Set::Scalar> tempnew = temp_mf.Patch(lev,mfi);
473 10550 : Set::Patch<const Set::Scalar> temp = temp_old_mf.Patch(lev,mfi);
474 10550 : Set::Patch<const Set::Scalar> alpha = alpha_mf.Patch(lev,mfi);
475 :
476 10550 : Set::Patch<Set::Scalar> temps = temps_mf.Patch(lev,mfi);
477 :
478 :
479 : // Phase field
480 10550 : Set::Patch<Set::Scalar> etanew = (*eta_mf[lev]).array(mfi);
481 10550 : Set::Patch<const Set::Scalar> eta = (*eta_old_mf[lev]).array(mfi);
482 : // Diagnostic fields
483 10550 : Set::Patch<const Set::Scalar> heatflux = heatflux_mf.Patch(lev,mfi);
484 :
485 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
486 : {
487 1747840 : auto sten = Numeric::GetStencil(i, j, k, bx);
488 1747840 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
489 1747840 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
490 1747840 : Set::Scalar lap_temp = Numeric::Laplacian(temp, i, j, k, 0, DX);
491 1747840 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
492 1747840 : Set::Vector grad_alpha = Numeric::Gradient(alpha, i, j, k, 0, DX, sten);
493 1747840 : Set::Scalar dTdt = 0.0;
494 3495680 : dTdt += grad_eta.dot(grad_temp * alpha(i, j, k));
495 3495680 : dTdt += grad_alpha.dot(eta(i, j, k) * grad_temp);
496 3495680 : dTdt += eta(i, j, k) * alpha(i, j, k) * lap_temp;
497 3495680 : dTdt += alpha(i, j, k) * heatflux(i, j, k) * grad_eta_mag;
498 :
499 5243520 : Set::Scalar Tsolid = dTdt + temps(i, j, k) * (etanew(i, j, k) - eta(i, j, k)) / dt;
500 3495680 : temps(i, j, k) = temps(i, j, k) + dt * Tsolid;
501 6991360 : tempnew(i, j, k) = etanew(i, j, k) * temps(i, j, k) + (1.0 - etanew(i, j, k)) * thermal.Tfluid;
502 1747840 : });
503 2550 : }
504 : }
505 :
506 8206 : } //Function
507 :
508 :
509 190 : void Flame::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
510 : {
511 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
512 190 : Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
513 :
514 190 : const Set::Scalar* DX = geom[lev].CellSize();
515 190 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
516 :
517 : // Eta criterion for refinement
518 602 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
519 : {
520 412 : const amrex::Box& bx = mfi.tilebox();
521 412 : amrex::Array4<char> const& tags = a_tags.array(mfi);
522 412 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
523 :
524 412 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
525 : {
526 519920 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
527 541390 : if (gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
528 41936 : tags(i, j, k) = amrex::TagBox::SET;
529 519920 : });
530 190 : }
531 :
532 : // Phi criterion for refinement
533 190 : if (elastic.phirefinement) {
534 602 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
535 : {
536 412 : const amrex::Box& bx = mfi.tilebox();
537 412 : amrex::Array4<char> const& tags = a_tags.array(mfi);
538 412 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
539 :
540 412 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
541 : {
542 519920 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
543 519920 : if (gradphi.lpNorm<2>() * dr >= phi_refinement_criterion)
544 0 : tags(i, j, k) = amrex::TagBox::SET;
545 519920 : });
546 190 : }
547 : }
548 :
549 :
550 : // Thermal criterion for refinement
551 190 : if (thermal.on) {
552 507 : for (amrex::MFIter mfi(*temp_mf[lev], true); mfi.isValid(); ++mfi)
553 : {
554 345 : const amrex::Box& bx = mfi.tilebox();
555 345 : amrex::Array4<char> const& tags = a_tags.array(mfi);
556 345 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
557 345 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
558 345 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
559 : {
560 50928 : Set::Vector tempgrad = Numeric::Gradient(temp, i, j, k, 0, DX);
561 50928 : if (tempgrad.lpNorm<2>() * dr > t_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
562 0 : tags(i, j, k) = amrex::TagBox::SET;
563 50928 : });
564 162 : }
565 : }
566 190 : }
567 :
568 2 : void Flame::Regrid(int lev, Set::Scalar time)
569 : {
570 : BL_PROFILE("Integrator::Flame::Regrid");
571 : //if (lev < finest_level) return;
572 : //phi_mf[lev]->setVal(0.0);
573 2 : ic_phi->Initialize(lev, phi_mf, time);
574 : //ic_phicell->Initialize(lev, phi_mf, time);
575 2 : }
576 :
577 220 : void Flame::Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
578 : const amrex::MFIter& mfi, const amrex::Box& box)
579 : {
580 : BL_PROFILE("Flame::Integrate");
581 220 : const Set::Scalar* DX = geom[amrlev].CellSize();
582 220 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
583 220 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(amrlev,mfi);
584 220 : Set::Patch<const Set::Scalar> mdot = mdot_mf.Patch(amrlev,mfi);
585 220 : if (variable_pressure) {
586 0 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
587 : {
588 0 : chamber.volume += eta(i, j, k, 0) * dv;
589 0 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
590 0 : Set::Scalar normgrad = grad.lpNorm<2>();
591 0 : Set::Scalar da = normgrad * dv;
592 0 : chamber.area += da;
593 :
594 0 : Set::Vector mgrad = Numeric::Gradient(mdot, i, j, k, 0, DX);
595 0 : Set::Scalar mnormgrad = mgrad.lpNorm<2>();
596 0 : Set::Scalar dm = mnormgrad * dv;
597 0 : chamber.massflux += dm;
598 :
599 0 : });
600 : }
601 : else {
602 220 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
603 : {
604 15560 : chamber.volume += eta(i, j, k, 0) * dv;
605 15560 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
606 15560 : Set::Scalar normgrad = grad.lpNorm<2>();
607 15560 : Set::Scalar da = normgrad * dv;
608 15560 : chamber.area += da;
609 15560 : });
610 : }
611 : // time dependent pressure data from experimenta -> p = 0.0954521220950523 * exp(15.289993148880678 * t)
612 220 : }
613 : } // namespace Integrator
614 :
615 :
|
