Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "InitialConditionTempl.h"
11 : #include "FEProblem.h"
12 : #include "Assembly.h"
13 : #include "MooseVariableFE.h"
14 : #include "SystemBase.h"
15 :
16 : #include "libmesh/fe_interface.h"
17 : #include "libmesh/quadrature.h"
18 :
19 : using namespace libMesh;
20 :
21 : template <typename T>
22 29188 : InitialConditionTempl<T>::InitialConditionTempl(const InputParameters & parameters)
23 : : InitialConditionBase(parameters),
24 29188 : _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
25 29188 : _tid(getParam<THREAD_ID>("_tid")),
26 29188 : _t(_fe_problem.time()),
27 29188 : _var(_sys.getActualFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
28 29188 : parameters.get<VariableName>("variable"))),
29 29188 : _fe_var(dynamic_cast<MooseVariableFE<T> *>(&_var)),
30 29188 : _assembly(
31 29188 : _fe_problem.assembly(_tid, _var.kind() == Moose::VAR_SOLVER ? _var.sys().number() : 0)),
32 29188 : _coord_sys(_assembly.coordSystem()),
33 29188 : _current_elem(_var.currentElem()),
34 29188 : _current_elem_volume(_assembly.elemVolume()),
35 29188 : _current_node(nullptr),
36 29188 : _qp(0),
37 58376 : _fe_type(_var.feType())
38 : {
39 29188 : }
40 :
41 : template <typename T>
42 25939 : InitialConditionTempl<T>::~InitialConditionTempl()
43 : {
44 25939 : }
45 :
46 : template <typename T>
47 : void
48 2216820 : InitialConditionTempl<T>::compute()
49 : {
50 : // -- NOTE ----
51 : // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
52 : // can couple variable values
53 : // and we also do not call any callbacks, but we use our initial condition system directly.
54 : // ------------
55 :
56 : // The dimension of the current element
57 2216820 : _dim = _current_elem->dim();
58 : // The number of nodes on the new element
59 2216820 : const unsigned int n_nodes = _current_elem->n_nodes();
60 :
61 : // Get FE objects of the appropriate type
62 : // We cannot use the FE object in Assembly, since the following code is messing with the
63 : // quadrature rules
64 : // for projections and would screw it up. However, if we implement projections from one mesh to
65 : // another,
66 : // this code should use that implementation.
67 2216820 : std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));
68 :
69 : // Prepare variables for projection
70 2216820 : std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
71 2216820 : std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
72 2216820 : std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));
73 :
74 : // The values of the shape functions at the quadrature points
75 2216820 : _phi = &fe->get_phi();
76 :
77 : // The gradients of the shape functions at the quadrature points on the child element.
78 2216820 : _dphi = nullptr;
79 :
80 2216820 : _cont = fe->get_continuity();
81 :
82 2216820 : if (_cont == C_ONE)
83 : {
84 3391 : const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
85 3391 : _dphi = &ref_dphi;
86 : }
87 :
88 : // The Jacobian * quadrature weight at the quadrature points
89 2216820 : _JxW = &fe->get_JxW();
90 : // The XYZ locations of the quadrature points
91 2216820 : _xyz_values = &fe->get_xyz();
92 :
93 : // Update the DOF indices for this element based on the current mesh
94 2216820 : _var.prepareIC();
95 2216820 : _dof_indices = _var.dofIndices();
96 :
97 : // The number of DOFs on the element
98 2216820 : const unsigned int n_dofs = _dof_indices.size();
99 2216820 : if (n_dofs == 0)
100 2254 : return;
101 :
102 : // Fixed vs. free DoFs on edge/face projections
103 2214566 : _dof_is_fixed.clear();
104 2214566 : _dof_is_fixed.resize(n_dofs, false);
105 2214566 : _free_dof.clear();
106 2214566 : _free_dof.resize(n_dofs, 0);
107 :
108 : // Zero the interpolated values
109 2214566 : _Ue.resize(n_dofs);
110 2214566 : _Ue.zero();
111 :
112 2214566 : DenseVector<char> mask(n_dofs, true);
113 :
114 : // In general, we need a series of
115 : // projections to ensure a unique and continuous
116 : // solution. We start by interpolating nodes, then
117 : // hold those fixed and project edges, then
118 : // hold those fixed and project faces, then
119 : // hold those fixed and project interiors
120 :
121 : // Interpolate node values first
122 2214566 : _current_dof = 0;
123 :
124 15681329 : for (_n = 0; _n != n_nodes; ++_n)
125 : {
126 13466767 : _nc = FEInterface::n_dofs_at_node(_fe_type, _current_elem, _n);
127 :
128 : // for nodes that are in more than one subdomain, only compute the initial
129 : // condition once on the lowest numbered block
130 13466767 : auto curr_node = _current_elem->node_ptr(_n);
131 13466767 : const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);
132 :
133 13466767 : auto priority_block = *(block_ids.begin());
134 13474202 : for (auto id : block_ids)
135 13474202 : if (_var.hasBlocks(id))
136 : {
137 13466767 : priority_block = id;
138 13466767 : break;
139 : }
140 :
141 13466767 : if (!hasBlocks(priority_block) && _var.isNodal())
142 : {
143 71056 : for (decltype(_nc) i = 0; i < _nc; ++i)
144 : {
145 36812 : mask(_current_dof) = false;
146 36812 : _current_dof++;
147 : }
148 34244 : continue;
149 34244 : }
150 :
151 : // FIXME: this should go through the DofMap,
152 : // not duplicate _dof_indices code badly!
153 13432523 : if (!_current_elem->is_vertex(_n))
154 : {
155 3256794 : _current_dof += _nc;
156 3256794 : continue;
157 : }
158 :
159 10175729 : if (_cont == DISCONTINUOUS || _cont == H_CURL || _cont == H_DIV)
160 : libmesh_assert(_nc == 0);
161 5778855 : else if (_cont == C_ZERO)
162 5765051 : setCZeroVertices();
163 13804 : else if (_fe_type.family == HERMITE)
164 13548 : setHermiteVertices();
165 256 : else if (_cont == C_ONE)
166 0 : setOtherCOneVertices();
167 256 : else if (_cont == SIDE_DISCONTINUOUS)
168 256 : continue;
169 : else
170 0 : libmesh_error();
171 : } // loop over nodes
172 :
173 : // From here on out we won't be sampling at nodes anymore
174 2214562 : _current_node = nullptr;
175 :
176 2214562 : auto & dof_map = _var.dofMap();
177 : const bool add_p_level =
178 2214562 : dof_map.should_p_refine(dof_map.var_group_from_var_number(_var.number()));
179 :
180 : // In 3D, project any edge values next
181 2214562 : if (_dim > 2 && _cont != DISCONTINUOUS)
182 4176518 : for (unsigned int e = 0; e != _current_elem->n_edges(); ++e)
183 : {
184 3854856 : FEInterface::dofs_on_edge(_current_elem, _dim, _fe_type, e, _side_dofs, add_p_level);
185 :
186 : // Some edge dofs are on nodes and already
187 : // fixed, others are free to calculate
188 3854856 : _free_dofs = 0;
189 11842488 : for (unsigned int i = 0; i != _side_dofs.size(); ++i)
190 7987632 : if (!_dof_is_fixed[_side_dofs[i]])
191 219936 : _free_dof[_free_dofs++] = i;
192 :
193 : // There may be nothing to project
194 3854856 : if (!_free_dofs)
195 3635432 : continue;
196 :
197 : // Initialize FE data on the edge
198 219424 : fe->attach_quadrature_rule(qedgerule.get());
199 219424 : fe->edge_reinit(_current_elem, e);
200 219424 : _n_qp = qedgerule->n_points();
201 :
202 219424 : choleskySolve(false);
203 : }
204 :
205 : // Project any side values (edges in 2D, faces in 3D)
206 2214562 : if (_dim > 1 && _cont != DISCONTINUOUS)
207 6221491 : for (unsigned int s = 0; s != _current_elem->n_sides(); ++s)
208 : {
209 5102952 : FEInterface::dofs_on_side(_current_elem, _dim, _fe_type, s, _side_dofs, add_p_level);
210 :
211 : // Some side dofs are on nodes/edges and already
212 : // fixed, others are free to calculate
213 5102952 : _free_dofs = 0;
214 20115772 : for (unsigned int i = 0; i != _side_dofs.size(); ++i)
215 15012820 : if (!_dof_is_fixed[_side_dofs[i]])
216 403576 : _free_dof[_free_dofs++] = i;
217 :
218 : // There may be nothing to project
219 5102952 : if (!_free_dofs)
220 4707809 : continue;
221 :
222 : // Initialize FE data on the side
223 395143 : fe->attach_quadrature_rule(qsiderule.get());
224 395143 : fe->reinit(_current_elem, s);
225 395143 : _n_qp = qsiderule->n_points();
226 :
227 395143 : choleskySolve(false);
228 : }
229 :
230 : // Project the interior values, finally
231 :
232 : // Some interior dofs are on nodes/edges/sides and
233 : // already fixed, others are free to calculate
234 2214562 : _free_dofs = 0;
235 10002772 : for (unsigned int i = 0; i != n_dofs; ++i)
236 7788210 : if (!_dof_is_fixed[i])
237 1322459 : _free_dof[_free_dofs++] = i;
238 :
239 : // There may be nothing to project
240 2214562 : if (_free_dofs)
241 : {
242 : // Initialize FE data
243 1143287 : fe->attach_quadrature_rule(qrule.get());
244 1143287 : fe->reinit(_current_elem);
245 1143287 : _n_qp = qrule->n_points();
246 :
247 1143287 : choleskySolve(true);
248 : } // if there are free interior dofs
249 :
250 : // Make sure every DoF got reached!
251 10002766 : for (unsigned int i = 0; i != n_dofs; ++i)
252 : libmesh_assert(_dof_is_fixed[i]);
253 :
254 : // Lock the new_vector since it is shared among threads.
255 : {
256 2214559 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
257 10002766 : for (size_t i = 0; i < mask.size(); i++)
258 7788207 : if (mask(i))
259 7751395 : _var.setDofValue(_Ue(i), i);
260 2214559 : }
261 2223575 : }
262 :
263 : template <typename T>
264 : void
265 5752547 : InitialConditionTempl<T>::setCZeroVertices()
266 : {
267 : // Assume that C_ZERO elements have a single nodal
268 : // value shape function
269 : libmesh_assert(_nc == 1);
270 5752547 : _qp = _n;
271 5752547 : _current_node = _current_elem->node_ptr(_n);
272 5752547 : _Ue(_current_dof) = value(*_current_node);
273 5752543 : _dof_is_fixed[_current_dof] = true;
274 5752543 : _current_dof++;
275 5752543 : }
276 :
277 : template <>
278 : void
279 12504 : InitialConditionTempl<RealVectorValue>::setCZeroVertices()
280 : {
281 12504 : _qp = _n;
282 12504 : _current_node = _current_elem->node_ptr(_n);
283 12504 : auto point_value = value(*_current_node);
284 48040 : for (decltype(_nc) i = 0; i < _nc; ++i)
285 : {
286 35536 : _Ue(_current_dof) = point_value(i);
287 35536 : _dof_is_fixed[_current_dof] = true;
288 35536 : _current_dof++;
289 : }
290 12504 : }
291 :
292 : template <typename T>
293 : T
294 54144 : InitialConditionTempl<T>::gradientComponent(GradientType grad, unsigned int i)
295 : {
296 54144 : return grad(i);
297 : }
298 :
299 : template <>
300 : RealVectorValue
301 0 : InitialConditionTempl<RealVectorValue>::gradientComponent(GradientType grad, unsigned int i)
302 : {
303 0 : return grad.row(i);
304 : }
305 :
306 : template <>
307 : RealEigenVector
308 0 : InitialConditionTempl<RealEigenVector>::gradientComponent(GradientType grad, unsigned int i)
309 : {
310 0 : return grad.col(i);
311 : }
312 :
313 : template <typename T>
314 : void
315 13548 : InitialConditionTempl<T>::setHermiteVertices()
316 : {
317 : // The hermite element vertex shape functions are weird
318 13548 : _qp = _n;
319 13548 : _current_node = _current_elem->node_ptr(_n);
320 13548 : _Ue(_current_dof) = value(*_current_node);
321 13548 : _dof_is_fixed[_current_dof] = true;
322 13548 : _current_dof++;
323 13548 : GradientType grad = gradient(*_current_node);
324 : // x derivative
325 13548 : _Ue(_current_dof) = gradientComponent(grad, 0);
326 13548 : _dof_is_fixed[_current_dof] = true;
327 13548 : _current_dof++;
328 13548 : if (_dim > 1)
329 : {
330 : // We'll finite difference mixed derivatives
331 13532 : Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n);
332 13532 : nxminus(0) -= TOLERANCE;
333 13532 : nxplus(0) += TOLERANCE;
334 13532 : GradientType gxminus = gradient(nxminus);
335 13532 : GradientType gxplus = gradient(nxplus);
336 : // y derivative
337 13532 : _Ue(_current_dof) = gradientComponent(grad, 1);
338 13532 : _dof_is_fixed[_current_dof] = true;
339 13532 : _current_dof++;
340 : // xy derivative
341 27064 : _Ue(_current_dof) =
342 13532 : (gradientComponent(gxplus, 1) - gradientComponent(gxminus, 1)) / 2. / TOLERANCE;
343 13532 : _dof_is_fixed[_current_dof] = true;
344 13532 : _current_dof++;
345 :
346 13532 : if (_dim > 2)
347 : {
348 : // z derivative
349 0 : _Ue(_current_dof) = gradientComponent(grad, 2);
350 0 : _dof_is_fixed[_current_dof] = true;
351 0 : _current_dof++;
352 : // xz derivative
353 0 : _Ue(_current_dof) =
354 0 : (gradientComponent(gxplus, 2) - gradientComponent(gxminus, 2)) / 2. / TOLERANCE;
355 0 : _dof_is_fixed[_current_dof] = true;
356 0 : _current_dof++;
357 : // We need new points for yz
358 0 : Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n);
359 0 : nyminus(1) -= TOLERANCE;
360 0 : nyplus(1) += TOLERANCE;
361 0 : GradientType gyminus = gradient(nyminus);
362 0 : GradientType gyplus = gradient(nyplus);
363 : // xz derivative
364 0 : _Ue(_current_dof) =
365 0 : (gradientComponent(gyplus, 2) - gradientComponent(gyminus, 2)) / 2. / TOLERANCE;
366 0 : _dof_is_fixed[_current_dof] = true;
367 0 : _current_dof++;
368 : // Getting a 2nd order xyz is more tedious
369 0 : Point nxmym = _current_elem->point(_n), nxmyp = _current_elem->point(_n),
370 0 : nxpym = _current_elem->point(_n), nxpyp = _current_elem->point(_n);
371 0 : nxmym(0) -= TOLERANCE;
372 0 : nxmym(1) -= TOLERANCE;
373 0 : nxmyp(0) -= TOLERANCE;
374 0 : nxmyp(1) += TOLERANCE;
375 0 : nxpym(0) += TOLERANCE;
376 0 : nxpym(1) -= TOLERANCE;
377 0 : nxpyp(0) += TOLERANCE;
378 0 : nxpyp(1) += TOLERANCE;
379 0 : GradientType gxmym = gradient(nxmym);
380 0 : GradientType gxmyp = gradient(nxmyp);
381 0 : GradientType gxpym = gradient(nxpym);
382 0 : GradientType gxpyp = gradient(nxpyp);
383 0 : DataType gxzplus =
384 0 : (gradientComponent(gxpyp, 2) - gradientComponent(gxmyp, 2)) / 2. / TOLERANCE;
385 0 : DataType gxzminus =
386 0 : (gradientComponent(gxpym, 2) - gradientComponent(gxmym, 2)) / 2. / TOLERANCE;
387 : // xyz derivative
388 0 : _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE;
389 0 : _dof_is_fixed[_current_dof] = true;
390 0 : _current_dof++;
391 0 : }
392 0 : }
393 13548 : }
394 :
395 : template <>
396 : void
397 0 : InitialConditionTempl<RealVectorValue>::setHermiteVertices()
398 : {
399 0 : }
400 :
401 : template <typename T>
402 : void
403 0 : InitialConditionTempl<T>::setOtherCOneVertices()
404 : {
405 : // Assume that other C_ONE elements have a single nodal
406 : // value shape function and nodal gradient component
407 : // shape functions
408 : libmesh_assert(_nc == 1 + _dim);
409 0 : _current_node = _current_elem->node_ptr(_n);
410 0 : _Ue(_current_dof) = value(*_current_node);
411 0 : _dof_is_fixed[_current_dof] = true;
412 0 : _current_dof++;
413 0 : GradientType grad = gradient(*_current_node);
414 0 : for (unsigned int i = 0; i != _dim; ++i)
415 : {
416 0 : _Ue(_current_dof) = gradientComponent(grad, i);
417 0 : _dof_is_fixed[_current_dof] = true;
418 0 : _current_dof++;
419 : }
420 0 : }
421 :
422 : template <>
423 : void
424 0 : InitialConditionTempl<RealVectorValue>::setOtherCOneVertices()
425 : {
426 0 : }
427 :
428 : template <typename T>
429 : void
430 1757854 : InitialConditionTempl<T>::choleskyAssembly(bool is_volume)
431 : {
432 : // Loop over the quadrature points
433 6632510 : for (_qp = 0; _qp < _n_qp; _qp++)
434 : {
435 : // solution at the quadrature point
436 4874659 : auto fineval = value((*_xyz_values)[_qp]);
437 : // solution grad at the quadrature point
438 4874656 : GradientType finegrad;
439 4874656 : if (_cont == C_ONE)
440 0 : finegrad = gradient((*_xyz_values)[_qp]);
441 :
442 4874656 : auto dofs_size = is_volume ? _dof_indices.size() : _side_dofs.size();
443 :
444 : // Form edge projection matrix
445 39326512 : for (decltype(dofs_size) geomi = 0, freei = 0; geomi != dofs_size; ++geomi)
446 : {
447 34451856 : auto i = is_volume ? geomi : _side_dofs[geomi];
448 :
449 : // fixed DoFs aren't test functions
450 34451856 : if (_dof_is_fixed[i])
451 27936282 : continue;
452 53113890 : for (decltype(dofs_size) geomj = 0, freej = 0; geomj != dofs_size; ++geomj)
453 : {
454 46598316 : auto j = is_volume ? geomj : _side_dofs[geomj];
455 46598316 : if (_dof_is_fixed[j])
456 27941914 : _Fe(freei) -= (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp] * _Ue(j);
457 : else
458 18656402 : _Ke(freei, freej) += (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp];
459 46598316 : if (_cont == C_ONE)
460 : {
461 0 : if (_dof_is_fixed[j])
462 0 : _Fe(freei) -= dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp] * _Ue(j);
463 : else
464 0 : _Ke(freei, freej) += dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp];
465 : }
466 46598316 : if (!_dof_is_fixed[j])
467 18656402 : freej++;
468 : }
469 6515574 : _Fe(freei) += (*_phi)[i][_qp] * fineval * (*_JxW)[_qp];
470 6515574 : if (_cont == C_ONE)
471 0 : _Fe(freei) += dotHelper(finegrad, (*_dphi)[i][_qp]) * (*_JxW)[_qp];
472 6515574 : freei++;
473 : }
474 : }
475 1757851 : }
476 :
477 : template <typename T>
478 : void
479 1729838 : InitialConditionTempl<T>::choleskySolve(bool is_volume)
480 : {
481 1729838 : _Ke.resize(_free_dofs, _free_dofs);
482 1729838 : _Ke.zero();
483 1729838 : _Fe.resize(_free_dofs);
484 1729838 : _Fe.zero();
485 :
486 1729838 : choleskyAssembly(is_volume);
487 :
488 : // The new edge coefficients
489 1729835 : DenseVector<DataType> U(_free_dofs);
490 :
491 1729835 : _Ke.cholesky_solve(_Fe, U);
492 :
493 : // Transfer new edge solutions to element
494 3640059 : for (unsigned int i = 0; i != _free_dofs; ++i)
495 : {
496 1910224 : auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
497 1910224 : DataType & ui = _Ue(the_dof);
498 : libmesh_assert(std::abs(ui) < TOLERANCE || std::abs(ui - U(i)) < TOLERANCE);
499 1910224 : ui = U(i);
500 1910224 : _dof_is_fixed[the_dof] = true;
501 : }
502 1729835 : }
503 :
504 : template <>
505 : void
506 28016 : InitialConditionTempl<RealEigenVector>::choleskySolve(bool is_volume)
507 : {
508 28016 : _Ke.resize(_free_dofs, _free_dofs);
509 28016 : _Ke.zero();
510 28016 : _Fe.resize(_free_dofs);
511 63760 : for (unsigned int i = 0; i < _free_dofs; ++i)
512 35744 : _Fe(i).setZero(_var.count());
513 :
514 28016 : choleskyAssembly(is_volume);
515 :
516 : // The new edge coefficients
517 28016 : DenseVector<DataType> U = _Fe;
518 :
519 85840 : for (unsigned int i = 0; i < _var.count(); ++i)
520 : {
521 57824 : DenseVector<Real> v(_free_dofs), x(_free_dofs);
522 133024 : for (unsigned int j = 0; j < _free_dofs; ++j)
523 75200 : v(j) = _Fe(j)(i);
524 :
525 57824 : _Ke.cholesky_solve(v, x);
526 :
527 133024 : for (unsigned int j = 0; j < _free_dofs; ++j)
528 75200 : U(j)(i) = x(j);
529 57824 : }
530 :
531 : // Transfer new edge solutions to element
532 63760 : for (unsigned int i = 0; i != _free_dofs; ++i)
533 : {
534 35744 : auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
535 35744 : DataType & ui = _Ue(the_dof);
536 : libmesh_assert(ui.matrix().norm() < TOLERANCE || (ui - U(i)).matrix().norm() < TOLERANCE);
537 35744 : ui = U(i);
538 35744 : _dof_is_fixed[the_dof] = true;
539 : }
540 28016 : }
541 :
542 : template <typename T>
543 : void
544 80 : InitialConditionTempl<T>::computeNodal(const Point & p)
545 : {
546 80 : _var.reinitNode();
547 80 : _var.computeNodalValues(); // has to call this to resize the internal array
548 80 : auto return_value = value(p);
549 :
550 80 : _var.setNodalValue(return_value); // update variable data, which is referenced by others, so the
551 : // value is up-to-date
552 :
553 : // We are done, so update the solution vector
554 : {
555 80 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
556 80 : _var.insert(_var.sys().solution());
557 80 : }
558 80 : }
559 :
560 : template class InitialConditionTempl<Real>;
561 : template class InitialConditionTempl<RealVectorValue>;
562 : template class InitialConditionTempl<RealEigenVector>;
|
