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 "Base/Mechanics.H"
12 :
13 : #include <cmath>
14 :
15 : namespace Integrator
16 : {
17 :
18 3 : Flame::Flame() : Base::Mechanics<model_type>() {}
19 :
20 3 : Flame::Flame(IO::ParmParse& pp) : Flame()
21 : {
22 3 : pp_queryclass(*this);
23 3 : }
24 :
25 : // [parser]
26 : void
27 3 : Flame::Parse(Flame& value, IO::ParmParse& pp)
28 : {
29 : BL_PROFILE("Integrator::Flame::Flame()");
30 : {
31 3 : pp_query_default("timestep", value.base_time, 1.0e-4); //Simulation timestep
32 3 : pp_query_default("pf.eps", value.pf.eps, 0.0); // Burn width thickness
33 3 : pp_query_default("pf.kappa", value.pf.kappa, 0.0); // Interface energy param
34 3 : pp_query_default("pf.gamma", value.pf.gamma, 1.0); // Scaling factor for mobility
35 3 : pp_query_default("pf.lambda", value.pf.lambda, 0.0); // Chemical potential multiplier
36 3 : pp_query_default("pf.w1", value.pf.w1, 0.0); // Unburned rest energy
37 3 : pp_query_default("pf.w12", value.pf.w12, 0.0); // Barrier energy
38 3 : pp_query_default("pf.w0", value.pf.w0, 0.0); // Burned rest energy
39 3 : pp_query_default("amr.ghost_cells", value.ghost_count, 2); // number of ghost cells in all fields
40 3 : pp_query_default("geometry.x_len", value.x_len, 0.001); // Domain x length
41 3 : pp_query_default("geometry.y_len", value.y_len, 0.001); // Domain y length
42 :
43 : //value.bc_eta = new BC::Constant(1);
44 :
45 :
46 3 : pp.select<BC::Constant>("pf.eta.bc", value.bc_eta, 1 );
47 3 : value.RegisterNewFab(value.eta_mf, value.bc_eta, 1, value.ghost_count, "eta", true);
48 3 : value.RegisterNewFab(value.eta_old_mf, value.bc_eta, 1, value.ghost_count, "eta_old", false);
49 :
50 : // phase field initial condition
51 3 : pp.select<IC::Constant,IC::Expression>("pf.eta.ic",value.ic_eta,value.geom); // after comment
52 :
53 : // phase field initial condition
54 3 : pp.select<IC::Laminate,IC::Constant,IC::Expression,IC::BMP,IC::PNG>("eta.ic",value.ic_eta,value.geom);
55 : }
56 :
57 : {
58 : //IO::ParmParse pp("thermal");
59 3 : pp_query_default("thermal.on", value.thermal.on, false); // Whether to use the Thermal Transport Model
60 3 : pp_query_default("elastic.on", value.elastic.on, 0); // Whether to use Neo-hookean Elastic model
61 3 : pp_query_default("thermal.bound", value.thermal.bound, 0.0); // System Initial Temperature
62 3 : pp_query_default("elastic.traction", value.elastic.traction, 0.0); // Body force
63 3 : pp_query_default("elastic.phirefinement", value.elastic.phirefinement, 1); // Phi refinement criteria
64 :
65 :
66 :
67 3 : if (value.thermal.on) {
68 1 : pp_query_required("thermal.rho_ap", value.thermal.rho_ap); // AP Density
69 1 : pp_query_required("thermal.rho_htpb", value.thermal.rho_htpb); // HTPB Density
70 1 : pp_query_required("thermal.k_ap", value.thermal.k_ap); // AP Thermal Conductivity
71 1 : pp_query_required("thermal.k_htpb", value.thermal.k_htpb); // HTPB Thermal Conductivity
72 1 : pp_query_required("thermal.cp_ap", value.thermal.cp_ap); // AP Specific Heat
73 1 : pp_query_required("thermal.cp_htpb", value.thermal.cp_htpb); //HTPB Specific Heat
74 :
75 1 : pp_query_default("thermal.q0", value.thermal.q0, 0.0); // Baseline heat flux
76 :
77 1 : pp_query_required("thermal.m_ap", value.thermal.m_ap); // AP Pre-exponential factor for Arrhenius Law
78 1 : pp_query_required("thermal.m_htpb", value.thermal.m_htpb); // HTPB Pre-exponential factor for Arrhenius Law
79 1 : pp_query_required("thermal.E_ap", value.thermal.E_ap); // AP Activation Energy for Arrhenius Law
80 1 : pp_query_required("thermal.E_htpb", value.thermal.E_htpb); // HTPB Activation Energy for Arrhenius Law
81 :
82 1 : pp_query_default("thermal.hc", value.thermal.hc, 1.0); // Used to change heat flux units
83 1 : pp_query_default("thermal.massfraction", value.thermal.massfraction, 0.8); // Systen AP mass fraction
84 1 : pp_query_default("thermal.mlocal_ap", value.thermal.mlocal_ap, 0.0); // AP mass flux reference value
85 1 : pp_query_default("thermal.mlocal_htpb", value.thermal.mlocal_htpb, 0.0); // HTPB mass flux reference value
86 1 : pp_query_default("thermal.mlocal_comb", value.thermal.mlocal_comb, 0.0); // AP/HTPB mass flux reference value
87 :
88 1 : pp_query_default("thermal.T_fluid", value.thermal.T_fluid, 300.0); // Temperature of the Standin Fluid
89 :
90 1 : pp_query("thermal.disperssion1", value.thermal.disperssion1); // K; dispersion variables are use to set the outter field properties for the void grain case.
91 1 : pp_query("thermal.disperssion2", value.thermal.disperssion2); // rho; dispersion variables are use to set the outter field properties for the void grain case.
92 1 : pp_query("thermal.disperssion3", value.thermal.disperssion3); // cp; dispersion variables are use to set the outter field properties for the void grain case.
93 :
94 1 : pp_query_default("thermal.modeling_ap", value.thermal.modeling_ap, 1.0); // Scaling factor for AP thermal conductivity (default = 1.0)
95 1 : pp_query_default("thermal.modeling_htpb", value.thermal.modeling_htpb, 1.0); // Scaling factor for HTPB thermal conductivity (default = 1.0)
96 :
97 :
98 : //value.bc_temp = new BC::Constant(1);
99 1 : pp.select<BC::Constant>("thermal.temp.bc", value.bc_temp, 1);
100 :
101 1 : value.RegisterNewFab(value.temp_mf, value.bc_temp, 1, value.ghost_count + 1, "temp", true);
102 1 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, 1, value.ghost_count + 1, "temp_old", false);
103 1 : value.RegisterNewFab(value.temps_mf, value.bc_temp, 1, value.ghost_count + 1, "temps", false);
104 1 : value.RegisterNewFab(value.temps_old_mf, value.bc_temp, 1, value.ghost_count + 1, "temps_old", false);
105 :
106 1 : value.RegisterNewFab(value.mdot_mf, value.bc_eta, 1, value.ghost_count + 1, "mdot", value.plot_field);
107 1 : value.RegisterNewFab(value.mob_mf, value.bc_eta, 1, value.ghost_count + 1, "mob", value.plot_field);
108 1 : value.RegisterNewFab(value.alpha_mf, value.bc_temp, 1, value.ghost_count + 1, "alpha", value.plot_field);
109 1 : value.RegisterNewFab(value.heatflux_mf, value.bc_temp, 1, value.ghost_count + 1, "heatflux", value.plot_field);
110 1 : value.RegisterNewFab(value.laser_mf, value.bc_temp, 1, value.ghost_count + 1, "laser", value.plot_field);
111 :
112 1 : value.RegisterIntegratedVariable(&value.volume, "total_area");
113 1 : value.RegisterIntegratedVariable(&value.area, "Interface_area");
114 1 : value.RegisterIntegratedVariable(&value.chamber_area, "chamber_area", false);
115 1 : value.RegisterIntegratedVariable(&value.massflux, "mass_flux");
116 1 : value.RegisterIntegratedVariable(&value.chamber_pressure, "Pressure", false);
117 :
118 : // laser initial condition
119 1 : pp.select<IC::Constant,IC::Expression>("laser.ic",value.ic_laser, value.geom);
120 :
121 : // thermal initial condition
122 1 : pp.select<IC::Constant,IC::Expression,IC::BMP,IC::PNG>("temp.ic",value.thermal.ic_temp,value.geom);
123 : }
124 : }
125 :
126 : {
127 3 : pp_query_default("pressure.P", value.pressure.P, 1.0); // Constant pressure value
128 3 : if (value.thermal.on)
129 : {
130 1 : pp_query_required("pressure.a1", value.pressure.arrhenius.a1); // Surgate heat flux model paramater - AP
131 1 : pp_query_required("pressure.a2", value.pressure.arrhenius.a2); // Surgate heat flux model paramater - HTPB
132 1 : pp_query_required("pressure.a3", value.pressure.arrhenius.a3); // Surgate heat flux model paramater - Total
133 1 : pp_query_required("pressure.b1", value.pressure.arrhenius.b1); // Surgate heat flux model paramater - AP
134 1 : pp_query_required("pressure.b2", value.pressure.arrhenius.b2); // Surgate heat flux model paramater - HTPB
135 1 : pp_query_required("pressure.b3", value.pressure.arrhenius.b3); // Surgate heat flux model paramater - Total
136 1 : pp_query_required("pressure.c1", value.pressure.arrhenius.c1); // Surgate heat flux model paramater - Total
137 1 : pp_query_default("pressure.mob_ap", value.pressure.arrhenius.mob_ap, 0); // Whether to include pressure to the arrhenius law
138 1 : pp_query_default("pressure.dependency", value.pressure.arrhenius.dependency, 1); // Whether to use pressure to determined the reference Zeta
139 1 : pp_query_default("pressure.h1", value.pressure.arrhenius.h1, 1.81); // Surgate heat flux model paramater - Homogenized
140 1 : pp_query_default("pressure.h2", value.pressure.arrhenius.h2, 1.34); // Surgate heat flux model paramater - Homogenized
141 :
142 : }
143 : else
144 : {
145 2 : pp_query_required("pressure.r_ap", value.pressure.power.r_ap); // AP power pressure law parameter (r*P^n)
146 2 : pp_query_required("pressure.r_htpb", value.pressure.power.r_htpb); // HTPB power pressure law parameter (r*P^n)
147 2 : pp_query_required("pressure.r_comb", value.pressure.power.r_comb); // AP/HTPB power pressure law parameter (r*P^n)
148 2 : pp_query_required("pressure.n_ap", value.pressure.power.n_ap); // AP power pressure law parameter (r*P^n)
149 2 : pp_query_required("pressure.n_htpb", value.pressure.power.n_htpb); // HTPB power pressure law parameter (r*P^n)
150 2 : pp_query_required("pressure.n_comb", value.pressure.power.n_comb); // AP/HTPB power pressure law parameter (r*P^n)
151 : }
152 3 : pp_query_default("variable_pressure", value.variable_pressure, 0); // Whether to compute the pressure evolution
153 3 : pp_query_default("homogeneousSystem", value.homogeneousSystem, 0); // Whether to initialize Phi with homogenized properties
154 : }
155 :
156 3 : pp_query_default("amr.refinement_criterion", value.m_refinement_criterion, 0.001);// Refinement criterion for eta field
157 3 : pp_query_default("amr.refinement_criterion_temp", value.t_refinement_criterion, 0.001);// Refinement criterion for temperature field
158 3 : pp_query_default("amr.refinament_restriction", value.t_refinement_restriction, 0.1);// Eta value to restrict the refinament for the temperature field
159 3 : pp_query_default("amr.phi_refinement_criterion", value.phi_refinement_criterion, 1.0e100);// Refinement criterion for phi field [infinity]
160 3 : pp_query_default("small", value.small, 1.0e-8); // Lowest value of Eta.
161 :
162 : {
163 : // The material field is referred to as :math:`\phi(\mathbf{x})` and is
164 : // specified using these parameters.
165 3 : std::string phi_ic_type = "packedspheres";
166 3 : pp_query_validate("phi.ic.type", phi_ic_type, {"psread", "laminate", "expression", "constant", "bmp", "png"}); // IC type (psread, laminate, constant)
167 3 : if (phi_ic_type == "psread") {
168 0 : value.ic_phi = new IC::PSRead(value.geom, pp, "phi.ic.psread");
169 : //value.ic_phicell = new IC::PSRead(value.geom, pp, "phi.ic.psread");
170 0 : pp_query_default("phi.ic.psread.eps", value.zeta, 1.0e-5); // AP/HTPB interface length
171 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
172 : }
173 3 : else if (phi_ic_type == "laminate") {
174 3 : value.ic_phi = new IC::Laminate(value.geom, pp, "phi.ic.laminate");
175 : //value.ic_phicell = new IC::Laminate(value.geom, pp, "phi.ic.laminate");
176 3 : pp_query_default("phi.ic.laminate.eps", value.zeta, 1.0e-5); // AP/HTPB interface length
177 3 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
178 : }
179 0 : else if (phi_ic_type == "expression") {
180 0 : value.ic_phi = new IC::Expression(value.geom, pp, "phi.ic.expression");
181 : //value.ic_phicell = new IC::Expression(value.geom, pp, "phi.ic.expression");
182 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
183 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
184 : }
185 0 : else if (phi_ic_type == "constant") {
186 0 : value.ic_phi = new IC::Constant(value.geom, pp, "phi.ic.constant");
187 : //value.ic_phicell = new IC::Constant(value.geom, pp, "phi.ic.constant");
188 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
189 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
190 : }
191 0 : else if (phi_ic_type == "bmp") {
192 0 : value.ic_phi = new IC::BMP(value.geom, pp, "phi.ic.bmp");
193 : //value.ic_phicell = new IC::BMP(value.geom, pp, "phi.ic.bmp");
194 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
195 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
196 : }
197 0 : else if (phi_ic_type == "png") {
198 0 : value.ic_phi = new IC::PNG(value.geom, pp, "phi.ic.png");
199 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
200 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
201 : }
202 0 : else Util::Abort(INFO, "Invalid IC type ", phi_ic_type);
203 : //value.RegisterNewFab(value.phi_mf, 1, "phi_cell", true);
204 3 : value.RegisterNodalFab(value.phi_mf, 1, value.ghost_count + 1, "phi", true);
205 : //value.RegisterNewFab(value.phicell_mf, value.bc_eta, 1, value.ghost_count + 1, "phi", true);
206 : }
207 :
208 3 : pp.queryclass<Base::Mechanics<model_type>>("elastic",value);
209 :
210 3 : if (value.m_type != Type::Disable)
211 : {
212 0 : value.elastic.Tref = value.thermal.bound;
213 0 : pp_query_default("Tref", value.elastic.Tref, 300.0); // Initial temperature for thermal expansion computation
214 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_ap", value.elastic.model_ap);
215 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_htpb", value.elastic.model_htpb);
216 :
217 0 : value.bc_psi = new BC::Nothing();
218 0 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, value.ghost_count, "psi", value.plot_psi);
219 0 : value.psi_on = true;
220 : }
221 3 : }
222 :
223 24 : void Flame::Initialize(int lev)
224 : {
225 : BL_PROFILE("Integrator::Flame::Initialize");
226 24 : Util::Message(INFO, m_type);
227 24 : Base::Mechanics<model_type>::Initialize(lev);
228 :
229 24 : ic_eta->Initialize(lev, eta_mf);
230 24 : ic_eta->Initialize(lev, eta_old_mf);
231 24 : ic_phi->Initialize(lev, phi_mf);
232 : //ic_phicell->Initialize(lev, phicell_mf);
233 :
234 24 : if (elastic.on) {
235 0 : psi_mf[lev]->setVal(1.0);
236 0 : rhs_mf[lev]->setVal(Set::Vector::Zero());
237 : }
238 24 : if (thermal.on) {
239 8 : if (thermal.ic_temp)
240 : {
241 8 : thermal.ic_temp->Initialize(lev,temp_mf);
242 8 : thermal.ic_temp->Initialize(lev,temp_old_mf);
243 8 : thermal.ic_temp->Initialize(lev,temps_mf);
244 8 : thermal.ic_temp->Initialize(lev,temps_old_mf);
245 : }
246 : else
247 : {
248 0 : temp_mf[lev]->setVal(thermal.bound);
249 0 : temp_old_mf[lev]->setVal(thermal.bound);
250 0 : temps_mf[lev]->setVal(thermal.bound);
251 0 : temps_old_mf[lev]->setVal(thermal.bound);
252 : }
253 8 : alpha_mf[lev]->setVal(0.0);
254 8 : mob_mf[lev]->setVal(0.0);
255 8 : mdot_mf[lev]->setVal(0.0);
256 8 : heatflux_mf[lev]->setVal(0.0);
257 : // pressure_mf[lev]->setVal(1.0); // [error]
258 8 : thermal.w1 = 0.2 * pressure.P + 0.9;
259 8 : thermal.T_fluid = thermal.bound;
260 8 : ic_laser->Initialize(lev, laser_mf);
261 8 : thermal.mlocal_htpb = 685000.0 - 850e3 * thermal.massfraction;
262 : }
263 24 : if (variable_pressure) pressure.P = 1.0;
264 24 : }
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 : //psi_mf[lev]->setVal(1.0);
277 0 : phi_mf[lev]->FillBoundary();
278 : //phicell_mf[lev]->FillBoundary();
279 0 : eta_mf[lev]->FillBoundary();
280 0 : temp_mf[lev]->FillBoundary();
281 :
282 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
283 : {
284 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
285 : //amrex::Box bx = mfi.nodaltilebox();
286 : //bx.grow(1);
287 0 : amrex::Array4<model_type> const& model = model_mf[lev]->array(mfi);
288 0 : amrex::Array4<const Set::Scalar> const& phi = phi_mf[lev]->array(mfi);
289 0 : amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
290 0 : amrex::Array4<Set::Vector> const& rhs = rhs_mf[lev]->array(mfi);
291 : // amrex::Array4<const Set::Scalar> const& Pressure = pressure_mf[lev]->array(mfi); // [error]
292 :
293 0 : if (elastic.on)
294 : {
295 0 : amrex::Array4<const Set::Scalar> const& temp = temp_mf[lev]->array(mfi);
296 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
297 : {
298 : //auto sten = Numeric::GetStencil(i, j, k, bx);
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 : Set::Vector grad_eta = Numeric::CellGradientOnNode(eta, i, j, k, 0, DX);
302 0 : rhs(i, j, k) = elastic.traction * grad_eta;
303 :
304 0 : model_type model_ap = elastic.model_ap;
305 0 : model_ap.F0 -= Set::Matrix::Identity();
306 0 : model_ap.F0 *= (temp_avg - elastic.Tref);
307 0 : model_ap.F0 += Set::Matrix::Identity();
308 0 : model_type model_htpb = elastic.model_htpb;
309 0 : model_htpb.F0 -= Set::Matrix::Identity();
310 0 : model_htpb.F0 *= (temp_avg - elastic.Tref);
311 0 : model_htpb.F0 += Set::Matrix::Identity();
312 :
313 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
314 0 : });
315 : }
316 : else
317 : {
318 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
319 : {
320 0 : Set::Scalar phi_avg = Numeric::Interpolate::CellToNodeAverage(phi, i, j, k, 0);
321 : //phi_avg = phi(i,j,k,0);
322 0 : model_type model_ap = elastic.model_ap;
323 0 : model_ap.F0 *= Set::Matrix::Zero();
324 0 : model_type model_htpb = elastic.model_htpb;
325 0 : model_htpb.F0 *= Set::Matrix::Zero();
326 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
327 0 : });
328 : }
329 : }
330 0 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
331 :
332 0 : psi_mf[lev]->setVal(1.0);
333 0 : amrex::MultiFab::Copy(*psi_mf[lev], *eta_mf[lev], 0, 0, 1, psi_mf[lev]->nGrow());
334 : }
335 : }
336 :
337 212 : void Flame::TimeStepBegin(Set::Scalar a_time, int a_iter)
338 : {
339 : BL_PROFILE("Integrator::Flame::TimeStepBegin");
340 212 : Base::Mechanics<model_type>::TimeStepBegin(a_time, a_iter);
341 212 : if (thermal.on) {
342 90 : for (int lev = 0; lev <= finest_level; ++lev)
343 80 : ic_laser->Initialize(lev, laser_mf, a_time);
344 : }
345 212 : }
346 :
347 212 : void Flame::TimeStepComplete(Set::Scalar /*a_time*/, int /*a_iter*/)
348 : {
349 : BL_PROFILE("Integrator::Flame::TimeStepComplete");
350 212 : if (variable_pressure) {
351 0 : Set::Scalar domain_area = x_len * y_len;
352 0 : chamber_pressure = pressure.P;
353 0 : chamber_area = domain_area - volume;
354 0 : Util::Message(INFO, "Mass = ", massflux);
355 0 : Util::Message(INFO, "Pressure = ", pressure.P);
356 : }
357 212 : }
358 :
359 54060 : void Flame::Advance(int lev, Set::Scalar time, Set::Scalar dt)
360 : {
361 : BL_PROFILE("Integrador::Flame::Advance");
362 54060 : Base::Mechanics<model_type>::Advance(lev, time, dt);
363 54060 : const Set::Scalar* DX = geom[lev].CellSize();
364 :
365 : if (true) //lev == finest_level) //(true)
366 : {
367 54060 : if (thermal.on)
368 : {
369 2550 : std::swap(eta_old_mf[lev], eta_mf[lev]);
370 2550 : std::swap(temp_old_mf[lev], temp_mf[lev]);
371 2550 : std::swap(temps_old_mf[lev], temps_mf[lev]);
372 :
373 : Numeric::Function::Polynomial<4> w(pf.w0,
374 : 0.0,
375 2550 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
376 2550 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
377 2550 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0);
378 2550 : Numeric::Function::Polynomial<3> dw = w.D();
379 :
380 13100 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
381 : {
382 10550 : const amrex::Box& bx = mfi.tilebox();
383 : // Phase fields
384 10550 : amrex::Array4<Set::Scalar> const& etanew = (*eta_mf[lev]).array(mfi);
385 10550 : amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
386 10550 : amrex::Array4<const Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
387 : //amrex::Array4<const Set::Scalar> const& phicell = (*phicell_mf[lev]).array(mfi);
388 : // Heat transfer fields
389 10550 : amrex::Array4<Set::Scalar> const& tempnew = (*temp_mf[lev]).array(mfi);
390 10550 : amrex::Array4<const Set::Scalar> const& temp = (*temp_old_mf[lev]).array(mfi);
391 10550 : amrex::Array4<Set::Scalar> const& alpha = (*alpha_mf[lev]).array(mfi);
392 10550 : amrex::Array4<Set::Scalar> const& tempsnew = (*temps_mf[lev]).array(mfi);
393 10550 : amrex::Array4<Set::Scalar> const& temps = (*temps_old_mf[lev]).array(mfi);
394 10550 : amrex::Array4<Set::Scalar> const& laser = (*laser_mf[lev]).array(mfi);
395 : // Diagnostic fields
396 10550 : amrex::Array4<Set::Scalar> const& mob = (*mob_mf[lev]).array(mfi);
397 10550 : amrex::Array4<Set::Scalar> const& mdot = (*mdot_mf[lev]).array(mfi);
398 10550 : amrex::Array4<Set::Scalar> const& heatflux = (*heatflux_mf[lev]).array(mfi);
399 :
400 10550 : if (variable_pressure) {
401 0 : pressure.P = exp(0.00075 * massflux);
402 0 : if (pressure.P > 10.0) {
403 0 : pressure.P = 10.0;
404 : }
405 0 : else if (pressure.P <= 0.99) {
406 0 : pressure.P = 0.99;
407 : }
408 0 : elastic.traction = pressure.P;
409 :
410 : }
411 :
412 10550 : Set::Scalar zeta_2 = 0.000045 - pressure.P * 6.42e-6;
413 : Set::Scalar zeta_1;
414 10550 : if (pressure.arrhenius.dependency == 1) zeta_1 = zeta_2;
415 0 : else zeta_1 = zeta_0;
416 :
417 10550 : Set::Scalar k1 = pressure.arrhenius.a1 * pressure.P + pressure.arrhenius.b1 - zeta_1 / zeta;
418 10550 : Set::Scalar k2 = pressure.arrhenius.a2 * pressure.P + pressure.arrhenius.b2 - zeta_1 / zeta;
419 10550 : Set::Scalar k3 = 4.0 * log((pressure.arrhenius.c1 * pressure.P * pressure.P + pressure.arrhenius.a3 * pressure.P + pressure.arrhenius.b3) - k1 / 2.0 - k2 / 2.0);
420 10550 : Set::Scalar k4 = pressure.arrhenius.h1 * pressure.P + pressure.arrhenius.h2;
421 :
422 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
423 : {
424 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
425 1747840 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
426 : Set::Scalar K; // Calculate effective thermal conductivity
427 : Set::Scalar rho; // No special interface mixure rule is needed here.
428 : Set::Scalar cp;
429 1747840 : if (homogeneousSystem) {
430 0 : K = (thermal.modeling_ap * thermal.k_ap * thermal.massfraction + thermal.modeling_htpb * thermal.k_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion1 * (1. - phi_avg); // Calculate effective thermal conductivity
431 0 : rho = (thermal.rho_ap * thermal.massfraction + thermal.rho_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion2 * (1. - phi_avg); // No special interface mixure rule is needed here.
432 0 : cp = (thermal.cp_ap * thermal.massfraction + thermal.cp_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion3 * (1. - phi_avg);
433 : }
434 : else {
435 1747840 : K = thermal.modeling_ap * thermal.k_ap * phi_avg + thermal.modeling_htpb * thermal.k_htpb * (1.0 - phi_avg); // Calculate effective thermal conductivity
436 1747840 : rho = thermal.rho_ap * phi_avg + thermal.rho_htpb * (1.0 - phi_avg); // No special interface mixure rule is needed here.
437 1747840 : cp = thermal.cp_ap * phi_avg + thermal.cp_htpb * (1.0 - phi_avg);
438 : }
439 3495680 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
440 5243520 : etanew(i, j, k) = eta(i, j, k) - mob(i, j, k) * dt * df_deta;
441 3495680 : if (etanew(i, j, k) <= small) etanew(i, j, k) = small;
442 :
443 1747840 : alpha(i, j, k) = K / rho / cp; // Calculate thermal diffusivity and store in fiel
444 5243520 : mdot(i, j, k) = rho * fabs(eta(i, j, k) - etanew(i, j, k)) / dt; // deta/dt
445 :
446 6991360 : if (isnan(etanew(i, j, k)) || isnan(alpha(i, j, k)) || isnan(mdot(i, j, k))) {
447 0 : Util::Message(INFO, etanew(i, j, k), "etanew contains nan (i=", i, " j= ", j, ")");
448 0 : Util::Message(INFO, mdot(i, j, k), "mdot contains nan (i=", i, " j= ", j, ")");
449 0 : Util::Message(INFO, alpha(i, j, k), "alpha contains nan (i=", i, " j= ", j, ")");
450 0 : Util::Abort(INFO);
451 : }
452 1747840 : });
453 :
454 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
455 : {
456 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
457 : //phi_avg = phicell(i,j,k);
458 : Set::Scalar K;
459 : Set::Scalar qflux;
460 : Set::Scalar mlocal;
461 : Set::Scalar mdota;
462 : Set::Scalar mbase;
463 :
464 1747840 : if (homogeneousSystem) {
465 0 : qflux = k4 * phi_avg;
466 0 : mlocal = (thermal.mlocal_ap) * thermal.massfraction + (thermal.mlocal_htpb) * (1.0 - thermal.massfraction);
467 0 : mdota = fabs(mdot(i, j, k));
468 0 : mbase = tanh(4.0 * mdota / (mlocal));
469 0 : K = (thermal.modeling_ap * thermal.k_ap * thermal.massfraction + thermal.modeling_htpb * thermal.k_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion1 * (1. - phi_avg);
470 0 : heatflux(i, j, k) = (laser(i, j, k) * phi_avg + thermal.hc * mbase * qflux) / K;
471 : }
472 : else {
473 1747840 : qflux = k1 * phi_avg +
474 1747840 : k2 * (1.0 - phi_avg) +
475 1747840 : (zeta_1 / zeta) * exp(k3 * phi_avg * (1.0 - phi_avg));
476 1747840 : mlocal = (thermal.mlocal_ap) * phi_avg + (thermal.mlocal_htpb) * (1.0 - phi_avg) + thermal.mlocal_comb * phi_avg * (1.0 - phi_avg);
477 1747840 : mdota = fabs(mdot(i, j, k));
478 1747840 : mbase = tanh(4.0 * mdota / (mlocal));
479 :
480 1747840 : K = thermal.modeling_ap * thermal.k_ap * phi_avg + thermal.modeling_htpb * thermal.k_htpb * (1.0 - phi_avg);
481 5243520 : heatflux(i, j, k) = (thermal.hc * mbase * qflux + laser(i, j, k)) / K;
482 : }
483 :
484 3495680 : if (isnan(heatflux(i, j, k))) {
485 0 : Util::Message(INFO, heatflux(i, j, k), "heat contains nan (i=", i, " j= ", j, ")");
486 0 : Util::Abort(INFO);
487 : }
488 1747840 : });
489 :
490 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
491 : {
492 1747840 : auto sten = Numeric::GetStencil(i, j, k, bx);
493 1747840 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
494 1747840 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
495 1747840 : Set::Scalar lap_temp = Numeric::Laplacian(temp, i, j, k, 0, DX);
496 1747840 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
497 1747840 : Set::Vector grad_alpha = Numeric::Gradient(alpha, i, j, k, 0, DX, sten);
498 1747840 : Set::Scalar dTdt = 0.0;
499 3495680 : dTdt += grad_eta.dot(grad_temp * alpha(i, j, k));
500 3495680 : dTdt += grad_alpha.dot(eta(i, j, k) * grad_temp);
501 3495680 : dTdt += eta(i, j, k) * alpha(i, j, k) * lap_temp;
502 3495680 : dTdt += alpha(i, j, k) * heatflux(i, j, k) * grad_eta_mag;
503 : Set::Scalar Tsolid;
504 5243520 : Tsolid = dTdt + temps(i, j, k) * (etanew(i, j, k) - eta(i, j, k)) / dt;
505 3495680 : tempsnew(i, j, k) = temps(i, j, k) + dt * Tsolid;
506 6991360 : tempnew(i, j, k) = etanew(i, j, k) * tempsnew(i, j, k) + (1.0 - etanew(i, j, k)) * thermal.T_fluid;
507 5243520 : if (isnan(tempsnew(i, j, k)) || isnan(temps(i, j, k))) {
508 0 : Util::Message(INFO, tempsnew(i, j, k), "tempsnew contains nan (i=", i, " j= ", j, ")");
509 0 : Util::Message(INFO, temps(i, j, k), "temps contains nan (i=", i, " j= ", j, ")");
510 0 : Util::Abort(INFO);
511 : }
512 :
513 1747840 : });
514 :
515 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
516 : {
517 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
518 : //phi_avg = phicell(i,j,k);
519 :
520 : Set::Scalar L;
521 1747840 : if (pressure.arrhenius.mob_ap == 1) L = thermal.m_ap * pressure.P * exp(-thermal.E_ap / tempnew(i, j, k)) * phi_avg;
522 3495680 : else L = thermal.m_ap * exp(-thermal.E_ap / tempnew(i, j, k)) * phi_avg;
523 1747840 : L += thermal.m_htpb * exp(-thermal.E_htpb / tempnew(i, j, k)) * (1.0 - phi_avg);
524 : //L += thermal.m_comb * (0.5 * tempnew(i, j, k) / thermal.bound) * phi(i, j, k) * (1.0 - phi(i, j, k));
525 3501040 : if (tempnew(i, j, k) <= thermal.bound) mob(i, j, k) = 0;
526 3484960 : else mob(i, j, k) = L;
527 : //mob(i,j,k) = L;
528 3495680 : if (isnan(mob(i, j, k))) {
529 0 : Util::Message(INFO, mob(i, j, k), "mob contains nan (i=", i, " j= ", j, ")");
530 0 : Util::Abort(INFO);
531 : }
532 1747840 : });
533 : } // MFi For loop
534 :
535 : } // thermal IF
536 : else
537 : {
538 51510 : std::swap(eta_old_mf[lev], eta_mf[lev]);
539 : Numeric::Function::Polynomial<4> w(pf.w0,
540 : 0.0,
541 51510 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
542 51510 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
543 51510 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0
544 51510 : );
545 51510 : Numeric::Function::Polynomial<3> dw = w.D();
546 :
547 154876 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
548 : {
549 103366 : const amrex::Box& bx = mfi.tilebox();
550 : //Phase Fields
551 103366 : amrex::Array4<Set::Scalar> const& etanew = (*eta_mf[lev]).array(mfi);
552 103366 : amrex::Array4<Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
553 103366 : amrex::Array4<Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
554 :
555 103366 : Set::Scalar fmod_ap = pressure.power.r_ap * pow(pressure.P, pressure.power.n_ap);
556 103366 : Set::Scalar fmod_htpb = pressure.power.r_htpb * pow(pressure.P, pressure.power.n_htpb);
557 103366 : Set::Scalar fmod_comb = pressure.power.r_comb * pow(pressure.P, pressure.power.n_comb);
558 :
559 103366 : pressure.power.a_fit = -1.16582 * sin(pressure.P) - 0.681788 * cos(pressure.P) + 3.3563;
560 103366 : pressure.power.b_fit = -0.708225 * sin(pressure.P) + 0.548067 * cos(pressure.P) + 1.55985;
561 103366 : pressure.power.c_fit = -0.0130849 * sin(pressure.P) - 0.03597 * cos(pressure.P) + 0.00725694;
562 :
563 :
564 103366 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
565 : {
566 : //Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
567 : Set::Scalar L;
568 95704000 : if (homogeneousSystem == 1) {
569 0 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
570 0 : Set::Scalar angle = acos(grad_eta[0] / grad_eta.lpNorm<2>()) * 180 / 3.1415;
571 0 : if (angle > 90) angle = angle - 90.0;
572 0 : if (angle > 180) angle = angle - 180.0;
573 0 : if (angle > 270) angle = angle - 270.0;
574 0 : L = pressure.power.a_fit + pressure.power.b_fit * exp(-pressure.power.c_fit * angle);
575 : }
576 : else {
577 : Set::Scalar fs_actual;
578 95704000 : fs_actual = fmod_ap * phi(i, j, k)
579 95704000 : + fmod_htpb * (1.0 - phi(i, j, k))
580 191408000 : + 4.0 * fmod_comb * phi(i, j, k) * (1.0 - phi(i, j, k));
581 95704000 : L = fs_actual / pf.gamma / (pf.w1 - pf.w0);
582 : }
583 95704000 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
584 191408000 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
585 191408000 : etanew(i, j, k) = eta(i, j, k) - L * dt * df_deta;
586 :
587 311333000 : if (etanew(i, j, k) > eta(i, j, k)) etanew(i, j, k) = eta(i, j, k);
588 95704000 : });
589 : } // MFI For Loop
590 : } // thermal ELSE
591 : }// First IF
592 54060 : } //Function
593 :
594 :
595 263 : void Flame::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
596 : {
597 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
598 263 : Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
599 :
600 263 : const Set::Scalar* DX = geom[lev].CellSize();
601 263 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
602 :
603 : // Eta criterion for refinement
604 709 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
605 : {
606 446 : const amrex::Box& bx = mfi.tilebox();
607 446 : amrex::Array4<char> const& tags = a_tags.array(mfi);
608 446 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
609 :
610 446 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
611 : {
612 78848 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
613 103214 : if (gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
614 47748 : tags(i, j, k) = amrex::TagBox::SET;
615 78848 : });
616 : }
617 :
618 : // Phi criterion for refinement
619 263 : if (elastic.phirefinement) {
620 709 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
621 : {
622 446 : const amrex::Box& bx = mfi.tilebox();
623 446 : amrex::Array4<char> const& tags = a_tags.array(mfi);
624 446 : amrex::Array4<const Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
625 :
626 446 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
627 : {
628 78848 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
629 78848 : if (gradphi.lpNorm<2>() * dr >= phi_refinement_criterion)
630 0 : tags(i, j, k) = amrex::TagBox::SET;
631 78848 : });
632 : }
633 : }
634 :
635 :
636 : // Thermal criterion for refinement
637 263 : if (thermal.on) {
638 507 : for (amrex::MFIter mfi(*temp_mf[lev], true); mfi.isValid(); ++mfi)
639 : {
640 345 : const amrex::Box& bx = mfi.tilebox();
641 345 : amrex::Array4<char> const& tags = a_tags.array(mfi);
642 345 : amrex::Array4<const Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
643 345 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
644 345 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
645 : {
646 50928 : Set::Vector tempgrad = Numeric::Gradient(temp, i, j, k, 0, DX);
647 50928 : if (tempgrad.lpNorm<2>() * dr > t_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
648 0 : tags(i, j, k) = amrex::TagBox::SET;
649 50928 : });
650 : }
651 : }
652 :
653 :
654 263 : }
655 :
656 15 : void Flame::Regrid(int lev, Set::Scalar time)
657 : {
658 : BL_PROFILE("Integrator::Flame::Regrid");
659 : //if (lev < finest_level) return;
660 : //phi_mf[lev]->setVal(0.0);
661 15 : ic_phi->Initialize(lev, phi_mf, time);
662 : //ic_phicell->Initialize(lev, phi_mf, time);
663 15 : }
664 :
665 220 : void Flame::Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
666 : const amrex::MFIter& mfi, const amrex::Box& box)
667 : {
668 : BL_PROFILE("Flame::Integrate");
669 220 : const Set::Scalar* DX = geom[amrlev].CellSize();
670 220 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
671 220 : amrex::Array4<amrex::Real> const& eta = (*eta_mf[amrlev]).array(mfi);
672 220 : amrex::Array4<amrex::Real> const& mdot = (*mdot_mf[amrlev]).array(mfi);
673 220 : if (variable_pressure) {
674 0 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
675 : {
676 0 : volume += eta(i, j, k, 0) * dv;
677 0 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
678 0 : Set::Scalar normgrad = grad.lpNorm<2>();
679 0 : Set::Scalar da = normgrad * dv;
680 0 : area += da;
681 :
682 0 : Set::Vector mgrad = Numeric::Gradient(mdot, i, j, k, 0, DX);
683 0 : Set::Scalar mnormgrad = mgrad.lpNorm<2>();
684 0 : Set::Scalar dm = mnormgrad * dv;
685 0 : massflux += dm;
686 :
687 0 : });
688 : }
689 : else {
690 220 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
691 : {
692 15560 : volume += eta(i, j, k, 0) * dv;
693 15560 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
694 15560 : Set::Scalar normgrad = grad.lpNorm<2>();
695 15560 : Set::Scalar da = normgrad * dv;
696 15560 : area += da;
697 15560 : });
698 : }
699 : // time dependent pressure data from experimenta -> p = 0.0954521220950523 * exp(15.289993148880678 * t)
700 220 : }
701 : } // namespace Integrator
702 :
703 :
|
