Line data Source code
1 : //
2 : // This class implements second and fourth order phase field brittle fracture with near singular solver.
3 : // This class inherits from the base class :code:`Base/Mechanics.H`.
4 : // See `this link <https://doi.org/10.1007/s00466-023-02325-8>`_
5 : // for more information on this implementation.
6 : //
7 : // The energy functional for a second order model is given by
8 : // .. math::
9 : //
10 : // \mathcal{L}= \int_\Omega \left[\left(g(c) + \eta \right)W_0^+ + W_0^-\right] dV + \int_\Omega G_c \left[ \frac{w(c)}{4\xi} + \xi |\nabla c|^2\right] dV
11 : //
12 : // where :math:`g(c)` and :math:`w(c)` are interpolation functions specified by the user.
13 : // The fracture energy :math:`G_c` and crack length scale :math:`xi` are also read from input file.
14 : // The tension compression asymmetry is accounted by splitting the strain energy :math:`W_0` as
15 : // .. math::
16 : //
17 : // W_0^\pm = \frac{1}{2} \lambda (\operatorname{tr}\bm{\varepsilon}_\pm)^2 + \mu \operatorname{tr}\left(\bm{\varepsilon}_\pm^2\right), \quad \bm{\varepsilon_\pm} = \sum_{i=1}^d \left(\varepsilon_i \right)_\pm \hat{\bm{v}}_i\otimes\hat{\bm{v}}_i
18 : //
19 : // where :math:`\lambda` and :math:`\mu` are material models for linear elastic isotropic material, :math:`\bm{\varepsilon}_\pm` are the positive and negative components of the strain tensor :math:`\bm{\varepsilon}` computed through eigenvalue decomposition.
20 : // The fourth order model adds a laplacian term to the free energy functional.
21 : // The code performs a staggered solve where it solves the elastic problem implicitly and crack problem explicitly.
22 : //
23 : // Class methods:
24 : //
25 : // #. :code:`Fracture()`:
26 : // Basic constructor. Does nothing, and leaves all values initiated as NAN.
27 : // #. :code:`Fracture(IO::ParmParse &pp)`:
28 : // Calls the parser.
29 : // #. :code:`static void Parse(Fracture &value, IO::ParmParse &pp)`
30 : // Parses input file, ICs, BCs, and sets up multifabs appropriately
31 : // #. :code:`void Initialize(int lev) override`
32 : // Calls IC and sets up the initial crack geometry.
33 : // #. :code:`virtual void UpdateModel(int a_step) override`
34 : // Performs degradation by updating the :math:`\psi` field using :math:`g(c)`.
35 : // #. :code:`void TimeStepBegin(Set::Scalar a_time, int a_step) override`
36 : // Solves the elastic problem, performs eigen value decomposition of strain, and computes the positive part of the strain energy.
37 : // #. :code:`void Advance(int a_lev, amrex::Real a_time, amrex::Real a_dt)`
38 : // Advances the crack field by computing the variational derivative of energy functional with :math:`c`.
39 : // #. :code:`void TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, Set::Scalar a_time, int a_ngrow) override`
40 : // Refines the grid based on the crack field.
41 : // #. :code:`void Integrate(int amrlev, Set::Scalar time, int step, const amrex::MFIter &mfi, const amrex::Box &a_box) override`
42 : // Performs spatial integration of the driving force to check for convergence of crack problem
43 : // #. :code:`void TimeStepComplete(Set::Scalar /*time*/, int /* iter*/)`
44 : // Checks whether the solver should work on crack problem or elastic problem.
45 : //
46 : #ifndef INTEGRATOR_FRACTURE_H
47 : #define INTEGRATOR_FRACTURE_H
48 :
49 : #include "Integrator/Base/Mechanics.H"
50 :
51 : #include "Model/Interface/Crack/PFCZM.H"
52 : #include "Model/Solid/Linear/Isotropic.H"
53 :
54 : #include "BC/Constant.H"
55 : #include "IC/BMP.H"
56 : #include "IC/Ellipse.H"
57 : // #include "IC/Ellipsoid.H"
58 : #include "IC/Expression.H"
59 : #include "IC/IC.H"
60 : #include "IC/Laminate.H"
61 : #include "IC/Notch.H"
62 : #include "IC/PNG.H"
63 : #include "IC/PerturbedInterface.H"
64 :
65 : #include "Numeric/Stencil.H"
66 :
67 : #include <cmath>
68 : #include <eigen3/Eigen/Dense>
69 :
70 : namespace Integrator
71 : {
72 : using brittle_model = Model::Solid::Linear::Isotropic;
73 : using pfczm_crack_type = Model::Interface::Crack::PFCZM;
74 :
75 : class Fracture : virtual public Base::Mechanics<brittle_model>
76 : {
77 : public:
78 : static constexpr const char *name = "fracture";
79 : Fracture() : Base::Mechanics<brittle_model>() {}
80 :
81 0 : Fracture(IO::ParmParse &pp) : Base::Mechanics<brittle_model>()
82 : {
83 0 : Parse(*this, pp);
84 0 : }
85 :
86 : static void
87 0 : Parse(Fracture &value, IO::ParmParse &pp)
88 : {
89 0 : Base::Mechanics<brittle_model>::Parse(value, pp);
90 :
91 : // Material field related parsing
92 0 : pp.query("material.refinement_threshold", value.material.m_eta_ref_threshold);
93 :
94 : // Driving force threshold
95 0 : pp.query_default("driving_force_refinement_threshold", value.crack.driving_force_refinement_threshold, 1E100);
96 :
97 : // Let's figure out if we are doing multi-material simulations.
98 0 : std::string mat_ic_type;
99 0 : if (pp.contains("material.ic.type"))
100 : {
101 : // IC determining material distribution
102 0 : pp.select<IC::Laminate, IC::Ellipse, IC::PerturbedInterface, IC::BMP, IC::PNG, IC::Expression>("material.ic", value.material.ic, value.geom);
103 0 : value.material.is_ic = true;
104 : }
105 : else
106 : {
107 0 : value.material.is_ic = false;
108 : }
109 :
110 : // Let's query material properties for different materials in the system
111 0 : pp.queryclass_enumerate<brittle_model>("material.model", value.material.models);
112 0 : value.material.num_mat = value.material.models.size();
113 0 : Util::Assert(INFO, TEST(value.material.models.size() > 0));
114 :
115 : // This is the fab that stores the material field.
116 : // For a simple, single material simulation, this will be a uniform field with value 1 everywhere
117 0 : value.RegisterNodalFab(value.material.eta_mf, value.material.num_mat, 2, "eta", true);
118 :
119 : // Crack related parsing
120 0 : pp.query_default("crack.refinement_threshold", value.crack.refinement_threshold, 0.0001); // mesh refinement criteria
121 0 : pp.query_default("crack.df.beta", value.crack.beta, 0.0); // constant multiplier for fourth orther model with bilaplacian
122 0 : pp.query_default("crack.df.tol_rel", value.crack.tol_rel, 1E-3); // relative tolerance for convergence of driving force
123 0 : pp.query_default("crack.df.tol_abs", value.crack.tol_abs, 1E-3); // absolute tolerance for convergence of driving force
124 0 : pp.query_default("crack.df.max_iter", value.crack.max_iter, 500.); // Maximum number of solver iterations
125 0 : pp.query_default("crack.df.mult_Gc", value.crack.mult_Gc, 1.0); // constant multiplier for controlling fracture energy of interface
126 0 : pp.query_default("crack.df.mult_lap", value.crack.mult_lap, 1.0); // constant multiplier for second order model with laplacian
127 0 : pp.query_default("crack.df.el_mult", value.crack.el_mult, 1.0); // unit conversion multiplier between elastic problem and crack problem
128 :
129 : // Let's figure out if there is an initial crack or void
130 0 : pp.select_enumerate<IC::Constant,IC::Notch, IC::Ellipse, IC::Expression>("crack.ic", value.crack.ic, value.geom);
131 :
132 0 : value.RegisterNodalFab(value.crack.c_mf, 1, 2, "crack", true);
133 0 : value.RegisterNodalFab(value.crack.c_old_mf, 1, 2, "crack_old", true);
134 0 : value.RegisterNodalFab(value.crack.driving_force_mf, 5, 2, "driving_force", true);
135 0 : value.RegisterNodalFab(value.crack.energy_pristine_mf, 3, 2, "energy_pristine", true);
136 0 : value.RegisterNodalFab(value.crack.history_var_mf, 3, 2, "history_variable", true);
137 0 : value.RegisterIntegratedVariable(&(value.crack.driving_force_norm), "driving_force_norm");
138 0 : value.RegisterIntegratedVariable(&(value.crack.crack_l2_err), "crack_l2_err");
139 0 : value.RegisterIntegratedVariable(&(value.crack.crack_norm), "crack_norm");
140 :
141 : // By default the psi variable will be on for fracture simulations
142 0 : value.psi_on = true;
143 0 : value.bc_psi = new BC::Constant(1, pp, "crack.bc"); // boundary conditions for crack
144 0 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, 2, "psi", true);
145 :
146 : // This is needed for the specific crack model we are implementing
147 0 : pp.queryclass_enumerate<pfczm_crack_type>("crack.model",value.crack.cracktype);
148 0 : Util::Assert(INFO, TEST(value.material.models.size() > 0));
149 0 : }
150 :
151 : void
152 0 : Initialize(int lev) override
153 : {
154 0 : Base::Mechanics<brittle_model>::Initialize(lev);
155 :
156 0 : crack.c_mf[lev]->setVal(1.0);
157 0 : crack.c_old_mf[lev]->setVal(1.0);
158 :
159 : // For now, we are setting psi field to 1. This prevents nans.
160 : // We will modify the psi field later.
161 0 : psi_mf[lev]->setVal(1.0);
162 :
163 0 : if (crack.is_ic)
164 : {
165 0 : for (unsigned int i = 0; i < crack.ic.size(); i++)
166 : {
167 : // This initializes the crack field.
168 0 : crack.ic[i]->Add(lev, crack.c_mf);
169 0 : crack.ic[i]->Add(lev, crack.c_old_mf);
170 : }
171 : }
172 :
173 : // This initializes the material field
174 0 : material.eta_mf[lev]->setVal(1.0);
175 0 : if (material.is_ic)
176 0 : material.ic->Initialize(lev, material.eta_mf);
177 : else
178 0 : material.eta_mf[lev]->setVal(1.0);
179 :
180 : // set history variable to zero
181 0 : crack.energy_pristine_mf[lev]->setVal(0.0);
182 0 : crack.history_var_mf[lev]->setVal(0.0);
183 0 : }
184 :
185 : virtual void
186 0 : UpdateModel(int a_step, Set::Scalar /*a_time*/) override
187 : {
188 0 : if (m_type == Base::Mechanics<brittle_model>::Type::Disable)
189 0 : return;
190 :
191 : // This sets the model field for the very first time step
192 : // We should not be resetting model field everytime.
193 0 : if (a_step == 0)
194 : {
195 0 : for (int lev = 0; lev <= finest_level; ++lev)
196 : {
197 0 : material.eta_mf[lev]->FillBoundary();
198 :
199 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
200 : {
201 0 : amrex::Box bx = mfi.grownnodaltilebox();
202 0 : amrex::Array4<brittle_model> const &model = model_mf[lev]->array(mfi);
203 0 : amrex::Array4<const Set::Scalar> const &eta = material.eta_mf[lev]->array(mfi);
204 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
205 0 : model(i, j, k) = brittle_model::Zero();
206 0 : for (int n = 0; n < material.num_mat; n++)
207 0 : model(i, j, k) += eta(i, j, k, n) * material.models[n];
208 0 : });
209 0 : }
210 0 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
211 : }
212 : }
213 :
214 : // This is where we perform "degradation"
215 : // Essentially we will set the psi field based on the c^2 value.
216 0 : for (int lev = 0; lev <= finest_level; ++lev)
217 : {
218 0 : crack.c_mf[lev]->FillBoundary();
219 0 : for (MFIter mfi(*psi_mf[lev], false); mfi.isValid(); ++mfi)
220 : {
221 0 : amrex::Box bx = mfi.growntilebox();
222 0 : amrex::Array4<Set::Scalar> const &psi = psi_mf[lev]->array(mfi);
223 0 : amrex::Array4<const Set::Scalar> const &c = crack.c_mf[lev]->array(mfi);
224 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
225 : // psi(i,j,k,0) = crack.cracktype[0].g_phi(c(i,j,k,0));
226 0 : psi(i, j, k, 0) = crack.cracktype[0].g_phi(Numeric::Interpolate::NodeToCellAverage(c, i, j, k, 0));
227 0 : });
228 0 : }
229 0 : psi_mf[lev]->FillBoundary();
230 : }
231 : }
232 :
233 : void
234 0 : TimeStepBegin(Set::Scalar /*a_time*/, int a_step) override
235 : {
236 : // Deciding whether to do an elastic solve or not
237 0 : Util::Message(INFO, crack.driving_force_norm, " ", crack.driving_force_reference, " ", crack.driving_force_reference_prev, " ", crack.tol_rel);
238 0 : if ((crack.driving_force_norm / crack.driving_force_reference < crack.tol_rel) || crack.crack_prop_iter > crack.max_iter)
239 : {
240 0 : crack.crack_prop_iter = 0;
241 0 : elastic_do_solve_now = true;
242 : }
243 0 : if ((crack.driving_force_norm < crack.tol_abs) || crack.crack_prop_iter > crack.max_iter)
244 : {
245 0 : elastic_do_solve_now = true;
246 0 : crack.crack_prop_iter = 0;
247 : }
248 0 : if (increase_load_step)
249 : {
250 0 : Util::Message(INFO, "Load step = ", loadstep + 1);
251 0 : elastic_do_solve_now = true;
252 0 : crack.crack_prop_iter = 0;
253 0 : increase_load_step = false;
254 0 : set_new_reference = true;
255 : // if (loadstep == 0) set_initial_df_once = true;
256 0 : loadstep++;
257 : }
258 :
259 : // This means that the crack field has not converged.
260 0 : if (!elastic_do_solve_now)
261 : {
262 0 : crack.crack_prop_iter++;
263 0 : return;
264 : }
265 :
266 : // Doing an elastic solve
267 0 : Base::Mechanics<brittle_model>::TimeStepBegin((double)loadstep, a_step);
268 :
269 : // Computing pristine energy based on eigen value decomposition
270 0 : for (int lev = 0; lev <= disp_mf.finest_level; lev++)
271 : {
272 0 : amrex::Box domain = geom[lev].Domain();
273 0 : domain.convert(amrex::IntVect::TheNodeVector());
274 0 : Set::Vector DX(geom[lev].CellSize());
275 :
276 0 : for (MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
277 : {
278 0 : amrex::Box bx = mfi.grownnodaltilebox(); // & domain;
279 0 : amrex::Array4<Set::Scalar> const &eta = material.eta_mf[lev]->array(mfi);
280 0 : amrex::Array4<Set::Scalar> const &c = crack.c_mf[lev]->array(mfi);
281 0 : amrex::Array4<Set::Matrix> const &stress = stress_mf[lev]->array(mfi);
282 0 : amrex::Array4<Set::Matrix> const &strain = strain_mf[lev]->array(mfi);
283 0 : amrex::Array4<Set::Scalar> const &energy = crack.energy_pristine_mf[lev]->array(mfi);
284 0 : amrex::Array4<Set::Scalar> const &history_var = crack.history_var_mf[lev]->array(mfi);
285 :
286 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
287 : //==========================================================
288 : // The formulation below uses decomposition of stress
289 : // Set::Matrix sig = stress(i,j,k);
290 : // Eigen::SelfAdjointEigenSolver<Set::Matrix> eigensolver(sig);
291 : // Set::Vector eValues = eigensolver.eigenvalues();
292 : // Set::Matrix eVectors = eigensolver.eigenvectors();
293 :
294 : // Set::Matrix sig_p = Set::Matrix::Zero();
295 : // Set::Matrix sig_n = Set::Matrix::Zero();
296 :
297 : // for (int n = 0; n < AMREX_SPACEDIM; n++)
298 : // {
299 : // if (eValues(n) > 0.0) sig_p += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
300 : // else sig_n += eValues(n)*(eVectors.col(n)*eVectors.col(n).transpose());
301 : // }
302 :
303 : // Set::Matrix sig_n_dev = sig_n - (1.0/3.0)*sig_n.trace()*Set::Matrix::Identity();
304 :
305 : // for (int n = 0; n < material.num_mat; n++)
306 : // {
307 : // // computing energy based on tensile stress
308 : // energy(i,j,k,0) += eta(i,j,k,n)*material.models[n].W2(sig_p);
309 :
310 : // // compute energy from deviatoric portion
311 : // energy(i,j,k,1) += eta(i,j,k,n) * material.models[n].W2(sig_n_dev);
312 : // }
313 :
314 : // // Only update energy if it is increasing. (H^+ according to Miehe)
315 : // if (energy(i,j,k,0) < energy_old(i,j,k,0)) energy(i,j,k,0) = energy_old(i,j,k,0);
316 : // if (energy(i,j,k,1) < energy_old(i,j,k,1)) energy(i,j,k,1) = energy_old(i,j,k,1);
317 :
318 : //==========================================================
319 : // The formulation below uses Ambati's decomposition of strain
320 : // Perform eigenvalue decomposition of strain
321 0 : Set::Matrix eps = strain(i, j, k);
322 0 : Set::Matrix sig = stress(i, j, k);
323 0 : energy(i, j, k, 0) = 0.0;
324 0 : energy(i, j, k, 1) = 0.0;
325 0 : energy(i, j, k, 2) = 0.0;
326 :
327 0 : if (!crack.cracktype[0].mixed_mode())
328 : {
329 : // perform eigenvalue decomposition of strain tensor.
330 0 : Eigen::SelfAdjointEigenSolver<Set::Matrix> eigensolver(eps);
331 0 : Set::Vector eValues = eigensolver.eigenvalues();
332 0 : Set::Matrix eVectors = eigensolver.eigenvectors();
333 :
334 : // Reconstruct positive and negative counterparts of strain
335 0 : Set::Matrix eps_p = Set::Matrix::Zero();
336 0 : Set::Matrix eps_n = Set::Matrix::Zero();
337 :
338 0 : for (int n = 0; n < AMREX_SPACEDIM; n++)
339 : {
340 0 : if (eValues(n) > 0.0)
341 0 : eps_p += eValues(n) * (eVectors.col(n) * eVectors.col(n).transpose());
342 : else
343 0 : eps_n += eValues(n) * (eVectors.col(n) * eVectors.col(n).transpose());
344 : }
345 :
346 0 : for (int n = 0; n < material.num_mat; n++)
347 0 : energy(i, j, k, 0) += eta(i, j, k, n) * material.models[n].W(eps_p);
348 :
349 0 : if (energy(i, j, k, 0) > history_var(i, j, k, 0))
350 0 : history_var(i, j, k, 0) = energy(i, j, k, 0);
351 : }
352 : else
353 : {
354 : // Here we implement the 2023 model by Wang et al (DOI: 10.1016/j.apm.2022.12.006)
355 0 : Set::Vector crack_grad = Numeric::Gradient(c, i, j, k, 0, DX.data());
356 0 : if (crack_grad.lpNorm<2>() > 1.e-4)
357 : {
358 :
359 : // perform eigenvalue decomposition of stress tensor.
360 0 : Eigen::SelfAdjointEigenSolver<Set::Matrix> eigensolver(sig);
361 0 : Set::Vector eValues = eigensolver.eigenvalues();
362 0 : Set::Matrix eVectors = eigensolver.eigenvectors();
363 :
364 : // Let's order the evalues
365 0 : Set::Scalar sig1 = 0., sig2 = 0.;
366 0 : Set::Vector vec1 = Set::Vector::Zero(), vec2 = Set::Vector::Zero();
367 0 : if (eValues(0) > eValues(1))
368 : {
369 0 : sig1 = eValues(0);
370 0 : sig2 = eValues(1);
371 0 : vec1 = eVectors.col(0);
372 0 : vec2 = eVectors.col(1);
373 : }
374 : else
375 : {
376 0 : sig1 = eValues(1);
377 0 : sig2 = eValues(0);
378 0 : vec1 = eVectors.col(1);
379 0 : vec2 = eVectors.col(0);
380 : }
381 :
382 0 : Set::Scalar irwing_length = crack.cracktype[0].l_w();
383 0 : Set::Scalar c_alpha = crack.cracktype[0].c_alpha();
384 :
385 0 : Set::Scalar mohr_center = 0.5 * (sig1 + sig2);
386 0 : Set::Scalar mohr_radius = 0.5 * std::abs(sig1 - sig2);
387 :
388 0 : Set::Scalar f_theta = 0.0;
389 :
390 0 : if (crack.cracktype[0].failure_surface() == Model::Interface::Crack::PFCZM::FSType::WANG2023)
391 : {
392 : // Storing useful quantities for later.
393 0 : Set::Scalar sig_t = crack.cracktype[0].sig_t();
394 0 : Set::Scalar tau_f = crack.cracktype[0].tau_s();
395 0 : Set::Scalar chi = crack.cracktype[0].chi();
396 0 : Set::Scalar beta_bar = crack.cracktype[0].beta_bar();
397 :
398 0 : Set::Scalar beta = 1.0;
399 :
400 : // Now let's compute a candidate theta
401 : // Ideally, theta = 0 is always a valid solution and corresponds to sig_nn = sig1, tau_nm = 0
402 : // Set::Scalar theta = 0.0; // removed because: set but not used
403 0 : beta = sig1 < 0 ? beta_bar : 1.0;
404 0 : f_theta = (beta * (crack.el_mult * crack.el_mult * sig1 * sig1 / (sig_t * sig_t))) - 1;
405 :
406 : // Now we check for other feasible solutions. First we check for beta = 1.
407 0 : Set::Scalar sig_nn1 = 0., tau_nm1 = 0.;
408 0 : Set::Scalar theta_candidate1 = 0.0, f_candidate1 = 0.;
409 0 : Set::Scalar temp = std::abs((mohr_center / mohr_radius) * (chi * chi / (1.0 - (chi * chi))));
410 0 : if (temp >= 0.0 && temp <= 1.0)
411 : {
412 0 : theta_candidate1 = 0.5 * std::acos(temp);
413 :
414 : // Now check if this actually leas to a positive sig_nn value
415 0 : sig_nn1 = mohr_center + (mohr_radius * std::cos(2.0 * theta_candidate1));
416 0 : tau_nm1 = mohr_radius * std::sin(2.0 * theta_candidate1);
417 :
418 0 : if (sig_nn1 >= 0.0)
419 : {
420 0 : f_candidate1 = (crack.el_mult * crack.el_mult * sig_nn1 * sig_nn1 / (sig_t * sig_t)) + (crack.el_mult * crack.el_mult * tau_nm1 * tau_nm1 / (tau_f * tau_f)) - 1.0;
421 :
422 : // We only use this solution if the value of the failure surface is larger.
423 0 : if (f_candidate1 > f_theta)
424 : {
425 0 : f_theta = f_candidate1;
426 0 : beta = 1.0;
427 : }
428 : }
429 : }
430 :
431 : // Now we check for the third feasible solution for beta = beta_bar.
432 0 : Set::Scalar sig_nn2 = 0., tau_nm2 = 0.;
433 0 : Set::Scalar theta_candidate2 = 0.0, f_candidate2 = 0.;
434 0 : temp = std::abs((mohr_center / mohr_radius) * (beta_bar * chi * chi / (1.0 - (beta_bar * chi * chi))));
435 0 : if (temp >= 0.0 && temp <= 1.0)
436 : {
437 0 : theta_candidate2 = 0.5 * std::acos(temp);
438 :
439 : // Now check if this actually leas to a negative sig_nn value
440 0 : sig_nn2 = mohr_center + (mohr_radius * std::cos(2.0 * theta_candidate2));
441 0 : tau_nm2 = mohr_radius * std::sin(2.0 * theta_candidate2);
442 :
443 0 : if (sig_nn2 <= 0.0)
444 : {
445 0 : f_candidate2 = (beta_bar * crack.el_mult * crack.el_mult * sig_nn2 * sig_nn2 / (sig_t * sig_t)) + (crack.el_mult * crack.el_mult * tau_nm2 * tau_nm2 / (tau_f * tau_f)) - 1.0;
446 :
447 : // We only use this solution if the value of the failure surface is larger.
448 0 : if (f_candidate2 > f_theta)
449 : {
450 0 : f_theta = f_candidate2;
451 0 : beta = beta_bar;
452 : }
453 : }
454 : }
455 : }
456 0 : else if (crack.cracktype[0].failure_surface() == Model::Interface::Crack::PFCZM::FSType::MC)
457 : {
458 0 : Set::Scalar cohesion = crack.cracktype[0].cohesion();
459 0 : Set::Scalar friction = crack.cracktype[0].friction();
460 : // In this case, the angle theta is straightforward.
461 : // theta = (pi / 4) - (\phi / 2)
462 0 : Set::Scalar sig_nn = 0., tau_nm = 0.;
463 0 : Set::Scalar theta = std::atan(1) - std::atan(friction);
464 :
465 0 : sig_nn = crack.el_mult * (mohr_center + (mohr_radius * std::cos(2.0 * theta)));
466 0 : tau_nm = crack.el_mult * (mohr_radius * std::sin(2.0 * theta));
467 :
468 0 : f_theta = (tau_nm + (friction * sig_nn)) * (tau_nm + (friction * sig_nn)) / (cohesion * cohesion);
469 0 : f_theta = f_theta - 1.0;
470 :
471 : // Util::Message(INFO, "sig1 = ", sig1, ", sig2 = ", sig2, ", f_theta = ", f_theta);
472 : }
473 :
474 0 : energy(i,j,k,0) = (f_theta > 0) ? (c_alpha / (2.0 * irwing_length)) * ((f_theta)) : 0.0;
475 0 : energy(i,j,k,1) = 0.0;
476 :
477 0 : if (energy(i, j, k, 0) + energy(i, j, k, 1) > history_var(i, j, k, 0) + history_var(i, j, k, 1))
478 : {
479 0 : history_var(i, j, k, 0) = energy(i, j, k, 0);
480 0 : history_var(i, j, k, 1) = energy(i, j, k, 1);
481 : }
482 : }
483 : }
484 : //=====================================================================
485 0 : });
486 0 : }
487 0 : Util::RealFillBoundary(*crack.energy_pristine_mf[lev], geom[lev]);
488 0 : Util::RealFillBoundary(*crack.history_var_mf[lev], geom[lev]);
489 : }
490 :
491 0 : integrate_variables_before_advance = false;
492 0 : integrate_variables_after_advance = true;
493 : }
494 :
495 : void
496 0 : Advance(int a_lev, amrex::Real a_time, amrex::Real a_dt) override
497 : {
498 : // advance for crack field
499 : // Util::RealFillBoundary(*crack.c_old_mf[a_lev],geom[a_lev]);
500 0 : crack.c_mf[a_lev]->FillBoundary();
501 : // return;
502 0 : std::swap(crack.c_old_mf[a_lev], crack.c_mf[a_lev]);
503 :
504 0 : const Set::Scalar *DX = geom[a_lev].CellSize();
505 0 : amrex::Box domain(geom[a_lev].Domain());
506 0 : domain.convert(amrex::IntVect::TheNodeVector());
507 0 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
508 :
509 : // Trying out the predictor corrector approach.
510 : //======================================================================
511 : // Predictor step
512 : //======================================================================
513 0 : for (amrex::MFIter mfi(*crack.c_mf[a_lev], true); mfi.isValid(); ++mfi)
514 : {
515 0 : amrex::Box bx = mfi.tilebox();
516 : // bx.grow(1);
517 : // bx = bx & domain;
518 :
519 0 : amrex::Array4<const Set::Scalar> const &c_old = crack.c_old_mf[a_lev]->array(mfi);
520 0 : amrex::Array4<const Set::Scalar> const &energy = crack.history_var_mf[a_lev]->array(mfi);
521 0 : amrex::Array4<const Set::Scalar> const &eta = material.eta_mf[a_lev]->array(mfi);
522 :
523 0 : amrex::Array4<Set::Scalar> const &c = crack.c_mf[a_lev]->array(mfi);
524 0 : amrex::Array4<Set::Scalar> const &df = crack.driving_force_mf[a_lev]->array(mfi);
525 :
526 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
527 : #if AMREX_SPACEDIM != 2
528 : Util::Abort(INFO, "This does not work for 1D or 3D yet.");
529 : #endif
530 :
531 0 : if (i == lo.x && j == lo.y)
532 0 : c(i, j, k, 0) = c(i + 1, j + 1, k, 0);
533 0 : else if (i == lo.x && j == hi.y)
534 0 : c(i, j, k, 0) = c(i + 1, j - 1, k, 0);
535 0 : else if (i == hi.x && j == lo.y)
536 0 : c(i, j, k, 0) = c(i - 1, j + 1, k, 0);
537 0 : else if (i == hi.x && j == hi.y)
538 0 : c(i, j, k, 0) = c(i - 1, j - 1, k, 0);
539 0 : else if (i == lo.x)
540 0 : c(i, j, k) = c(i + 1, j, k, 0);
541 0 : else if (j == lo.y)
542 0 : c(i, j, k) = c(i, j + 1, k, 0);
543 0 : else if (i == hi.x)
544 0 : c(i, j, k) = c(i - 1, j, k, 0);
545 0 : else if (j == hi.y)
546 0 : c(i, j, k) = c(i, j - 1, k, 0);
547 :
548 : // Next, we are doing crack evolution
549 : else
550 : {
551 0 : Set::Scalar rhs = 0.0;
552 0 : Set::Scalar bilap = 0.0;
553 0 : if (crack.beta > 0.0)
554 : {
555 0 : bilap = Numeric::Stencil<Set::Scalar, 4, 0, 0>::D(c_old, i, j, k, 0, DX) + Numeric::Stencil<Set::Scalar, 2, 2, 0>::D(c_old, i, j, k, 0, DX) * 2.0 + Numeric::Stencil<Set::Scalar, 0, 4, 0>::D(c_old, i, j, k, 0, DX);
556 : }
557 :
558 : // if (std::isnan(bilap)) Util::Message(INFO, "Bilaplacian is nan at (",i,", ",j,")");
559 :
560 0 : Set::Scalar laplacian = Numeric::Laplacian(c_old, i, j, k, 0, DX);
561 0 : if (std::isnan(laplacian))
562 0 : Util::Message(INFO, "Laplacian is nan at (", i, ", ", j, ")");
563 :
564 0 : Set::Scalar Gc = 0.0;
565 0 : Set::Scalar Zeta = 0.0;
566 0 : Set::Scalar Threshold = 0.0;
567 0 : Set::Scalar Mobility = 0.0;
568 0 : Set::Scalar _temp_product = 1.0, _temp_product2 = 1.0;
569 :
570 0 : for (int m = 0; m < material.num_mat; m++)
571 : {
572 0 : Gc += eta(i, j, k, m) * crack.cracktype[m].Gc(c_old(i, j, k, 0));
573 0 : Zeta += eta(i, j, k, m) * crack.cracktype[m].Zeta(c_old(i, j, k, 0));
574 0 : Threshold += eta(i, j, k, m) * crack.cracktype[m].DrivingForceThreshold(c_old(i, j, k, 0));
575 0 : Mobility += eta(i, j, k, m) * crack.cracktype[m].Mobility(c_old(i, j, k, 0));
576 0 : _temp_product *= eta(i, j, k, m);
577 0 : _temp_product2 *= 0.5;
578 : }
579 0 : if (material.num_mat > 1)
580 0 : Gc *= (1.0 - _temp_product * (1. - crack.mult_Gc) / _temp_product2);
581 :
582 : // =====================================================================
583 : // Set::Scalar en_cell = energy(i, j, k, 0); // set but not used
584 0 : if (!crack.cracktype[0].mixed_mode())
585 : {
586 0 : df(i, j, k, 0) = crack.cracktype[0].Dg_phi(c_old(i, j, k)) * (energy(i, j, k, 0) + energy(i, j, k, 1)) * crack.el_mult / Gc;
587 : }
588 : else
589 : {
590 0 : df(i, j, k, 0) = crack.cracktype[0].Dg_phi(c_old(i, j, k)) * (energy(i, j, k, 0) + energy(i, j, k, 1)); // * crack.el_mult;
591 : }
592 0 : rhs += df(i, j, k, 0);
593 :
594 0 : df(i, j, k, 1) = crack.cracktype[0].Dw_phi(c_old(i, j, k, 0), 0.0) / (Zeta);
595 0 : rhs += df(i, j, k, 1);
596 :
597 0 : df(i, j, k, 2) = 2.0 * Zeta * laplacian * crack.mult_lap;
598 0 : if (std::isnan(df(i, j, k, 2)))
599 0 : Util::Message(INFO, "DF(2) is nan at (", i, ",", j, ")");
600 0 : rhs -= df(i, j, k, 2);
601 :
602 0 : df(i, j, k, 3) = crack.beta * (0.5 * Zeta * Zeta * Zeta) * bilap;
603 0 : if (std::isnan(df(i, j, k, 3)))
604 0 : Util::Message(INFO, "DF(3) is nan at (", i, ",", j, ")");
605 0 : rhs += df(i, j, k, 3);
606 :
607 0 : df(i, j, k, 4) = std::max(0., rhs - Threshold);
608 0 : c(i, j, k, 0) = c_old(i, j, k, 0) - a_dt * df(i, j, k, 4) * Mobility; //*(4.*c_old(i,j,k,0) - 4.*c_old(i,j,k,0)*c_old(i,j,k,0));
609 :
610 0 : if (c(i, j, k, 0) < 0.0)
611 0 : c(i, j, k, 0) = 0.0;
612 0 : if (c(i, j, k, 0) > 1.0)
613 0 : c(i, j, k, 0) = 1.0;
614 : }
615 0 : });
616 0 : }
617 0 : crack.c_mf[a_lev]->FillBoundary();
618 0 : crack.driving_force_mf[a_lev]->FillBoundary();
619 :
620 0 : Base::Mechanics<brittle_model>::Advance(a_lev, a_time, a_dt);
621 0 : }
622 :
623 : void
624 0 : TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, Set::Scalar a_time, int a_ngrow) override
625 : {
626 0 : Base::Mechanics<brittle_model>::TagCellsForRefinement(lev, a_tags, a_time, a_ngrow);
627 :
628 0 : Set::Vector DX(geom[lev].CellSize());
629 0 : Set::Scalar DXnorm = DX.lpNorm<2>();
630 :
631 0 : for (amrex::MFIter mfi(*crack.c_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
632 : {
633 0 : amrex::Box bx = mfi.tilebox();
634 0 : bx.convert(amrex::IntVect::TheCellVector());
635 0 : amrex::Array4<char> const &tags = a_tags.array(mfi);
636 0 : amrex::Array4<Set::Scalar> const &c = crack.c_mf[lev]->array(mfi);
637 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
638 0 : Set::Vector grad = Numeric::NodeGradientOnCell(c, i, j, k, DX.data());
639 0 : if (grad.lpNorm<2>() * DXnorm > crack.refinement_threshold)
640 0 : tags(i, j, k) = amrex::TagBox::SET;
641 0 : });
642 0 : }
643 :
644 0 : for (amrex::MFIter mfi(*crack.driving_force_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
645 : {
646 0 : amrex::Box bx = mfi.tilebox();
647 0 : bx.convert(amrex::IntVect::TheCellVector());
648 0 : amrex::Array4<char> const &tags = a_tags.array(mfi);
649 0 : amrex::Array4<Set::Scalar> const &df = crack.driving_force_mf[lev]->array(mfi);
650 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
651 0 : Set::Vector grad = Numeric::Gradient(df, i, j, k, 0, DX.data());
652 0 : if (grad.lpNorm<2>() * DXnorm > crack.driving_force_refinement_threshold)
653 0 : tags(i, j, k) = amrex::TagBox::SET;
654 0 : });
655 0 : }
656 :
657 0 : for (amrex::MFIter mfi(*psi_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
658 : {
659 0 : amrex::Box bx = mfi.tilebox();
660 0 : bx.convert(amrex::IntVect::TheCellVector());
661 0 : amrex::Array4<char> const &tags = a_tags.array(mfi);
662 0 : amrex::Array4<Set::Scalar> const &psi = psi_mf[lev]->array(mfi);
663 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
664 0 : Set::Vector grad = Numeric::Gradient(psi, i, j, k, 0, DX.data());
665 0 : if (grad.lpNorm<2>() * DXnorm > crack.refinement_threshold)
666 0 : tags(i, j, k) = amrex::TagBox::SET;
667 0 : });
668 0 : }
669 :
670 0 : if (material.num_mat > 1)
671 : {
672 0 : for (amrex::MFIter mfi(*material.eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
673 : {
674 0 : amrex::Box bx = mfi.tilebox();
675 0 : bx.convert(amrex::IntVect::TheCellVector());
676 0 : amrex::Array4<char> const &tags = a_tags.array(mfi);
677 0 : amrex::Array4<Set::Scalar> const &eta = material.eta_mf[lev]->array(mfi);
678 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
679 0 : Set::Matrix grad = Numeric::Gradient(eta, i, j, k, DX.data());
680 0 : if (grad.lpNorm<2>() * DXnorm > material.m_eta_ref_threshold)
681 0 : tags(i, j, k) = amrex::TagBox::SET;
682 0 : });
683 0 : }
684 : }
685 0 : }
686 :
687 : void
688 0 : Integrate(int amrlev, Set::Scalar time, int step, const amrex::MFIter &mfi, const amrex::Box &a_box) override
689 : {
690 0 : Set::Vector DX(geom[amrlev].CellSize());
691 0 : const Set::Scalar DV = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
692 0 : amrex::Array4<const Set::Scalar> const &df = crack.driving_force_mf[amrlev]->array(mfi);
693 : // amrex::Array4<const Set::Scalar> const &c = crack.c_mf[amrlev]->array(mfi); // not used
694 : // amrex::Array4<const Set::Scalar> const &c_old = crack.c_old_mf[amrlev]->array(mfi); // not used
695 0 : amrex::ParallelFor(a_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
696 0 : crack.driving_force_norm += df(i, j, k, 4) * DV;
697 0 : });
698 :
699 0 : Base::Mechanics<brittle_model>::Integrate(amrlev, time, step, mfi, a_box);
700 0 : }
701 :
702 : void
703 0 : TimeStepComplete(Set::Scalar /*time*/, int /* iter*/) override
704 : {
705 0 : if (elastic_do_solve_now)
706 0 : crack.driving_force_reference = crack.driving_force_norm;
707 :
708 0 : if (set_initial_df_once)
709 : {
710 0 : Util::Message(INFO, "First setting the reference value");
711 0 : crack.driving_force_reference_prev = crack.driving_force_reference;
712 0 : set_initial_df_once = false;
713 : }
714 :
715 0 : if (crack.driving_force_reference / crack.driving_force_reference_prev < 0.1)
716 : {
717 0 : increase_load_step = true;
718 : // crack.driving_force_reference_prev = crack.driving_force_reference;
719 : }
720 0 : if (set_new_reference)
721 : {
722 0 : crack.driving_force_reference_prev = crack.driving_force_reference;
723 0 : set_new_reference = false;
724 : }
725 :
726 0 : elastic_do_solve_now = false;
727 0 : }
728 :
729 : // Member variables
730 : protected:
731 : struct
732 : {
733 : Set::Field<Set::Scalar> c_mf;
734 : Set::Field<Set::Scalar> c_old_mf;
735 : Set::Field<Set::Scalar> energy_pristine_mf;
736 : Set::Field<Set::Scalar> history_var_mf;
737 : Set::Field<Set::Scalar> driving_force_mf;
738 : std::vector<pfczm_crack_type> cracktype;
739 :
740 : std::vector<std::string> ic_type;
741 : std::vector<IC::IC<Set::Scalar> *> ic;
742 : bool is_ic = true;
743 :
744 : Set::Scalar scaleModulusMax = 0.02;
745 : Set::Scalar refinement_threshold = NAN;
746 : Set::Scalar mult_lap = NAN;
747 : Set::Scalar mult_Gc = NAN;
748 : Set::Scalar el_mult = NAN;
749 :
750 : Set::Scalar driving_force_reference = 1.0;
751 : Set::Scalar driving_force_reference_prev = 1.0;
752 : Set::Scalar driving_force_norm = 0.0;
753 : Set::Scalar driving_force_refinement_threshold = NAN;
754 : Set::Scalar tol_rel = NAN;
755 : Set::Scalar tol_abs = NAN;
756 : Set::Scalar max_iter = NAN;
757 :
758 : Set::Scalar crack_l2_err = 0.0;
759 : Set::Scalar crack_norm = 0.0;
760 : Set::Scalar crack_prop_iter = 0;
761 :
762 : Set::Scalar beta = NAN;
763 : } crack;
764 : struct
765 : {
766 : Set::Field<Set::Scalar> eta_mf;
767 : Set::Scalar m_eta_ref_threshold = 0.01;
768 : std::vector<brittle_model> models;
769 : IC::IC<Set::Scalar> *ic;
770 : bool is_ic = false;
771 : int num_mat = 1;
772 : } material;
773 :
774 : bool elastic_do_solve_now = true;
775 : bool increase_load_step = true;
776 : bool set_initial_df_once = false;
777 : bool set_new_reference = true;
778 : int loadstep = 0;
779 :
780 : BC::BC<Set::Scalar> *bc_psi = nullptr;
781 :
782 : using Base::Mechanics<brittle_model>::m_type;
783 : using Base::Mechanics<brittle_model>::finest_level;
784 : using Base::Mechanics<brittle_model>::geom;
785 : using Base::Mechanics<brittle_model>::model_mf;
786 : using Base::Mechanics<brittle_model>::psi_mf;
787 : using Base::Mechanics<brittle_model>::psi_on;
788 : };
789 : }
790 :
791 : #endif
|
