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/PointList.H"
8 : #include "IC/PSRead.H"
9 : #include "Numeric/Function.H"
10 : #include "IC/Expression.H"
11 : #include "IC/BMP.H"
12 : #include "IC/PNG.H"
13 : #include "Base/Mechanics.H"
14 : #include "Util/Util.H"
15 : #include "Model/Propellant/Propellant.H"
16 : #include "Model/Propellant/FullFeedback.H"
17 : #include "Model/Propellant/Homogenize.H"
18 : #include <cmath>
19 :
20 : namespace Integrator
21 : {
22 :
23 3 : Flame::Flame() :
24 3 : Base::Mechanics<model_type>() {}
25 :
26 3 : Flame::Flame(IO::ParmParse& pp) : Flame()
27 : {
28 3 : pp_queryclass(*this);
29 3 : }
30 :
31 :
32 : void
33 3 : Flame::Forbids(IO::ParmParse& pp)
34 : {
35 12 : pp.forbid("pressure.P","use chamber.pressure instead");
36 :
37 12 : pp.forbid("geometry.x_len","This is specified by geometry.prob_lo/hi");
38 12 : pp.forbid("geometry.y_len","This is specified by geometry.prob_lo/hi");
39 12 : pp.forbid("amr.ghost_cells", "This should not be adjustable ");
40 :
41 12 : pp.forbid("pf.gamma","use propellant.powerlaw.gamma");
42 :
43 12 : pp.forbid("pressure.r_ap", "use propellant.powerlaw.r_ap");
44 12 : pp.forbid("pressure.r_htpb", "use propellant.powerlaw.r_htpb");
45 12 : pp.forbid("pressure.r_comb", "use propellant.powerlaw.r_comb");
46 12 : pp.forbid("pressure.n_ap", "use propellant.powerlaw.n_ap");
47 12 : pp.forbid("pressure.n_htpb", "use propellant.powerlaw.n_htpb");
48 12 : pp.forbid("pressure.n_comb", "use propellant.powerlaw.n_comb");
49 :
50 12 : pp.forbid("thermal.bound", "use thermal.Tref");
51 12 : pp.forbid("thermal.T_fluid", "use thermal.Tfluid (or nothing)");
52 12 : pp.forbid("thermal.m_ap", "use propellant.fullfeedback.m_ap");
53 12 : pp.forbid("thermal.m_htpb", "use propellant.fullfeedback.m_htpb");
54 12 : pp.forbid("thermal.E_ap", "use propellant.fullfeedback.E_ap");
55 12 : pp.forbid("thermal.E_htpb", "use propellant.fullfeedback.E_htpb");
56 12 : pp.forbid("thermal.modeling_ap", "Old debug variable. Should equal 1 ");
57 12 : pp.forbid("thermal.modeling_htpb", "Old debug variable. Should equal 1");
58 :
59 12 : pp.forbid("pressure.a1", "use propellant.fullfeedback.a1 instead");
60 12 : pp.forbid("pressure.a2", "use propellant.fullfeedback.a2 instead");
61 12 : pp.forbid("pressure.a3", "use propellant.fullfeedback.a3 instead");
62 12 : pp.forbid("pressure.b1", "use propellant.fullfeedback.b1 instead");
63 12 : pp.forbid("pressure.b2", "use propellant.fullfeedback.b2 instead");
64 12 : pp.forbid("pressure.b3", "use propellant.fullfeedback.b3 instead");
65 12 : pp.forbid("pressure.c1", "use propellant.fullfeedback.c1 instead");
66 12 : pp.forbid("pressure.mob_ap", "no longer used");
67 12 : pp.forbid("pressure.dependency", "use propellant.fullfeedback.pressure_dependency");
68 12 : pp.forbid("pressure.h1", "use propellant.homogenize.h1 instead");
69 12 : pp.forbid("pressure.h2", "use propellant.homogenize.h2 instead");
70 12 : pp.forbid("thermal.mlocal_ap", "use propellant.homogenize.mlocal_ap");
71 12 : pp.forbid("thermal.mlocal_comb", "use propellant.homogenize.mlocal_comb");
72 12 : pp.forbid("thermal.mlocal_htpb", "this actually did **nothing** - it was overridden by a hard code using massfraction.");
73 :
74 12 : pp.forbid("thermal.disperssion1", "use propellant.homogenize.dispersion1");
75 12 : pp.forbid("thermal.disperssion2", "use propellant.homogenize.dispersion2");
76 12 : pp.forbid("thermal.disperssion3", "use propellant.homogenize.dispersion3");
77 :
78 12 : pp.forbid("thermal.rho_ap", "use propellant.fullfeedback/homogenize.rho_ap ");
79 12 : pp.forbid("thermal.rho_htpb","use propellant.fullfeedback/homogenize.rho_htpb ");
80 12 : pp.forbid("thermal.k_ap", "use propellant.fullfeedback/homogenize.k_ap ");
81 12 : pp.forbid("thermal.k_htpb", "use propellant.fullfeedback/homogenize.k_htpb ");
82 12 : pp.forbid("thermal.cp_ap", "use propellant.fullfeedback/homogenize.cp_ap ");
83 9 : pp.forbid("thermal.cp_htpb","use propellant.fullfeedback/homogenize.cp_htpb ");
84 3 : }
85 :
86 :
87 : // [parser]
88 : void
89 3 : Flame::Parse(Flame& value, IO::ParmParse& pp)
90 : {
91 : BL_PROFILE("Integrator::Flame::Flame()");
92 :
93 3 : Forbids(pp);
94 :
95 : // Whether to include extra fields (such as mdot, etc) in the plot output
96 6 : pp.query_default("plot_field",value.plot_field,true);
97 :
98 : //
99 : // PHASE FIELD VARIABLES
100 : //
101 :
102 : // Burn width thickness
103 12 : pp.query_default("pf.eps", value.pf.eps, "1.0_m", Unit::Length());
104 : // Interface energy param
105 12 : pp.query_default("pf.kappa", value.pf.kappa, "0.0_J/m^2", Unit::Energy() / Unit::Area());
106 : // Chemical potential multiplier
107 12 : pp.query_default("pf.lambda", value.pf.lambda, "0.0_J/m^2", Unit::Energy()/Unit::Area());
108 : // Unburned rest energy
109 12 : pp.query_default("pf.w1", value.pf.w1, "0.0",Unit::Less());
110 : // Barrier energy
111 12 : pp.query_default("pf.w12", value.pf.w12, "0.0", Unit::Less());
112 : // Burned rest energy
113 12 : pp.query_default("pf.w0", value.pf.w0, "0.0",Unit::Less());
114 :
115 : // Boundary conditions for phase field order params
116 6 : pp.select<BC::Constant>("pf.eta.bc", value.bc_eta, 1 );
117 9 : value.RegisterNewFab(value.eta_mf, value.bc_eta, 1, 2, "eta", true);
118 9 : value.RegisterNewFab(value.eta_old_mf, value.bc_eta, 1, 2, "eta_old", 0);
119 :
120 : // Inital value of eta that doesn't evolve and is used during refiment to set the updated values of eta with voids in the domain.
121 : // Used to fix a bug where duirn refinement, a void won't be updated correctly and would be a square, not a circle
122 9 : value.RegisterNewFab(value.eta_0_mf, value.bc_eta, 1, 2, "eta_0", 0);
123 :
124 : // value.RegisterNewFab(value.eta_mf_frozen, value.bc_eta, 1, 2, "eta_frozen", value.plot_field);
125 :
126 : // phase field initial condition
127 6 : pp.select<IC::Laminate,IC::Constant,IC::Expression,IC::BMP,IC::PNG, IC::PSRead>("pf.eta.ic",value.ic_eta,value.geom);
128 :
129 :
130 : // Select reduced order model to capture heat feedback
131 : pp.select< Model::Propellant::PowerLaw,
132 : Model::Propellant::FullFeedback,
133 : Model::Propellant::Homogenize>
134 6 : ("propellant",value.propellant);
135 :
136 :
137 : // Whether to use the Thermal Transport Model
138 6 : pp_query_default("thermal.on", value.thermal.on, false);
139 :
140 : // Reference temperature
141 : // Used to set all other reference temperatures by default.
142 12 : pp_query_default("thermal.Tref", value.thermal.Tref, "300.0_K",Unit::Temperature());
143 :
144 3 : if (value.thermal.on) {
145 :
146 : // Used to change heat flux units
147 4 : pp_query_default("thermal.hc", value.thermal.hc, "1.0", Unit::Power()/Unit::Area());
148 :
149 : // Effective fluid temperature, temp of the eta = 0 (fluid) region
150 2 : pp_query_default("thermal.Tfluid", value.thermal.Tfluid, value.thermal.Tref);
151 :
152 : // Cutoff value for regression, if T < Tcutoff eta won't evolve/regress
153 4 : pp.query_default("thermal.Tcutoff", value.thermal.Tcutoff, "0.0", Unit::Temperature());
154 :
155 : // Switch time of the improved regridding where eta and the temperature field are both used. It is recommended to make this time ~10x the timestep.
156 : // Before this the refinement is based on the gradient of eta which helps the laser IC start correctly. A regrid is forced when this time is reached.
157 4 : pp.query_default("thermal.end_initial_refine_time", value.thermal.end_initial_refine_time, "0.0", Unit::Time());
158 :
159 : // Inital refinement of the phi field based on phi gradient. After time > end_initial_refine_time stops refining at these phi values.
160 2 : pp.query_default("thermal.phi_refinement_criterion_inital", value.thermal.phi_refinement_criterion_inital, 1.0e100);
161 :
162 : //Temperature boundary condition
163 2 : pp.select_default<BC::Constant>("thermal.temp.bc", value.bc_temp, 1, Unit::Temperature());
164 :
165 3 : value.RegisterNewFab(value.temp_mf, value.bc_temp, 1, 3, "temp", true);
166 3 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, 1, 3, "temp_old", false);
167 3 : value.RegisterNewFab(value.temps_mf, value.bc_temp, 1, 0, "temps", false);
168 :
169 3 : value.RegisterNewFab(value.mdot_mf, value.bc_temp, 1, 0, "mdot", value.plot_field);
170 3 : value.RegisterNewFab(value.alpha_mf, value.bc_temp, 1, 0, "alpha", value.plot_field);
171 3 : value.RegisterNewFab(value.heatflux_mf, value.bc_temp, 1, 0, "heatflux", value.plot_field);
172 3 : value.RegisterNewFab(value.laser_mf, value.bc_temp, 1, 0, "laser", value.plot_field);
173 :
174 2 : value.RegisterIntegratedVariable(&value.chamber.volume, "volume");
175 2 : value.RegisterIntegratedVariable(&value.chamber.area, "area");
176 2 : value.RegisterIntegratedVariable(&value.chamber.massflux, "mass_flux");
177 :
178 3 : value.RegisterNewFab(value.thermal.has_exceeded_Tcutoff, value.bc_temp, 1, 2, "exceeded_Tcutoff", 0); // Used to determine where regression has started
179 :
180 : // laser initial condition
181 : pp.select_default< IC::Constant,
182 : IC::Expression >
183 2 : ("laser.ic",value.ic_laser, value.geom, Unit::Power()/Unit::Area());
184 :
185 : // thermal initial condition
186 : pp.select_default< IC::Constant,
187 : IC::Expression,
188 : IC::BMP,
189 : IC::PNG >
190 3 : ("temp.ic",value.thermal.ic_temp,value.geom, Unit::Temperature());
191 : }
192 :
193 :
194 : // Constant pressure value
195 12 : pp_query_default("chamber.pressure", value.chamber.pressure, "1.0_MPa", Unit::Pressure());
196 :
197 : // Whether to compute the pressure evolution
198 6 : pp_query_default("variable_pressure", value.variable_pressure, false);
199 :
200 : // Refinement criterion for eta field, if thermal is on, cells will only be tagged for refinement if T>0.9*TCutoff,
201 : // and the gradient of eta > m_refinement_criterion at each cell
202 12 : pp_query_default( "amr.refinement_criterion", value.m_refinement_criterion, "0.001",
203 : Unit::Less());
204 :
205 : // Refinement criterion for temperature field
206 12 : pp.query_default( "amr.refinement_criterion_temp", value.t_refinement_criterion, "0.001_K",
207 : Unit::Temperature());
208 :
209 : // Eta value to restrict the refinament for the temperature field
210 12 : pp.query_default( "amr.refinament_restriction", value.t_refinement_restriction, "0.1",
211 : Unit::Less());
212 :
213 : // Refinement criterion for phi field [infinity]
214 6 : pp_query_default("amr.phi_refinement_criterion", value.phi_refinement_criterion, 1.0e100);
215 :
216 : // Minimum allowable threshold for $\eta$
217 6 : pp_query_default("small", value.small, 1.0e-8);
218 :
219 : // Initial condition for $\phi$ field.
220 : pp.select_default<IC::Laminate,IC::Expression,IC::Constant,IC::BMP,IC::PNG, IC::PSRead>
221 6 : ("phi.ic",value.ic_phi,value.geom);
222 :
223 9 : value.RegisterNodalFab(value.phi_mf, 1, 2, "phi", true);
224 :
225 : // Whether to use Neo-hookean Elastic model
226 6 : pp_query_default("elastic.on", value.elastic.on, 0);
227 :
228 : // Body force
229 6 : pp_query_default("elastic.traction", value.elastic.traction, 0.0);
230 :
231 : // Phi refinement criteria
232 6 : pp_query_default("elastic.phirefinement", value.elastic.phirefinement, 1);
233 :
234 6 : pp.queryclass<Base::Mechanics<model_type>>("elastic",value);
235 :
236 3 : if (value.m_type != Type::Disable)
237 : {
238 : // Reference temperature for thermal expansion
239 : // (temperature at which the material is strain-free)
240 0 : pp_query_default("Telastic", value.elastic.Telastic, value.thermal.Tref);
241 : // elastic model of AP
242 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_ap", value.elastic.model_ap);
243 : // elastic model of HTPB
244 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_htpb", value.elastic.model_htpb);
245 :
246 : // eta cutoff value to stop applying the traction force and/or chamber pressure to the RHS.
247 : // Below this vlaue the RHS is set to 0.0
248 0 : pp.query_default("etacutoff",value.elastic.etacutoff, "0", Unit::Less());
249 :
250 : // Boolean value, when not 0 the chamber.pressure value is applied at the diffuse interface where t>Tcutoff
251 0 : pp.query_default("apply_chamber_pressure",value.elastic.apply_chamber_pressure, "0", Unit::Less());
252 :
253 : // Use our current eta field as the psi field for the solver
254 0 : value.psi_on = false;
255 0 : value.solver.setPsi(value.eta_mf);
256 : }
257 :
258 : bool allow_unused;
259 : // Set this to true to allow unused inputs without error.
260 : // (Not recommended.)
261 3 : pp.query_default("allow_unused",allow_unused,false);
262 3 : if (!allow_unused && pp.AnyUnusedInputs())
263 : {
264 0 : Util::Warning(INFO,"The following inputs were specified but not used:");
265 0 : pp.AllUnusedInputs();
266 0 : Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
267 : }
268 3 : }
269 :
270 14 : void Flame::Initialize(int lev)
271 : {
272 : BL_PROFILE("Integrator::Flame::Initialize");
273 14 : Base::Mechanics<model_type>::Initialize(lev);
274 :
275 14 : ic_eta->Initialize(lev, eta_mf);
276 14 : ic_eta->Initialize(lev, eta_old_mf);
277 14 : ic_phi->Initialize(lev, phi_mf);
278 : //ic_phicell->Initialize(lev, phicell_mf);
279 :
280 14 : if (elastic.on) {
281 0 : rhs_mf[lev]->setVal(Set::Vector::Zero());
282 : }
283 14 : if (thermal.on) {
284 8 : if (thermal.ic_temp)
285 : {
286 8 : thermal.ic_temp->Initialize(lev,temp_mf);
287 8 : thermal.ic_temp->Initialize(lev,temp_old_mf);
288 8 : thermal.ic_temp->Initialize(lev,temps_mf);
289 : }
290 : else
291 : {
292 0 : temp_mf[lev]->setVal(thermal.Tref);
293 0 : temp_old_mf[lev]->setVal(thermal.Tref);
294 0 : temps_mf[lev]->setVal(thermal.Tref);
295 : }
296 8 : alpha_mf[lev]->setVal(0.0);
297 8 : mdot_mf[lev]->setVal(0.0);
298 8 : heatflux_mf[lev]->setVal(0.0);
299 8 : ic_laser->Initialize(lev, laser_mf);
300 : }
301 14 : if (variable_pressure) chamber.pressure = 1.0;
302 14 : }
303 :
304 0 : void Flame::UpdateModel(int /*a_step*/, Set::Scalar /*a_time*/)
305 : {
306 0 : if (m_type == Base::Mechanics<model_type>::Type::Disable) return;
307 :
308 0 : for (int lev = 0; lev <= finest_level; ++lev)
309 : {
310 0 : amrex::Box domain = this->geom[lev].Domain();
311 0 : domain.convert(amrex::IntVect::TheNodeVector());
312 0 : const Set::Scalar* DX = geom[lev].CellSize();
313 :
314 0 : phi_mf[lev]->FillBoundary();
315 0 : eta_mf[lev]->FillBoundary();
316 0 : temp_mf[lev]->FillBoundary();
317 :
318 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
319 : {
320 0 : amrex::Box smallbox = mfi.nodaltilebox();
321 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
322 0 : Set::Patch<model_type> model = model_mf.Patch(lev,mfi);
323 0 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
324 0 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
325 0 : Set::Patch<Set::Vector> rhs = rhs_mf.Patch(lev,mfi);
326 0 : Set::Scalar Tcutoff = thermal.Tcutoff;
327 :
328 0 : if (elastic.on)
329 : {
330 0 : Set::Patch <const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
331 0 : amrex::ParallelFor(smallbox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
332 :
333 : {
334 0 : Set::Vector grad_eta = Numeric::CellGradientOnNode(eta, i, j, k, 0, DX);
335 :
336 0 : if (temp(i,j,k) > Tcutoff && eta(i,j,k) > elastic.etacutoff && elastic.apply_chamber_pressure)
337 : {
338 0 : rhs(i, j, k) = (elastic.traction) * grad_eta - chamber.pressure*grad_eta;
339 : // std::cout << "Applying chamber pressure" << std::endl;
340 : }
341 0 : else if (temp(i,j,k) > Tcutoff && eta(i,j,k) > elastic.etacutoff && !elastic.apply_chamber_pressure)
342 : {
343 0 : rhs(i, j, k) = (elastic.traction) * grad_eta;
344 : // std::cout << "Applying traction" << std::endl;
345 : }
346 : else
347 0 : rhs(i, j, k) = 0.0 * grad_eta;
348 : // std::cout << "Applying neither" << std::endl;
349 :
350 0 : });
351 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
352 : {
353 0 : Set::Scalar phi_avg = phi(i, j, k, 0);
354 0 : Set::Scalar temp_avg = Numeric::Interpolate::CellToNodeAverage(temp, i, j, k, 0);
355 0 : model_type model_ap = elastic.model_ap;
356 0 : model_ap.F0 -= Set::Matrix::Identity();
357 0 : model_ap.F0 *= (temp_avg - elastic.Telastic);
358 0 : model_ap.F0 += Set::Matrix::Identity();
359 0 : model_type model_htpb = elastic.model_htpb;
360 0 : model_htpb.F0 -= Set::Matrix::Identity();
361 0 : model_htpb.F0 *= (temp_avg - elastic.Telastic);
362 0 : model_htpb.F0 += Set::Matrix::Identity();
363 :
364 0 : model(i, j, k) = (model_ap * phi_avg + model_htpb * (1. - phi_avg));
365 0 : });
366 : }
367 : else
368 : {
369 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
370 : {
371 0 : Set::Scalar phi_avg = Numeric::Interpolate::CellToNodeAverage(phi, i, j, k, 0);
372 : //phi_avg = phi(i,j,k,0);
373 0 : model_type model_ap = elastic.model_ap;
374 0 : model_ap.F0 *= Set::Matrix::Zero();
375 0 : model_type model_htpb = elastic.model_htpb;
376 0 : model_htpb.F0 *= Set::Matrix::Zero();
377 0 : model(i, j, k) = (model_ap * phi_avg + model_htpb * (1. - phi_avg));
378 0 : });
379 : }
380 0 : }
381 0 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
382 :
383 : }
384 : }
385 :
386 818 : void Flame::TimeStepBegin(Set::Scalar a_time, int a_iter)
387 : {
388 : BL_PROFILE("Integrator::Flame::TimeStepBegin");
389 818 : Base::Mechanics<model_type>::TimeStepBegin(a_time, a_iter);
390 818 : if (thermal.on) {
391 90 : for (int lev = 0; lev <= finest_level; ++lev)
392 80 : ic_laser->Initialize(lev, laser_mf, a_time);
393 : }
394 :
395 818 : if (a_time > thermal.end_initial_refine_time)
396 : {
397 9 : if (!end_initial_refine) {
398 9 : for (int lev = 0; lev <= finest_level; ++lev)
399 8 : Flame::Regrid(lev, a_time);
400 1 : end_initial_refine = 1;
401 : }
402 :
403 9 : prev_finest_ba = grids[finest_level];
404 9 : prev_finest_level = finest_level;
405 : }
406 818 : }
407 :
408 818 : void Flame::TimeStepComplete(Set::Scalar /*a_time*/, int /*a_iter*/)
409 : {
410 : BL_PROFILE("Integrator::Flame::TimeStepComplete");
411 818 : if (variable_pressure) {
412 : //Set::Scalar x_len = geom[0].ProbDomain().length(0);
413 : //Set::Scalar y_len = geom[0].ProbDomain().length(1);
414 : // Set::Scalar domain_area = x_len * y_len;
415 0 : Util::Message(INFO, "Mass = ", chamber.massflux);
416 0 : Util::Message(INFO, "Pressure = ", chamber.pressure);
417 : }
418 818 : }
419 :
420 8206 : void Flame::Advance(int lev, Set::Scalar time, Set::Scalar dt)
421 : {
422 : BL_PROFILE("Integrador::Flame::Advance");
423 8206 : Base::Mechanics<model_type>::Advance(lev, time, dt);
424 8206 : const Set::Scalar* DX = geom[lev].CellSize();
425 :
426 8206 : std::swap(eta_old_mf[lev], eta_mf[lev]);
427 :
428 : //
429 : // Chamber pressure update
430 : //
431 8206 : if (variable_pressure) {
432 0 : chamber.pressure = exp(0.00075 * chamber.massflux);
433 0 : if (chamber.pressure > 10.0) {
434 0 : chamber.pressure = 10.0;
435 : }
436 0 : else if (chamber.pressure <= 0.99) {
437 0 : chamber.pressure = 0.99;
438 : }
439 0 : elastic.traction = chamber.pressure;
440 : }
441 :
442 :
443 : //
444 : // Multi-well chemical potential
445 : //
446 : Numeric::Function::Polynomial<4> w( pf.w0,
447 : 0.0,
448 8206 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
449 8206 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
450 8206 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0);
451 8206 : Numeric::Function::Polynomial<3> dw = w.D();
452 :
453 8206 : propellant.set_pressure(chamber.pressure);
454 :
455 30068 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
456 : {
457 21862 : const amrex::Box& bx = mfi.tilebox();
458 : // Phase fields
459 21862 : Set::Patch<Set::Scalar> etanew = eta_mf.Patch(lev,mfi);
460 21862 : Set::Patch<const Set::Scalar> eta = eta_old_mf.Patch(lev,mfi);
461 21862 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
462 : // Heat transfer fields
463 21862 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
464 21862 : Set::Patch<Set::Scalar> alpha = alpha_mf.Patch(lev,mfi);
465 21862 : Set::Patch<Set::Scalar> laser = laser_mf.Patch(lev,mfi);
466 : // Diagnostic fields
467 21862 : Set::Patch<Set::Scalar> mdot = mdot_mf.Patch(lev,mfi);
468 21862 : Set::Patch<Set::Scalar> heatflux = heatflux_mf.Patch(lev,mfi);
469 :
470 21862 : Set::Patch<Set::Scalar> exceeded_Tcutoff = thermal.has_exceeded_Tcutoff.Patch(lev, mfi);
471 21862 : Set::Scalar Tcutoff = thermal.Tcutoff;
472 :
473 :
474 21862 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
475 : {
476 : //
477 : // CALCULATE PHI-AVERAGED QUANTITIES
478 : //
479 59108224 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
480 60856064 : Set::Scalar T = thermal.on ? temp(i,j,k) : NAN;
481 :
482 59108224 : Set::Scalar K = propellant.get_K(phi_avg);
483 :
484 59108224 : Set::Scalar rho = propellant.get_rho(phi_avg);
485 :
486 59108224 : Set::Scalar cp = propellant.get_cp(phi_avg);
487 :
488 : //
489 : // CALCULATE MOBILITY
490 : //
491 59108224 : Set::Scalar L = propellant.get_L( phi_avg, T);
492 :
493 : //
494 : // EVOLVE PHASE FIELD (ETA)
495 : //
496 :
497 59108224 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
498 118216448 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
499 :
500 59108224 : if (df_deta < 0) {
501 : // Prevent eta from increasing/healing. A bug was found where if the diffuse thickness was too large compared to a void
502 : // (region of eta = 0), eta would heal/increase in a non-physcial way, this statement stops that behavior
503 196236 : df_deta = 0.0;
504 : }
505 59108224 : if (thermal.on && T < thermal.Tcutoff) {
506 : // If the temperature is lower then the cutoff temperature don't evolve the eta field
507 0 : df_deta = 0.0;
508 : }
509 118216448 : etanew(i, j, k) = eta(i, j, k) - L * dt * df_deta;
510 :
511 118216448 : if (etanew(i, j, k) <= small) etanew(i, j, k) = small;
512 :
513 59108224 : if (thermal.on)
514 : {
515 : //
516 : // Calculate thermal diffisivity and store for later gradient
517 : //
518 :
519 1747840 : alpha(i, j, k) = K / rho / cp;
520 :
521 : //
522 : // CALCULATE MASS FLUX BASED ON EVOLVING ETA
523 : //
524 :
525 5243520 : mdot(i, j, k) = rho * fabs(eta(i, j, k) - etanew(i, j, k)) / dt;
526 :
527 : //
528 : // CALCULATE HEAT FLUX BASED ON THE CALCULATED MASS FLUX
529 : //
530 :
531 3495680 : Set::Scalar q0 = propellant.get_qdot(mdot(i,j,k), phi_avg);
532 3495680 : heatflux(i,j,k) = ( thermal.hc*q0 + laser(i,j,k) ) / K;
533 :
534 3495680 : if (temp(i,j,k) > Tcutoff)
535 : {
536 3495680 : exceeded_Tcutoff(i,j,k) = 1;
537 : }
538 :
539 : }
540 :
541 59108224 : });
542 :
543 8206 : } // MFi For loop
544 :
545 :
546 : //
547 : // THERMAL TRANSPORT
548 : //
549 8206 : if (thermal.on)
550 : {
551 2550 : std::swap(temp_old_mf[lev], temp_mf[lev]);
552 :
553 13100 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
554 : {
555 10550 : const amrex::Box& bx = mfi.tilebox();
556 :
557 10550 : Set::Patch<Set::Scalar> tempnew = temp_mf.Patch(lev,mfi);
558 10550 : Set::Patch<const Set::Scalar> temp = temp_old_mf.Patch(lev,mfi);
559 10550 : Set::Patch<const Set::Scalar> alpha = alpha_mf.Patch(lev,mfi);
560 :
561 10550 : Set::Patch<Set::Scalar> temps = temps_mf.Patch(lev,mfi);
562 :
563 :
564 : // Phase field
565 10550 : Set::Patch<Set::Scalar> etanew = (*eta_mf[lev]).array(mfi);
566 10550 : Set::Patch<const Set::Scalar> eta = (*eta_old_mf[lev]).array(mfi);
567 : // Diagnostic fields
568 10550 : Set::Patch<const Set::Scalar> heatflux = heatflux_mf.Patch(lev,mfi);
569 :
570 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
571 : {
572 1747840 : auto sten = Numeric::GetStencil(i, j, k, bx);
573 1747840 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
574 1747840 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
575 1747840 : Set::Scalar lap_temp = Numeric::Laplacian(temp, i, j, k, 0, DX);
576 1747840 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
577 1747840 : Set::Vector grad_alpha = Numeric::Gradient(alpha, i, j, k, 0, DX, sten);
578 1747840 : Set::Scalar dTdt = 0.0;
579 3495680 : dTdt += grad_eta.dot(grad_temp * alpha(i, j, k));
580 3495680 : dTdt += grad_alpha.dot(eta(i, j, k) * grad_temp);
581 3495680 : dTdt += eta(i, j, k) * alpha(i, j, k) * lap_temp;
582 3495680 : dTdt += alpha(i, j, k) * heatflux(i, j, k) * grad_eta_mag;
583 :
584 5243520 : Set::Scalar Tsolid = dTdt + temps(i, j, k) * (etanew(i, j, k) - eta(i, j, k)) / dt;
585 3495680 : temps(i, j, k) = temps(i, j, k) + dt * Tsolid;
586 6991360 : tempnew(i, j, k) = etanew(i, j, k) * temps(i, j, k) + (1.0 - etanew(i, j, k)) * thermal.Tfluid;
587 1747840 : });
588 2550 : }
589 : }
590 :
591 8206 : } //Function
592 :
593 :
594 190 : void Flame::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
595 : {
596 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
597 190 : Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
598 :
599 190 : const Set::Scalar* DX = geom[lev].CellSize();
600 190 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
601 :
602 : // Eta criterion for refinement
603 602 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
604 : {
605 412 : const amrex::Box& bx = mfi.tilebox();
606 412 : amrex::Array4<char> const& tags = a_tags.array(mfi);
607 412 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
608 412 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
609 :
610 412 : if (thermal.on) {
611 345 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
612 : {
613 50928 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
614 76356 : if (gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion && eta(i, j, k) >= t_refinement_restriction && temp(i,j,k) > thermal.Tcutoff*0.9)
615 25428 : tags(i, j, k) = amrex::TagBox::SET;
616 50928 : });
617 :
618 : } else {
619 67 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
620 : {
621 468992 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
622 477748 : if (gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
623 16508 : tags(i, j, k) = amrex::TagBox::SET;
624 468992 : });
625 : }
626 190 : }
627 :
628 : // Phi criterion for refinement
629 190 : if (elastic.phirefinement) {
630 602 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
631 : {
632 412 : const amrex::Box& bx = mfi.tilebox();
633 412 : amrex::Array4<char> const& tags = a_tags.array(mfi);
634 412 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
635 :
636 412 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
637 : {
638 519920 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
639 519920 : if (gradphi.lpNorm<2>() * dr >= phi_refinement_criterion)
640 0 : tags(i, j, k) = amrex::TagBox::SET;
641 519920 : });
642 190 : }
643 : }
644 :
645 :
646 : // Thermal criterion for refinement
647 190 : if (thermal.on) {
648 507 : for (amrex::MFIter mfi(*temp_mf[lev], true); mfi.isValid(); ++mfi)
649 : {
650 345 : const amrex::Box& bx = mfi.tilebox();
651 345 : amrex::Array4<char> const& tags = a_tags.array(mfi);
652 345 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev,mfi);
653 345 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
654 345 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
655 : {
656 50928 : Set::Vector tempgrad = Numeric::Gradient(temp, i, j, k, 0, DX);
657 50928 : if (tempgrad.lpNorm<2>() * dr > t_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
658 0 : tags(i, j, k) = amrex::TagBox::SET;
659 50928 : });
660 162 : }
661 : }
662 :
663 : // Refine at start
664 602 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
665 : {
666 412 : const amrex::Box& bx = mfi.tilebox();
667 412 : amrex::Array4<char> const& tags = a_tags.array(mfi);
668 412 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
669 412 : Set::Patch<const Set::Scalar> phi = phi_mf.Patch(lev,mfi);
670 :
671 412 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
672 : {
673 519920 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
674 519920 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
675 519920 : if ((gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion || gradphi.lpNorm<2>() * dr >= thermal.phi_refinement_criterion_inital) && time < thermal.end_initial_refine_time)
676 0 : tags(i, j, k) = amrex::TagBox::SET;
677 519920 : });
678 190 : }
679 190 : }
680 :
681 10 : void Flame::Regrid(int lev, Set::Scalar time)
682 : {
683 : BL_PROFILE("Integrator::Flame::Regrid");
684 :
685 10 : ic_phi->Initialize(lev, phi_mf, time);
686 10 : ic_eta->Initialize(lev, eta_0_mf, time);
687 :
688 10 : if (thermal.on) {
689 : /*
690 : This regrid function works by using the "has_exceeded_Tcutoff" field. If the temperature in a cell is greater than Tcutoff,
691 : eta will change and when regridding won't use the initial eta field. If T < T_cutoff, when regriding happens it applies the inital
692 : eta field condition. This gives at leat a 4x speed improvement in 2D when doing regression with voids. This is because orgionally
693 : there was a bug where when regridding, the orgional eta field wouldn't be applied, so there would be "squares" of voids instead of
694 : circles/spheres when using .xyzr files as the inital condition.
695 : */
696 24 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
697 : {
698 16 : const amrex::Box &bx = mfi.tilebox();
699 16 : Set::Patch<Set::Scalar> eta = eta_mf.Patch(lev, mfi);
700 16 : Set::Patch<const Set::Scalar> temp = temp_mf.Patch(lev, mfi);
701 16 : Set::Patch<const Set::Scalar> eta_0 = eta_0_mf.Patch(lev, mfi);
702 16 : Set::Patch<Set::Scalar> exceeded_Tcutoff = thermal.has_exceeded_Tcutoff.Patch(lev, mfi);
703 :
704 16 : Set::Scalar Tcutoff = thermal.Tcutoff;
705 :
706 16 : amrex::BoxList boxes_to_update;
707 16 : if (lev == finest_level && prev_finest_level == finest_level)
708 0 : boxes_to_update = amrex::complementIn(bx, prev_finest_ba).boxList();
709 : else
710 16 : boxes_to_update.push_back(bx);
711 :
712 32 : for (const amrex::Box &box : boxes_to_update)
713 16 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
714 : {
715 :
716 4128 : if (!exceeded_Tcutoff(i,j,k) && temp(i,j,k) < Tcutoff)
717 : {
718 0 : eta(i, j, k) = eta_0(i, j, k);
719 : }
720 2064 : });
721 24 : }
722 :
723 8 : if (lev == finest_level)
724 : {
725 1 : prev_finest_ba = grids[finest_level];
726 1 : prev_finest_level = finest_level;
727 : }
728 : }
729 10 : }
730 :
731 220 : void Flame::Integrate(int amrlev, Set::Scalar time, int /*step*/,
732 : const amrex::MFIter& mfi, const amrex::Box& box)
733 : {
734 : BL_PROFILE("Flame::Integrate");
735 :
736 220 : Base::Mechanics<model_type>::Integrate(amrlev,time,timestep,mfi,box);
737 :
738 220 : const Set::Scalar* DX = geom[amrlev].CellSize();
739 220 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
740 220 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(amrlev,mfi);
741 220 : Set::Patch<const Set::Scalar> mdot = mdot_mf.Patch(amrlev,mfi);
742 220 : if (variable_pressure) {
743 0 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
744 : {
745 0 : chamber.volume += eta(i, j, k, 0) * dv;
746 0 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
747 0 : Set::Scalar normgrad = grad.lpNorm<2>();
748 0 : Set::Scalar da = normgrad * dv;
749 0 : chamber.area += da;
750 :
751 0 : Set::Vector mgrad = Numeric::Gradient(mdot, i, j, k, 0, DX);
752 0 : Set::Scalar mnormgrad = mgrad.lpNorm<2>();
753 0 : Set::Scalar dm = mnormgrad * dv;
754 0 : chamber.massflux += dm;
755 :
756 0 : });
757 : }
758 : else {
759 220 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
760 : {
761 15560 : chamber.volume += eta(i, j, k, 0) * dv;
762 15560 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
763 15560 : Set::Scalar normgrad = grad.lpNorm<2>();
764 15560 : Set::Scalar da = normgrad * dv;
765 15560 : chamber.area += da;
766 15560 : });
767 : }
768 : // time dependent pressure data from experimenta -> p = 0.0954521220950523 * exp(15.289993148880678 * t)
769 220 : }
770 : } // namespace Integrator
771 :
|
