MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
spde_solver.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details
11
12#include <algorithm>
13#include <ctime>
14
16#include "spde_solver.hpp"
17
18namespace mfem
19{
20namespace spde
21{
22
23// Helper function that determines if output should be printed to the console.
24// The output is printed if the rank is 0 and if the print level is greater than
25// 0. The rank is retrieved via the fespace.
26bool PrintOutput(const ParFiniteElementSpace *fespace_ptr, int print_level)
27{
28 return (fespace_ptr->GetMyRank() == 0 && print_level > 0);
29}
30
31void Boundary::PrintInfo(std::ostream &os) const
32{
33 os << "\n<Boundary Info>\n"
34 << " Boundary Conditions:\n";
35 for (const auto &it : boundary_attributes)
36 {
37 os << " Boundary " << it.first << ": ";
38 switch (it.second)
39 {
41 os << "Neumann";
42 break;
44 os << "Dirichlet";
45 break;
47 os << "Robin, coefficient: " << robin_coefficient;
48 break;
50 os << "Periodic";
51 break;
52 default:
53 os << "Undefined";
54 break;
55 }
56 os << "\n";
57 }
58 bool first_print_statement = true;
59 // If the map is not empty
60 if (!dirichlet_coefficients.empty())
61 {
62 os << " Inhomogeneous Dirichlet defined on ";
63 for (const auto &it : dirichlet_coefficients)
64 {
65 if (!first_print_statement)
66 {
67 os << ", ";
68 }
69 else
70 {
71 first_print_statement = false;
72 }
73 os << it.first << "(=" << it.second << ")";
74 }
75 os << "\n";
76 }
77 os << "<Boundary Info>\n\n";
78}
79
81{
82 // Verify that all defined boundaries are actually defined on the
83 // mesh, i.e. if the keys of boundary attributes appear in the boundary
84 // attributes of the mesh.
85 mfem::out << "\n<Boundary Verify>\n";
86 const Array<int> boundary(mesh.bdr_attributes);
87 for (const auto &it : boundary_attributes)
88 {
89 if (boundary.Find(it.first) == -1)
90 {
91 MFEM_ABORT(" Boundary "
92 << it.first
93 << " is not defined on the mesh but in Boundary class."
94 << "Exiting...")
95 }
96 }
97
98 /// Verify if all boundary attributes appear as keys in the
99 /// boundary attributes, if not let the user know that we use Neumann by
100 /// default.
101 std::vector<int> boundary_attributes_keys;
102 for (int i = 0; i < boundary.Size(); i++)
103 {
104 if (boundary_attributes.find(boundary[i]) == boundary_attributes.end())
105 {
106 boundary_attributes_keys.push_back(boundary[i]);
107 }
108 }
109 if (!boundary_attributes_keys.empty())
110 {
111 mfem::out << " Boundaries (";
112 for (const auto &it : boundary_attributes_keys)
113 {
114 mfem::out << it << ", ";
115 }
116 mfem::out << ") are defined on the mesh but not in the";
117 mfem::out << " boundary attributes (Use Neumann).";
118 }
119
120 /// Check if any periodic boundary is registered
121 for (const auto &it : boundary_attributes)
122 {
123 if (it.second == BoundaryType::kPeriodic)
124 {
125 MFEM_ABORT(" Periodic boundaries must be defined on the mesh"
126 << ", not in Boundaries. Exiting...")
127 }
128 }
129
130 mfem::out << "\n<Boundary Verify>\n\n";
131}
132
134{
135 const ParFiniteElementSpace &fes = *solution.ParFESpace();
136 const ParMesh &pmesh = *fes.GetParMesh();
137
138 if (PrintOutput(&fes, 1))
139 {
140 mfem::out << "<Boundary::ComputeBoundaryError>"
141 << "\n";
142 mfem::out << " GetVDim: " << fes.GetVDim() << "\n";
143 }
144
145 real_t alpha{0.0};
146 real_t beta{1.0};
147 real_t gamma{0.0};
148
149 // Index i needs to be incremented by one to map to the boundary attributes
150 // in the mesh.
151 for (int i = 0; i < pmesh.bdr_attributes.Max(); i++)
152 {
153 real_t error{0};
154 real_t avg{0};
155 Array<int> bdr(pmesh.bdr_attributes.Max());
156 bdr = 0;
157 bdr[i] = 1;
158
160 avg = IntegrateBC(solution, bdr, alpha, beta, gamma, error);
161 if (PrintOutput(&fes, 1))
162 {
163 mfem::out << "->Boundary " << i + 1 << "\n";
164 mfem::out << " Alpha : " << alpha << "\n";
165 mfem::out << " Beta : " << beta << "\n";
166 mfem::out << " Gamma : " << gamma << "\n";
167 mfem::out << " Average : " << avg << "\n";
168 mfem::out << " Error : " << error << "\n\n";
169 }
170 }
171
172 if (PrintOutput(&fes, 1))
173 {
174 mfem::out << "<Boundary::ComputeBoundaryError>" << std::endl;
175 }
176}
177
179 real_t &gamma)
180{
181 // Check if i is a key in boundary_attributes
182 if (boundary_attributes.find(i) != boundary_attributes.end())
183 {
184 switch (boundary_attributes[i])
185 {
187 alpha = 1.0;
188 beta = 0.0;
189 gamma = 0.0;
190 break;
192 alpha = 0.0;
193 beta = 1.0;
195 {
196 gamma = dirichlet_coefficients[i];
197 }
198 else
199 {
200 gamma = 0.0;
201 }
202 break;
204 alpha = 1.0;
206 gamma = 0.0;
207 break;
208 default:
209 alpha = 1.0;
210 beta = 0.0;
211 gamma = 0.0;
212 break;
213 }
214 }
215 else
216 {
217 // If i is not a key in boundary_attributes, it corresponds to Neumann.
218 alpha = 1.0;
219 beta = 0.0;
220 gamma = 0.0;
221 }
222}
223
225 BoundaryType type)
226{
227 boundary_attributes[boundary] = type;
228}
229
231 real_t coefficient)
232{
234 dirichlet_coefficients[boundary] = coefficient;
235}
236
238{
239 robin_coefficient = coefficient;
240}
241
243 real_t alpha, real_t beta, real_t gamma, real_t &glb_err)
244{
245 real_t loc_vals[3];
246 real_t &nrm = loc_vals[0];
247 real_t &avg = loc_vals[1];
248 real_t &error = loc_vals[2];
249
250 nrm = 0.0;
251 avg = 0.0;
252 error = 0.0;
253
254 const bool a_is_zero = alpha == 0.0;
255 const bool b_is_zero = beta == 0.0;
256
257 const ParFiniteElementSpace &fes = *x.ParFESpace();
258 MFEM_ASSERT(fes.GetVDim() == 1, "");
259 ParMesh &mesh = *fes.GetParMesh();
260 Vector shape;
261 Vector loc_dofs;
262 Vector w_nor;
263 DenseMatrix dshape;
264 Array<int> dof_ids;
265 for (int i = 0; i < mesh.GetNBE(); i++)
266 {
267 if (bdr[mesh.GetBdrAttribute(i) - 1] == 0)
268 {
269 continue;
270 }
271
273 if (FTr == nullptr)
274 {
275 continue;
276 }
277
278 const FiniteElement &fe = *fes.GetFE(FTr->Elem1No);
279 MFEM_ASSERT(fe.GetMapType() == FiniteElement::VALUE, "");
280 const int int_order = 2 * fe.GetOrder() + 3;
281 const IntegrationRule &ir = IntRules.Get(FTr->FaceGeom, int_order);
282
283 fes.GetElementDofs(FTr->Elem1No, dof_ids);
284 x.GetSubVector(dof_ids, loc_dofs);
285 if (!a_is_zero)
286 {
287 const int sdim = FTr->Face->GetSpaceDim();
288 w_nor.SetSize(sdim);
289 dshape.SetSize(fe.GetDof(), sdim);
290 }
291 if (!b_is_zero)
292 {
293 shape.SetSize(fe.GetDof());
294 }
295 for (int j = 0; j < ir.GetNPoints(); j++)
296 {
297 const IntegrationPoint &ip = ir.IntPoint(j);
299 FTr->Loc1.Transform(ip, eip);
300 FTr->Face->SetIntPoint(&ip);
301 real_t face_weight = FTr->Face->Weight();
302 real_t val = 0.0;
303 if (!a_is_zero)
304 {
305 FTr->Elem1->SetIntPoint(&eip);
306 fe.CalcPhysDShape(*FTr->Elem1, dshape);
307 CalcOrtho(FTr->Face->Jacobian(), w_nor);
308 val += alpha * dshape.InnerProduct(w_nor, loc_dofs) / face_weight;
309 }
310 if (!b_is_zero)
311 {
312 fe.CalcShape(eip, shape);
313 val += beta * (shape * loc_dofs);
314 }
315
316 // Measure the length of the boundary
317 nrm += ip.weight * face_weight;
318
319 // Integrate alpha * n.Grad(x) + beta * x
320 avg += val * ip.weight * face_weight;
321
322 // Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
323 val -= gamma;
324 error += (val * val) * ip.weight * face_weight;
325 }
326 }
327
328 real_t glb_vals[3];
329 MPI_Allreduce(loc_vals, glb_vals, 3, MPITypeMap<real_t>::mpi_type, MPI_SUM,
330 fes.GetComm());
331
332 real_t glb_nrm = glb_vals[0];
333 real_t glb_avg = glb_vals[1];
334 glb_err = glb_vals[2];
335
336 // Normalize by the length of the boundary
337 if (std::abs(glb_nrm) > 0.0)
338 {
339 glb_err /= glb_nrm;
340 glb_avg /= glb_nrm;
341 }
342
343 // Compute l2 norm of the error in the boundary condition (negative
344 // quadrature weights may produce negative 'error')
345 glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
346
347 // Return the average value of alpha * n.Grad(x) + beta * x
348 return glb_avg;
349}
350
352 ParFiniteElementSpace *fespace, real_t l1, real_t l2,
353 real_t l3, real_t e1, real_t e2, real_t e3)
354 : k_(fespace),
355 m_(fespace),
356 fespace_ptr_(fespace),
357 bc_(bc),
358 nu_(nu),
359 l1_(l1),
360 l2_(l2),
361 l3_(l3),
362 e1_(e1),
363 e2_(e2),
364 e3_(e3)
365{
366 if (PrintOutput(fespace_ptr_, print_level_))
367 {
368 mfem::out << "<SPDESolver> Initialize Solver .." << std::endl;
369 }
370 StopWatch sw;
371 sw.Start();
372
373 // Resize the marker arrays for the boundary conditions
374 // Number of boundary attributes in the mesh
375 int nbc{0};
376 const auto &bdr_attributes = fespace_ptr_->GetParMesh()->bdr_attributes;
377 if (bdr_attributes.Size() > 0)
378 {
379 // Assumes a contiguous range of boundary attributes (1, 2, 3, ...)
380 nbc = bdr_attributes.Max() - bdr_attributes.Min() + 1;
381 }
382 dbc_marker_.SetSize(nbc);
383 rbc_marker_.SetSize(nbc);
384 dbc_marker_ = 0;
385 rbc_marker_ = 0;
386
387 // Fill the marker arrays for the boundary conditions. We decrement the number
388 // it.first by one because the boundary attributes in the mesh start at 1 and
389 // the marker arrays start at 0.
390 for (const auto &it : bc_.boundary_attributes)
391 {
392 switch (it.second)
393 {
395 dbc_marker_[it.first - 1] = 1;
396 break;
398 rbc_marker_[it.first - 1] = 1;
399 break;
400 default:
401 break;
402 }
403 }
404
405 // Handle homogeneous Dirichlet boundary conditions
406 // Note: for non zero DBC we usually need to project the boundary onto the
407 // solution. This is not necessary in this case since the boundary is
408 // homogeneous. For inhomogeneous Dirichlet we consider a lifting scheme.
409 fespace_ptr_->GetEssentialTrueDofs(dbc_marker_, ess_tdof_list_);
410
411 // Compute the rational approximation coefficients.
412 int dim = fespace_ptr_->GetParMesh()->Dimension();
413 int space_dim = fespace_ptr_->GetParMesh()->SpaceDimension();
414 alpha_ = (nu_ + dim / 2.0) / 2.0; // fractional exponent
415 integer_order_of_exponent_ = static_cast<int>(std::floor(alpha_));
416 real_t exponent_to_approximate = alpha_ - integer_order_of_exponent_;
417
418 // Compute the rational approximation coefficients.
419 ComputeRationalCoefficients(exponent_to_approximate);
420
421 // Set the bilinear forms.
422
423 // Assemble stiffness matrix
424 auto diffusion_tensor =
425 ConstructMatrixCoefficient(l1_, l2_, l3_, e1_, e2_, e3_, nu_, space_dim);
426 MatrixConstantCoefficient diffusion_coefficient(diffusion_tensor);
427 k_.AddDomainIntegrator(new DiffusionIntegrator(diffusion_coefficient));
428 ConstantCoefficient robin_coefficient(bc_.robin_coefficient);
429 k_.AddBoundaryIntegrator(new MassIntegrator(robin_coefficient), rbc_marker_);
430 k_.Assemble(0);
431
432 // Assemble mass matrix
433 ConstantCoefficient one(1.0);
435 m_.Assemble(0);
436
437 // Form matrices for the linear system
438 Array<int> empty;
439 k_.FormSystemMatrix(empty, stiffness_);
440 m_.FormSystemMatrix(empty, mass_bc_);
441
442 // Get the restriction and prolongation matrix for transformations
443 restriction_matrix_ = fespace->GetRestrictionMatrix();
444 prolongation_matrix_ = fespace->GetProlongationMatrix();
445
446 // Resize the vectors B and X to the appropriate size
447 if (prolongation_matrix_)
448 {
449 B_.SetSize(prolongation_matrix_->Width());
450 }
451 else
452 {
453 mfem::err << "<SPDESolver> prolongation matrix is not defined" << std::endl;
454 }
455 if (restriction_matrix_)
456 {
457 X_.SetSize(restriction_matrix_->Height());
458 }
459 else
460 {
461 mfem::err << "<SPDESolver> restriction matrix is not defined" << std::endl;
462 }
463
464 sw.Stop();
465 if (PrintOutput(fespace_ptr_, print_level_))
466 {
467 mfem::out << "<SPDESolver::Timing> matrix assembly " << sw.RealTime()
468 << " [s]" << std::endl;
469 }
470}
471
473{
474 // ------------------------------------------------------------------------
475 // Solve the PDE (A)^N g = f, i.e. compute g = (A)^{-1}^N f iteratively.
476 // ------------------------------------------------------------------------
477
478 StopWatch sw;
479 sw.Start();
480
481 // Zero initialize x to avoid touching uninitialized memory
482 x = 0.0;
483
484 ParGridFunction helper_gf(fespace_ptr_);
485 helper_gf = 0.0;
486
487 if (integer_order_of_exponent_ > 0)
488 {
489 if (PrintOutput(fespace_ptr_, print_level_))
490 {
491 mfem::out << "<SPDESolver> Solving PDE (A)^" << integer_order_of_exponent_
492 << " u = f" << std::endl;
493 }
494 ActivateRepeatedSolve();
495 Solve(b, helper_gf, 1.0, 1.0, integer_order_of_exponent_);
496 if (integer_order_)
497 {
498 // If the exponent is an integer, we can directly add the solution to the
499 // final solution and return.
500 x += helper_gf;
501 if (!bc_.dirichlet_coefficients.empty())
502 {
503 LiftSolution(x);
504 }
505 return;
506 }
507 UpdateRHS(b);
508 DeactivateRepeatedSolve();
509 }
510
511 // ------------------------------------------------------------------------
512 // Solve the (remaining) fractional PDE by solving M integer order PDEs and
513 // adding up the solutions.
514 // ------------------------------------------------------------------------
515 if (!integer_order_)
516 {
517 // Iterate over all expansion coefficient that contribute to the
518 // solution.
519 for (int i = 0; i < coeffs_.Size(); i++)
520 {
521 if (PrintOutput(fespace_ptr_, print_level_))
522 {
523 mfem::out << "\n<SPDESolver> Solving PDE -Δ u + " << -poles_[i]
524 << " u = " << coeffs_[i] << " g " << std::endl;
525 }
526 helper_gf = 0.0;
527 Solve(b, helper_gf, 1.0 - poles_[i], coeffs_[i]);
528 x += helper_gf;
529 }
530 }
531
532 // Apply the inhomogeneous Dirichlet boundary conditions.
533 if (!bc_.dirichlet_coefficients.empty())
534 {
535 LiftSolution(x);
536 }
537
538 sw.Stop();
539 if (PrintOutput(fespace_ptr_, print_level_))
540 {
541 mfem::out << "<SPDESolver::Timing> all PCG solves " << sw.RealTime()
542 << " [s]" << std::endl;
543 }
544}
545
547{
548 delete b_wn;
549 integ =
550 new WhiteGaussianNoiseDomainLFIntegrator(fespace_ptr_->GetComm(), seed);
551 b_wn = new ParLinearForm(fespace_ptr_);
552 b_wn->AddDomainIntegrator(integ);
553}
554
556{
557 if (!b_wn)
558 {
559 MFEM_ABORT("Need to call SPDESolver::SetupRandomFieldGenerator(...) first");
560 }
561 // Create stochastic load
562 b_wn->Assemble();
564 nu_, l1_, l2_, l3_, fespace_ptr_->GetParMesh()->Dimension());
565 (*b_wn) *= normalization;
566
567 // Call back to solve to generate the random field
568 Solve(*b_wn, x);
569}
570
572 real_t l2, real_t l3,
573 int dim)
574{
575 // Computation considers squaring components, computing determinant, and
576 // squaring
577 real_t det = 0;
578 if (dim == 1)
579 {
580 det = l1;
581 }
582 else if (dim == 2)
583 {
584 det = l1 * l2;
585 }
586 else if (dim == 3)
587 {
588 det = l1 * l2 * l3;
589 }
590 const real_t gamma1 = tgamma(nu + static_cast<real_t>(dim) / 2.0);
591 const real_t gamma2 = tgamma(nu);
592 return sqrt(pow(2 * M_PI, dim / 2.0) * det * gamma1 /
593 (gamma2 * pow(nu, dim / 2.0)));
594}
595
597 real_t l3, real_t e1,
598 real_t e2, real_t e3,
599 real_t nu, int dim)
600{
601 if (dim == 3)
602 {
603 // Compute cosine and sine of the angles e1, e2, e3
604 const real_t c1 = cos(e1);
605 const real_t s1 = sin(e1);
606 const real_t c2 = cos(e2);
607 const real_t s2 = sin(e2);
608 const real_t c3 = cos(e3);
609 const real_t s3 = sin(e3);
610
611 // Fill the rotation matrix R with the Euler angles.
612 DenseMatrix R(3, 3);
613 R(0, 0) = c1 * c3 - c2 * s1 * s3;
614 R(0, 1) = -c1 * s3 - c2 * c3 * s1;
615 R(0, 2) = s1 * s2;
616 R(1, 0) = c3 * s1 + c1 * c2 * s3;
617 R(1, 1) = c1 * c2 * c3 - s1 * s3;
618 R(1, 2) = -c1 * s2;
619 R(2, 0) = s2 * s3;
620 R(2, 1) = c3 * s2;
621 R(2, 2) = c2;
622
623 // Multiply the rotation matrix R with the translation vector.
624 Vector l(3);
625 l(0) = std::pow(l1, 2);
626 l(1) = std::pow(l2, 2);
627 l(2) = std::pow(l3, 2);
628 l *= (1 / (2.0 * nu));
629
630 // Compute result = R^t diag(l) R
631 DenseMatrix res(3, 3);
632 R.Transpose();
633 MultADBt(R, l, R, res);
634 return res;
635 }
636 else if (dim == 2)
637 {
638 const real_t c1 = cos(e1);
639 const real_t s1 = sin(e1);
640 DenseMatrix Rt(2, 2);
641 Rt(0, 0) = c1;
642 Rt(0, 1) = s1;
643 Rt(1, 0) = -s1;
644 Rt(1, 1) = c1;
645 Vector l(2);
646 l(0) = std::pow(l1, 2);
647 l(1) = std::pow(l2, 2);
648 l *= (1 / (2.0 * nu));
649 DenseMatrix res(2, 2);
650 MultADAt(Rt,l,res);
651 return res;
652 }
653 else
654 {
655 DenseMatrix res(1, 1);
656 res(0, 0) = std::pow(l1, 2) / (2.0 * nu);
657 return res;
658 }
659}
660
662 real_t beta, int exponent)
663{
664 // Form system of equations. This is less general than
665 // BilinearForm::FormLinearSystem and kind of resembles the necessary subset
666 // of instructions that we need in this case.
667 if (prolongation_matrix_)
668 {
669 prolongation_matrix_->MultTranspose(b, B_);
670 }
671 else
672 {
673 B_ = b;
674 }
675 B_ *= beta;
676
677 if (!apply_lift_)
678 {
679 // Initialize X_ to zero. Important! Might contain nan/inf -> crash.
680 X_ = 0.0;
681 }
682 else
683 {
684 restriction_matrix_->Mult(x, X_);
685 }
686
687 HypreParMatrix *Op =
688 Add(1.0, stiffness_, alpha, mass_bc_); // construct Operator
689 HypreParMatrix *Ae = Op->EliminateRowsCols(ess_tdof_list_);
690 Op->EliminateBC(*Ae, ess_tdof_list_, X_, B_); // only for homogeneous BC
691
692 for (int i = 0; i < exponent; i++)
693 {
694 // Solve the linear system Op X_ = B_
695 SolveLinearSystem(Op);
696 k_.RecoverFEMSolution(X_, b, x);
697 if (repeated_solve_)
698 {
699 // Prepare for next iteration. X is a primal and B is a dual vector. B_
700 // must be updated to represent X_ in the next step. Instead of copying
701 // it, we must transform it appropriately.
702 GridFunctionCoefficient gfc(&x);
703 ParLinearForm previous_solution(fespace_ptr_);
704 previous_solution.AddDomainIntegrator(new DomainLFIntegrator(gfc));
705 previous_solution.Assemble();
706 prolongation_matrix_->MultTranspose(previous_solution, B_);
707 Op->EliminateBC(*Ae, ess_tdof_list_, X_, B_);
708 }
709 }
710 delete Ae;
711 delete Op;
712}
713
714void SPDESolver::LiftSolution(ParGridFunction &x)
715{
716 // Set lifting flag
717 apply_lift_ = true;
718
719 // Lifting of the solution takes care of inhomogeneous boundary conditions.
720 // See doi:10.1016/j.jcp.2019.109009; section 2.6
721 if (PrintOutput(fespace_ptr_, print_level_))
722 {
723 mfem::out << "\n<SPDESolver> Applying inhomogeneous DBC" << std::endl;
724 }
725
726 // Define temporary grid function for lifting.
727 ParGridFunction helper_gf(fespace_ptr_);
728 helper_gf = 0.0;
729
730 // Project the boundary conditions onto the solution space.
731 for (const auto &bc : bc_.dirichlet_coefficients)
732 {
733 Array<int> marker(fespace_ptr_->GetParMesh()->bdr_attributes.Max());
734 marker = 0;
735 marker[bc.first - 1] = 1;
736 ConstantCoefficient cc(bc.second);
737 helper_gf.ProjectBdrCoefficient(cc, marker);
738 }
739
740 // Create linear form for the right hand side.
741 ParLinearForm b(fespace_ptr_);
742 ConstantCoefficient zero(0.0);
743 b.AddDomainIntegrator(new DomainLFIntegrator(zero));
744 b.Assemble();
745
746 // Solve the PDE for the lifting.
747 Solve(b, helper_gf, 1.0, 1.0);
748
749 // Add the lifting to the solution.
750 x += helper_gf;
751
752 // Reset the lifting flag.
753 apply_lift_ = false;
754}
755
756void SPDESolver::UpdateRHS(ParLinearForm &b) const
757{
758 if (!repeated_solve_)
759 {
760 // This function is only relevant for repeated solves.
761 return;
762 }
763 if (restriction_matrix_)
764 {
765 // This effectively writes the solution of the previous iteration X_ to the
766 // linear form b. Note that at the end of solve we update B_ = Mass * X_.
767 restriction_matrix_->MultTranspose(B_, b);
768 }
769 else
770 {
771 b = B_;
772 }
773}
774
775void SPDESolver::SolveLinearSystem(const HypreParMatrix *Op)
776{
777 HypreBoomerAMG prec(*Op);
778 prec.SetPrintLevel(-1);
779 CGSolver cg(fespace_ptr_->GetComm());
780 cg.SetRelTol(1e-12);
781 cg.SetMaxIter(2000);
782 cg.SetPrintLevel(3);
783 cg.SetPreconditioner(prec);
784 cg.SetOperator(*Op);
785 cg.SetPrintLevel(std::max(0, print_level_ - 1));
786 cg.Mult(B_, X_);
787}
788
789void SPDESolver::ComputeRationalCoefficients(real_t exponent)
790{
791 if (abs(exponent) > 1e-12)
792 {
793 if (PrintOutput(fespace_ptr_, print_level_))
794 {
795 mfem::out << "<SPDESolver> Approximating the fractional exponent "
796 << exponent << std::endl;
797 }
798 ComputePartialFractionApproximation(exponent, coeffs_, poles_);
799
800 // If the example is build without LAPACK, the exponent
801 // might be modified by the function call above.
802 alpha_ = exponent + integer_order_of_exponent_;
803 }
804 else
805 {
806 integer_order_ = true;
807 if (PrintOutput(fespace_ptr_, print_level_))
808 {
809 mfem::out << "<SPDESolver> Treating integer order PDE." << std::endl;
810 }
811 }
812}
813
814SPDESolver::~SPDESolver() { delete b_wn; }
815
816} // namespace spde
817} // namespace mfem
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
Definition array.hpp:828
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator. Assumes ownership of bfi.
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Transpose()
(*this) = (*this)^t
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
Definition densemat.cpp:383
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:119
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:484
ElementTransformation * Elem1
Definition eltrans.hpp:525
ElementTransformation * Face
Definition eltrans.hpp:526
IntegrationPointTransformation Loc1
Definition eltrans.hpp:527
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Abstract class for all finite elements.
Definition fe_base.hpp:239
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:333
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition fe_base.cpp:192
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:355
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition eltrans.cpp:543
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
A matrix coefficient that is constant in space and time.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1339
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1099
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition operator.hpp:93
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
const Operator * GetProlongationMatrix() const override
The returned Operator is owned by the FiniteElementSpace.
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
Definition pfespace.cpp:468
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
Definition pfespace.hpp:389
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
const FiniteElement * GetFE(int i) const override
Definition pfespace.cpp:534
Class for parallel grid function.
Definition pgridfunc.hpp:33
ParFiniteElementSpace * ParFESpace() const
Class for parallel linear form.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Timing object.
Definition tic_toc.hpp:36
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
void Solve(ParLinearForm &b, ParGridFunction &x)
static DenseMatrix ConstructMatrixCoefficient(real_t l1, real_t l2, real_t l3, real_t e1, real_t e2, real_t e3, real_t nu, int dim)
static real_t ConstructNormalizationCoefficient(real_t nu, real_t l1, real_t l2, real_t l3, int dim)
void GenerateRandomField(ParGridFunction &x)
void SetupRandomFieldGenerator(int seed=0)
SPDESolver(real_t nu, const Boundary &bc, ParFiniteElementSpace *fespace, real_t l1=0.1, real_t l2=0.1, real_t l3=0.1, real_t e1=0.0, real_t e2=0.0, real_t e3=0.0)
Vector beta
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
void ComputePartialFractionApproximation(real_t &alpha, Array< real_t > &coeffs, Array< real_t > &poles, real_t lmax=1000., real_t tol=1e-10, int npoints=1000, int max_order=100)
Definition ex33.hpp:295
real_t b
Definition lissajous.cpp:42
real_t IntegrateBC(const ParGridFunction &x, const Array< int > &bdr, real_t alpha, real_t beta, real_t gamma, real_t &glb_err)
bool PrintOutput(const ParFiniteElementSpace *fespace_ptr, int print_level)
void CalcOrtho(const DenseMatrix &J, Vector &n)
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition globals.hpp:71
float real_t
Definition config.hpp:43
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
Helper struct to convert a C++ type to an MPI type.
std::map< int, BoundaryType > boundary_attributes
Map to assign homogeneous boundary conditions to defined boundary types.
std::map< int, real_t > dirichlet_coefficients
Coefficient for inhomogeneous Dirichlet boundary conditions.
void PrintInfo(std::ostream &os=mfem::out) const
Print the information specifying the boundary conditions.
void ComputeBoundaryError(const ParGridFunction &solution)
void AddInhomogeneousDirichletBoundaryCondition(int boundary, real_t coefficient)
Add a inhomogeneous Dirichlet boundary condition to the boundary.
void SetRobinCoefficient(real_t coefficient)
Set the robin coefficient for the boundary.
void UpdateIntegrationCoefficients(int i, real_t &alpha, real_t &beta, real_t &gamma)
void VerifyDefinedBoundaries(const Mesh &mesh) const
void AddHomogeneousBoundaryCondition(int boundary, BoundaryType type)
Add a homogeneous boundary condition to the boundary.