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 33324 : InitialConditionTempl<T>::InitialConditionTempl(const InputParameters & parameters)
23 : : InitialConditionBase(parameters),
24 99972 : _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
25 66648 : _tid(getParam<THREAD_ID>("_tid")),
26 33324 : _t(_fe_problem.time()),
27 33324 : _var(_sys.getActualFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
28 33324 : parameters.get<VariableName>("variable"))),
29 33324 : _fe_var(dynamic_cast<MooseVariableFE<T> *>(&_var)),
30 33324 : _assembly(
31 33324 : _fe_problem.assembly(_tid, _var.kind() == Moose::VAR_SOLVER ? _var.sys().number() : 0)),
32 33324 : _coord_sys(_assembly.coordSystem()),
33 33324 : _current_elem(_var.currentElem()),
34 33324 : _current_elem_volume(_assembly.elemVolume()),
35 33324 : _current_node(nullptr),
36 33324 : _qp(0),
37 33324 : _fe_type(_var.feType()),
38 99972 : _dof_indices(_var.dofIndices())
39 : {
40 33324 : }
41 :
42 : template <typename T>
43 30306 : InitialConditionTempl<T>::~InitialConditionTempl()
44 : {
45 30306 : }
46 :
47 : template <typename T>
48 : void
49 2312231 : InitialConditionTempl<T>::compute()
50 : {
51 : // -- NOTE ----
52 : // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
53 : // can couple variable values
54 : // and we also do not call any callbacks, but we use our initial condition system directly.
55 : // Eventually we should try to fix things so we're not duplicating code
56 : // ------------
57 :
58 : // The dimension of the current element
59 2312231 : _dim = _current_elem->dim();
60 : // The number of nodes on the new element
61 2312231 : const unsigned int n_nodes = _current_elem->n_nodes();
62 :
63 : // Get FE objects of the appropriate type
64 : // We cannot use the FE object in Assembly, since the following code is messing with the
65 : // quadrature rules
66 : // for projections and would screw it up. However, if we implement projections from one mesh to
67 : // another,
68 : // this code should use that implementation.
69 2312231 : std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));
70 :
71 : // Prepare variables for projection
72 2312231 : std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
73 2312231 : std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
74 2312231 : std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));
75 :
76 : // The values of the shape functions at the quadrature points
77 2312231 : _phi = &fe->get_phi();
78 :
79 : // The gradients of the shape functions at the quadrature points on the child element.
80 2312231 : _dphi = nullptr;
81 :
82 2312231 : _cont = fe->get_continuity();
83 :
84 2312231 : if (_cont == C_ONE)
85 : {
86 3366 : const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
87 3366 : _dphi = &ref_dphi;
88 : }
89 :
90 : // The Jacobian * quadrature weight at the quadrature points
91 2312231 : _JxW = &fe->get_JxW();
92 : // The XYZ locations of the quadrature points
93 2312231 : _xyz_values = &fe->get_xyz();
94 :
95 : // Update the DOF indices for this element based on the current mesh
96 2312231 : _var.prepareIC();
97 :
98 : // The number of DOFs on the element for this finite element type
99 2312231 : const unsigned int n_dofs = _dof_indices.size() / _var.count();
100 : mooseAssert(_dof_indices.size() % _var.count() == 0,
101 : "The number of degrees of freedom should be cleanly divisible by the variable count");
102 :
103 2312231 : if (n_dofs == 0)
104 3155 : return;
105 :
106 : // Fixed vs. free DoFs on edge/face projections
107 2309076 : _dof_is_fixed.clear();
108 2309076 : _dof_is_fixed.resize(n_dofs, false);
109 2309076 : _free_dof.clear();
110 2309076 : _free_dof.resize(n_dofs, 0);
111 :
112 : // Zero the interpolated values
113 2309076 : _Ue.resize(n_dofs);
114 2309076 : _Ue.zero();
115 :
116 2309076 : DenseVector<char> mask(n_dofs, true);
117 :
118 : // In general, we need a series of
119 : // projections to ensure a unique and continuous
120 : // solution. We start by interpolating nodes, then
121 : // hold those fixed and project edges, then
122 : // hold those fixed and project faces, then
123 : // hold those fixed and project interiors
124 :
125 : // Interpolate node values first
126 2309076 : _current_dof = 0;
127 :
128 17564465 : for (_n = 0; _n != n_nodes; ++_n)
129 : {
130 15255392 : _nc = FEInterface::n_dofs_at_node(_fe_type, _current_elem, _n);
131 :
132 : // for nodes that are in more than one subdomain, only compute the initial
133 : // condition once on the lowest numbered block
134 15255392 : auto curr_node = _current_elem->node_ptr(_n);
135 15255392 : const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);
136 :
137 15255392 : auto priority_block = *(block_ids.begin());
138 15263998 : for (auto id : block_ids)
139 15263998 : if (_var.hasBlocks(id))
140 : {
141 15255392 : priority_block = id;
142 15255392 : break;
143 : }
144 :
145 15255392 : if (!hasBlocks(priority_block) && _var.isNodal())
146 : {
147 68564 : for (decltype(_nc) i = 0; i < _nc; ++i)
148 : {
149 35506 : mask(_current_dof) = false;
150 35506 : _current_dof++;
151 : }
152 33058 : continue;
153 33058 : }
154 :
155 15222334 : if (!_current_elem->is_vertex(_n))
156 : {
157 4144441 : _current_dof += _nc;
158 4144441 : continue;
159 : }
160 :
161 11077893 : if (_cont == DISCONTINUOUS || _cont == H_CURL || _cont == H_DIV)
162 : libmesh_assert(_nc == 0);
163 6359477 : else if (_cont == C_ZERO)
164 6345773 : setCZeroVertices();
165 13704 : else if (_fe_type.family == HERMITE)
166 13448 : setHermiteVertices();
167 256 : else if (_cont == C_ONE)
168 0 : setOtherCOneVertices();
169 256 : else if (_cont == SIDE_DISCONTINUOUS)
170 256 : continue;
171 : else
172 0 : libmesh_error();
173 : } // loop over nodes
174 :
175 : // From here on out we won't be sampling at nodes anymore
176 2309073 : _current_node = nullptr;
177 :
178 2309073 : auto & dof_map = _var.dofMap();
179 : const bool add_p_level =
180 2309073 : dof_map.should_p_refine(dof_map.var_group_from_var_number(_var.number()));
181 :
182 : // In 3D, project any edge values next
183 2309073 : if (_dim > 2 && _cont != DISCONTINUOUS)
184 4682868 : for (unsigned int e = 0; e != _current_elem->n_edges(); ++e)
185 : {
186 4322256 : FEInterface::dofs_on_edge(_current_elem, _dim, _fe_type, e, _side_dofs, add_p_level);
187 :
188 : // Some edge dofs are on nodes and already
189 : // fixed, others are free to calculate
190 4322256 : _free_dofs = 0;
191 13304688 : for (unsigned int i = 0; i != _side_dofs.size(); ++i)
192 8982432 : if (!_dof_is_fixed[_side_dofs[i]])
193 279936 : _free_dof[_free_dofs++] = i;
194 :
195 : // There may be nothing to project
196 4322256 : if (!_free_dofs)
197 4042832 : continue;
198 :
199 : // Initialize FE data on the edge
200 279424 : fe->attach_quadrature_rule(qedgerule.get());
201 279424 : fe->edge_reinit(_current_elem, e);
202 279424 : _n_qp = qedgerule->n_points();
203 :
204 279424 : choleskySolve(false);
205 : }
206 :
207 : // Project any side values (edges in 2D, faces in 3D)
208 2309073 : if (_dim > 1 && _cont != DISCONTINUOUS)
209 6817944 : for (unsigned int s = 0; s != _current_elem->n_sides(); ++s)
210 : {
211 5595320 : FEInterface::dofs_on_side(_current_elem, _dim, _fe_type, s, _side_dofs, add_p_level);
212 :
213 : // Some side dofs are on nodes/edges and already
214 : // fixed, others are free to calculate
215 5595320 : _free_dofs = 0;
216 22220140 : for (unsigned int i = 0; i != _side_dofs.size(); ++i)
217 16624820 : if (!_dof_is_fixed[_side_dofs[i]])
218 442992 : _free_dof[_free_dofs++] = i;
219 :
220 : // There may be nothing to project
221 5595320 : if (!_free_dofs)
222 5160347 : continue;
223 :
224 : // Initialize FE data on the side
225 434973 : fe->attach_quadrature_rule(qsiderule.get());
226 434973 : fe->reinit(_current_elem, s);
227 434973 : _n_qp = qsiderule->n_points();
228 :
229 434973 : choleskySolve(false);
230 : }
231 :
232 : // Project the interior values, finally
233 :
234 : // Some interior dofs are on nodes/edges/sides and
235 : // already fixed, others are free to calculate
236 2309073 : _free_dofs = 0;
237 10780385 : for (unsigned int i = 0; i != n_dofs; ++i)
238 8471312 : if (!_dof_is_fixed[i])
239 1325830 : _free_dof[_free_dofs++] = i;
240 :
241 : // There may be nothing to project
242 2309073 : if (_free_dofs)
243 : {
244 : // Initialize FE data
245 1136851 : fe->attach_quadrature_rule(qrule.get());
246 1136851 : fe->reinit(_current_elem);
247 1136851 : _n_qp = qrule->n_points();
248 :
249 1136851 : choleskySolve(true);
250 : } // if there are free interior dofs
251 :
252 : // Make sure every DoF got reached!
253 10780377 : for (unsigned int i = 0; i != n_dofs; ++i)
254 : libmesh_assert(_dof_is_fixed[i]);
255 :
256 10780377 : for (size_t i = 0; i < mask.size(); i++)
257 8471308 : if (mask(i))
258 8435802 : _var.setDofValue(_Ue(i), i);
259 2321689 : }
260 :
261 : template <typename T>
262 : void
263 6333277 : InitialConditionTempl<T>::setCZeroVertices()
264 : {
265 : // Assume that C_ZERO elements have a single nodal
266 : // value shape function
267 : libmesh_assert(_nc == 1);
268 6333277 : _qp = _n;
269 6333277 : _current_node = _current_elem->node_ptr(_n);
270 6333277 : _Ue(_current_dof) = value(*_current_node);
271 6333274 : _dof_is_fixed[_current_dof] = true;
272 6333274 : _current_dof++;
273 6333274 : }
274 :
275 : template <>
276 : void
277 12496 : InitialConditionTempl<RealVectorValue>::setCZeroVertices()
278 : {
279 12496 : _qp = _n;
280 12496 : _current_node = _current_elem->node_ptr(_n);
281 12496 : auto point_value = value(*_current_node);
282 48016 : for (decltype(_nc) i = 0; i < _nc; ++i)
283 : {
284 35520 : _Ue(_current_dof) = point_value(i);
285 35520 : _dof_is_fixed[_current_dof] = true;
286 35520 : _current_dof++;
287 : }
288 12496 : }
289 :
290 : template <typename T>
291 : T
292 53744 : InitialConditionTempl<T>::gradientComponent(GradientType grad, unsigned int i)
293 : {
294 53744 : return grad(i);
295 : }
296 :
297 : template <>
298 : RealVectorValue
299 0 : InitialConditionTempl<RealVectorValue>::gradientComponent(GradientType grad, unsigned int i)
300 : {
301 0 : return grad.row(i);
302 : }
303 :
304 : template <>
305 : RealEigenVector
306 0 : InitialConditionTempl<RealEigenVector>::gradientComponent(GradientType grad, unsigned int i)
307 : {
308 0 : return grad.col(i);
309 : }
310 :
311 : template <typename T>
312 : void
313 13448 : InitialConditionTempl<T>::setHermiteVertices()
314 : {
315 : // The hermite element vertex shape functions are weird
316 13448 : _qp = _n;
317 13448 : _current_node = _current_elem->node_ptr(_n);
318 13448 : _Ue(_current_dof) = value(*_current_node);
319 13448 : _dof_is_fixed[_current_dof] = true;
320 13448 : _current_dof++;
321 13448 : GradientType grad = gradient(*_current_node);
322 : // x derivative
323 13448 : _Ue(_current_dof) = gradientComponent(grad, 0);
324 13448 : _dof_is_fixed[_current_dof] = true;
325 13448 : _current_dof++;
326 13448 : if (_dim > 1)
327 : {
328 : // We'll finite difference mixed derivatives
329 13432 : Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n);
330 13432 : nxminus(0) -= TOLERANCE;
331 13432 : nxplus(0) += TOLERANCE;
332 13432 : GradientType gxminus = gradient(nxminus);
333 13432 : GradientType gxplus = gradient(nxplus);
334 : // y derivative
335 13432 : _Ue(_current_dof) = gradientComponent(grad, 1);
336 13432 : _dof_is_fixed[_current_dof] = true;
337 13432 : _current_dof++;
338 : // xy derivative
339 26864 : _Ue(_current_dof) =
340 13432 : (gradientComponent(gxplus, 1) - gradientComponent(gxminus, 1)) / 2. / TOLERANCE;
341 13432 : _dof_is_fixed[_current_dof] = true;
342 13432 : _current_dof++;
343 :
344 13432 : if (_dim > 2)
345 : {
346 : // z derivative
347 0 : _Ue(_current_dof) = gradientComponent(grad, 2);
348 0 : _dof_is_fixed[_current_dof] = true;
349 0 : _current_dof++;
350 : // xz derivative
351 0 : _Ue(_current_dof) =
352 0 : (gradientComponent(gxplus, 2) - gradientComponent(gxminus, 2)) / 2. / TOLERANCE;
353 0 : _dof_is_fixed[_current_dof] = true;
354 0 : _current_dof++;
355 : // We need new points for yz
356 0 : Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n);
357 0 : nyminus(1) -= TOLERANCE;
358 0 : nyplus(1) += TOLERANCE;
359 0 : GradientType gyminus = gradient(nyminus);
360 0 : GradientType gyplus = gradient(nyplus);
361 : // xz derivative
362 0 : _Ue(_current_dof) =
363 0 : (gradientComponent(gyplus, 2) - gradientComponent(gyminus, 2)) / 2. / TOLERANCE;
364 0 : _dof_is_fixed[_current_dof] = true;
365 0 : _current_dof++;
366 : // Getting a 2nd order xyz is more tedious
367 0 : Point nxmym = _current_elem->point(_n), nxmyp = _current_elem->point(_n),
368 0 : nxpym = _current_elem->point(_n), nxpyp = _current_elem->point(_n);
369 0 : nxmym(0) -= TOLERANCE;
370 0 : nxmym(1) -= TOLERANCE;
371 0 : nxmyp(0) -= TOLERANCE;
372 0 : nxmyp(1) += TOLERANCE;
373 0 : nxpym(0) += TOLERANCE;
374 0 : nxpym(1) -= TOLERANCE;
375 0 : nxpyp(0) += TOLERANCE;
376 0 : nxpyp(1) += TOLERANCE;
377 0 : GradientType gxmym = gradient(nxmym);
378 0 : GradientType gxmyp = gradient(nxmyp);
379 0 : GradientType gxpym = gradient(nxpym);
380 0 : GradientType gxpyp = gradient(nxpyp);
381 0 : DataType gxzplus =
382 0 : (gradientComponent(gxpyp, 2) - gradientComponent(gxmyp, 2)) / 2. / TOLERANCE;
383 0 : DataType gxzminus =
384 0 : (gradientComponent(gxpym, 2) - gradientComponent(gxmym, 2)) / 2. / TOLERANCE;
385 : // xyz derivative
386 0 : _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE;
387 0 : _dof_is_fixed[_current_dof] = true;
388 0 : _current_dof++;
389 0 : }
390 0 : }
391 13448 : }
392 :
393 : template <>
394 : void
395 0 : InitialConditionTempl<RealVectorValue>::setHermiteVertices()
396 : {
397 0 : }
398 :
399 : template <typename T>
400 : void
401 0 : InitialConditionTempl<T>::setOtherCOneVertices()
402 : {
403 : // Assume that other C_ONE elements have a single nodal
404 : // value shape function and nodal gradient component
405 : // shape functions
406 : libmesh_assert(_nc == 1 + _dim);
407 0 : _current_node = _current_elem->node_ptr(_n);
408 0 : _Ue(_current_dof) = value(*_current_node);
409 0 : _dof_is_fixed[_current_dof] = true;
410 0 : _current_dof++;
411 0 : GradientType grad = gradient(*_current_node);
412 0 : for (unsigned int i = 0; i != _dim; ++i)
413 : {
414 0 : _Ue(_current_dof) = gradientComponent(grad, i);
415 0 : _dof_is_fixed[_current_dof] = true;
416 0 : _current_dof++;
417 : }
418 0 : }
419 :
420 : template <>
421 : void
422 0 : InitialConditionTempl<RealVectorValue>::setOtherCOneVertices()
423 : {
424 0 : }
425 :
426 : template <typename T>
427 : void
428 1851248 : InitialConditionTempl<T>::choleskyAssembly(bool is_volume)
429 : {
430 : // Loop over the quadrature points
431 7364985 : for (_qp = 0; _qp < _n_qp; _qp++)
432 : {
433 : // solution at the quadrature point
434 5513741 : auto fineval = value((*_xyz_values)[_qp]);
435 : // solution grad at the quadrature point
436 5513737 : GradientType finegrad;
437 5513737 : if (_cont == C_ONE)
438 0 : finegrad = gradient((*_xyz_values)[_qp]);
439 :
440 5513737 : auto dofs_size = is_volume ? (_dof_indices.size() / _var.count()) : _side_dofs.size();
441 :
442 : // Form edge projection matrix
443 46928377 : for (decltype(dofs_size) geomi = 0, freei = 0; geomi != dofs_size; ++geomi)
444 : {
445 41414640 : auto i = is_volume ? geomi : _side_dofs[geomi];
446 :
447 : // fixed DoFs aren't test functions
448 41414640 : if (_dof_is_fixed[i])
449 34222196 : continue;
450 60789786 : for (decltype(dofs_size) geomj = 0, freej = 0; geomj != dofs_size; ++geomj)
451 : {
452 53597342 : auto j = is_volume ? geomj : _side_dofs[geomj];
453 53597342 : if (_dof_is_fixed[j])
454 34227572 : _Fe(freei) -= (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp] * _Ue(j);
455 : else
456 19369770 : _Ke(freei, freej) += (*_phi)[i][_qp] * (*_phi)[j][_qp] * (*_JxW)[_qp];
457 53597342 : if (_cont == C_ONE)
458 : {
459 0 : if (_dof_is_fixed[j])
460 0 : _Fe(freei) -= dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp] * _Ue(j);
461 : else
462 0 : _Ke(freei, freej) += dotHelper((*_dphi)[i][_qp], (*_dphi)[j][_qp]) * (*_JxW)[_qp];
463 : }
464 53597342 : if (!_dof_is_fixed[j])
465 19369770 : freej++;
466 : }
467 7192444 : _Fe(freei) += (*_phi)[i][_qp] * fineval * (*_JxW)[_qp];
468 7192444 : if (_cont == C_ONE)
469 0 : _Fe(freei) += dotHelper(finegrad, (*_dphi)[i][_qp]) * (*_JxW)[_qp];
470 7192444 : freei++;
471 : }
472 : }
473 1851244 : }
474 :
475 : template <typename T>
476 : void
477 1826376 : InitialConditionTempl<T>::choleskySolve(const bool is_volume)
478 : {
479 1826376 : _Ke.resize(_free_dofs, _free_dofs);
480 1826376 : _Ke.zero();
481 1826376 : _Fe.resize(_free_dofs);
482 1826376 : _Fe.zero();
483 :
484 1826376 : choleskyAssembly(is_volume);
485 :
486 : // The new edge coefficients
487 1826372 : DenseVector<DataType> U(_free_dofs);
488 :
489 1826372 : _Ke.cholesky_solve(_Fe, U);
490 :
491 : // Transfer new edge solutions to element
492 3842382 : for (unsigned int i = 0; i != _free_dofs; ++i)
493 : {
494 2016010 : auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
495 2016010 : DataType & ui = _Ue(the_dof);
496 : libmesh_assert(std::abs(ui) < TOLERANCE || std::abs(ui - U(i)) < TOLERANCE);
497 2016010 : ui = U(i);
498 2016010 : _dof_is_fixed[the_dof] = true;
499 : }
500 1826372 : }
501 :
502 : template <>
503 : void
504 24872 : InitialConditionTempl<RealEigenVector>::choleskySolve(const bool is_volume)
505 : {
506 24872 : _Ke.resize(_free_dofs, _free_dofs);
507 24872 : _Ke.zero();
508 24872 : _Fe.resize(_free_dofs);
509 57616 : for (unsigned int i = 0; i < _free_dofs; ++i)
510 32744 : _Fe(i).setZero(_var.count());
511 :
512 24872 : choleskyAssembly(is_volume);
513 :
514 : // The new edge coefficients
515 24872 : DenseVector<DataType> U = _Fe;
516 :
517 76320 : for (unsigned int i = 0; i < _var.count(); ++i)
518 : {
519 51448 : DenseVector<Real> v(_free_dofs), x(_free_dofs);
520 120416 : for (unsigned int j = 0; j < _free_dofs; ++j)
521 68968 : v(j) = _Fe(j)(i);
522 :
523 51448 : _Ke.cholesky_solve(v, x);
524 :
525 120416 : for (unsigned int j = 0; j < _free_dofs; ++j)
526 68968 : U(j)(i) = x(j);
527 51448 : }
528 :
529 : // Transfer new edge solutions to element
530 57616 : for (unsigned int i = 0; i != _free_dofs; ++i)
531 : {
532 32744 : auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]];
533 32744 : DataType & ui = _Ue(the_dof);
534 : libmesh_assert(ui.matrix().norm() < TOLERANCE || (ui - U(i)).matrix().norm() < TOLERANCE);
535 32744 : ui = U(i);
536 32744 : _dof_is_fixed[the_dof] = true;
537 : }
538 24872 : }
539 :
540 : template <typename T>
541 : void
542 80 : InitialConditionTempl<T>::computeNodal(const Point & p)
543 : {
544 80 : _var.reinitNode();
545 80 : _var.computeNodalValues(); // has to call this to resize the internal array
546 80 : auto return_value = value(p);
547 :
548 80 : _var.setNodalValue(return_value); // update variable data, which is referenced by others, so the
549 : // value is up-to-date
550 :
551 : // We are done, so update the solution vector
552 80 : _var.insert(_var.sys().solution());
553 80 : }
554 :
555 : template class InitialConditionTempl<Real>;
556 : template class InitialConditionTempl<RealVectorValue>;
557 : template class InitialConditionTempl<RealEigenVector>;
|
