Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "errorestimator.hpp"
5 :
6 : #include <limits>
7 : #include "fem/bilinearform.hpp"
8 : #include "fem/integrator.hpp"
9 : #include "fem/libceed/ceed.hpp"
10 : #include "fem/libceed/coefficient.hpp"
11 : #include "fem/libceed/integrator.hpp"
12 : #include "linalg/amg.hpp"
13 : #include "linalg/densematrix.hpp"
14 : #include "linalg/gmg.hpp"
15 : #include "linalg/iterative.hpp"
16 : #include "linalg/jacobi.hpp"
17 : #include "linalg/rap.hpp"
18 : #include "models/materialoperator.hpp"
19 : #include "utils/communication.hpp"
20 : #include "utils/diagnostic.hpp"
21 : #include "utils/omp.hpp"
22 : #include "utils/timer.hpp"
23 :
24 : PalacePragmaDiagnosticPush
25 : PalacePragmaDiagnosticDisableUnused
26 :
27 : #include "fem/qfunctions/hcurlhdiv_error_qf.h"
28 :
29 : PalacePragmaDiagnosticPop
30 :
31 : namespace palace
32 : {
33 :
34 : namespace
35 : {
36 :
37 : template <typename OperType>
38 : auto BuildLevelParOperator(std::unique_ptr<Operator> &&a,
39 : const FiniteElementSpace &trial_fespace,
40 : const FiniteElementSpace &test_fespace);
41 :
42 : template <>
43 : auto BuildLevelParOperator<Operator>(std::unique_ptr<Operator> &&a,
44 : const FiniteElementSpace &trial_fespace,
45 : const FiniteElementSpace &test_fespace)
46 : {
47 0 : return std::make_unique<ParOperator>(std::move(a), trial_fespace, test_fespace, false);
48 : }
49 :
50 : template <>
51 : auto BuildLevelParOperator<ComplexOperator>(std::unique_ptr<Operator> &&a,
52 : const FiniteElementSpace &trial_fespace,
53 : const FiniteElementSpace &test_fespace)
54 : {
55 0 : return std::make_unique<ComplexParOperator>(std::move(a), nullptr, trial_fespace,
56 0 : test_fespace, false);
57 : }
58 :
59 : template <typename OperType>
60 : auto BuildLevelParOperator(std::unique_ptr<Operator> &&a, const FiniteElementSpace &fespace)
61 : {
62 : return BuildLevelParOperator<OperType>(std::move(a), fespace, fespace);
63 : }
64 :
65 : template <typename OperType>
66 0 : auto ConfigureLinearSolver(const FiniteElementSpaceHierarchy &fespaces, double tol,
67 : int max_it, int print, bool use_mg)
68 : {
69 : // The system matrix for the projection is real, SPD and diagonally dominant.
70 0 : std::unique_ptr<Solver<OperType>> pc;
71 0 : if (!use_mg)
72 : {
73 : // Use eigenvalue estimate to compute optimal Jacobi damping parameter.
74 0 : pc = std::make_unique<JacobiSmoother<OperType>>(fespaces.GetFinestFESpace().GetComm(),
75 0 : 0.0);
76 : }
77 : else
78 : {
79 0 : auto amg = std::make_unique<BoomerAmgSolver>(1, 1, true, 0);
80 : amg->SetStrengthThresh(0.8); // More coarsening to save memory
81 0 : if (fespaces.GetNumLevels() > 1)
82 : {
83 0 : const int mg_smooth_order = 2; // Smooth order independent of FE space order
84 0 : pc = std::make_unique<GeometricMultigridSolver<OperType>>(
85 0 : fespaces.GetFinestFESpace().GetComm(),
86 0 : std::make_unique<MfemWrapperSolver<OperType>>(std::move(amg)),
87 0 : fespaces.GetProlongationOperators(), nullptr, 1, 1, mg_smooth_order, 1.0, 0.0,
88 0 : true);
89 : }
90 : else
91 : {
92 0 : pc = std::make_unique<MfemWrapperSolver<OperType>>(std::move(amg));
93 : }
94 : }
95 0 : auto pcg =
96 0 : std::make_unique<CgSolver<OperType>>(fespaces.GetFinestFESpace().GetComm(), print);
97 0 : pcg->SetInitialGuess(false);
98 : pcg->SetRelTol(tol);
99 : pcg->SetAbsTol(std::numeric_limits<double>::epsilon());
100 : pcg->SetMaxIter(max_it);
101 0 : return std::make_unique<BaseKspSolver<OperType>>(std::move(pcg), std::move(pc));
102 : }
103 :
104 : } // namespace
105 :
106 : template <typename VecType>
107 0 : FluxProjector<VecType>::FluxProjector(const MaterialPropertyCoefficient &coeff,
108 : const FiniteElementSpaceHierarchy &smooth_fespaces,
109 : const FiniteElementSpace &rhs_fespace, double tol,
110 0 : int max_it, int print, bool use_mg)
111 : {
112 0 : BlockTimer bt(Timer::CONSTRUCT_ESTIMATOR);
113 : const auto &smooth_fespace = smooth_fespaces.GetFinestFESpace();
114 : {
115 : constexpr bool skip_zeros = false;
116 : BilinearForm m(smooth_fespace);
117 0 : m.AddDomainIntegrator<VectorFEMassIntegrator>();
118 0 : if (!use_mg)
119 : {
120 0 : M = BuildLevelParOperator<OperType>(m.Assemble(skip_zeros), smooth_fespace);
121 : }
122 : else
123 : {
124 0 : auto m_vec = m.Assemble(smooth_fespaces, skip_zeros);
125 0 : auto M_mg =
126 0 : std::make_unique<BaseMultigridOperator<OperType>>(smooth_fespaces.GetNumLevels());
127 0 : for (std::size_t l = 0; l < smooth_fespaces.GetNumLevels(); l++)
128 : {
129 : const auto &fespace_l = smooth_fespaces.GetFESpaceAtLevel(l);
130 0 : M_mg->AddOperator(BuildLevelParOperator<OperType>(std::move(m_vec[l]), fespace_l));
131 : }
132 : M = std::move(M_mg);
133 0 : }
134 : }
135 : {
136 : // Flux operator is always partially assembled.
137 : BilinearForm flux(rhs_fespace, smooth_fespace);
138 0 : flux.AddDomainIntegrator<VectorFEMassIntegrator>(coeff);
139 0 : Flux = BuildLevelParOperator<OperType>(flux.PartialAssemble(), rhs_fespace,
140 : smooth_fespace);
141 : }
142 0 : ksp = ConfigureLinearSolver<OperType>(smooth_fespaces, tol, max_it, print, use_mg);
143 0 : ksp->SetOperators(*M, *M);
144 0 : rhs.SetSize(smooth_fespace.GetTrueVSize());
145 0 : rhs.UseDevice(true);
146 0 : }
147 :
148 : template <typename VecType>
149 0 : void FluxProjector<VecType>::Mult(const VecType &x, VecType &y) const
150 : {
151 0 : BlockTimer bt(Timer::SOLVE_ESTIMATOR);
152 : MFEM_ASSERT(x.Size() == Flux->Width() && y.Size() == rhs.Size(),
153 : "Invalid vector dimensions for FluxProjector::Mult!");
154 : // Mpi::Print(" Computing smooth flux recovery (projection) for error estimation\n");
155 0 : Flux->Mult(x, rhs);
156 0 : ksp->Mult(rhs, y);
157 0 : }
158 :
159 : namespace
160 : {
161 :
162 : template <typename VecType>
163 0 : Vector ComputeErrorEstimates(const VecType &F, VecType &F_gf, VecType &G, VecType &G_gf,
164 : const FiniteElementSpace &fespace,
165 : const FiniteElementSpace &smooth_fespace,
166 : const FluxProjector<VecType> &projector,
167 : const ceed::Operator &integ_op)
168 : {
169 : // Compute the projection of the discontinuous flux onto the smooth finite element space
170 : // (recovery) and populate the corresponding grid functions.
171 0 : BlockTimer bt(Timer::ESTIMATION);
172 0 : projector.Mult(F, G);
173 : if constexpr (std::is_same<VecType, ComplexVector>::value)
174 : {
175 0 : fespace.GetProlongationMatrix()->Mult(F.Real(), F_gf.Real());
176 0 : fespace.GetProlongationMatrix()->Mult(F.Imag(), F_gf.Imag());
177 0 : smooth_fespace.GetProlongationMatrix()->Mult(G.Real(), G_gf.Real());
178 0 : smooth_fespace.GetProlongationMatrix()->Mult(G.Imag(), G_gf.Imag());
179 : }
180 : else
181 : {
182 0 : fespace.GetProlongationMatrix()->Mult(F, F_gf);
183 0 : smooth_fespace.GetProlongationMatrix()->Mult(G, G_gf);
184 : }
185 :
186 : // Use libCEED operators to perform the error estimate integration over each element.
187 : const auto &mesh = fespace.GetMesh();
188 : Vector estimates(mesh.GetNE());
189 : estimates.UseDevice(true);
190 0 : estimates = 0.0;
191 0 : PalacePragmaOmp(parallel if (ceed::internal::NumCeeds() > 1))
192 : {
193 : Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
194 :
195 : // We need to update the state of the underlying libCEED vectors to indicate that the
196 : // data has changed. Each thread has it's own vector, referencing the same underlying
197 : // data.
198 : CeedVector F_gf_vec, G_gf_vec;
199 : {
200 : CeedInt nsub_ops;
201 : CeedOperator *sub_ops;
202 : PalaceCeedCall(
203 : ceed, CeedOperatorCompositeGetNumSub(integ_op[utils::GetThreadNum()], &nsub_ops));
204 : PalaceCeedCall(
205 : ceed, CeedOperatorCompositeGetSubList(integ_op[utils::GetThreadNum()], &sub_ops));
206 : MFEM_ASSERT(nsub_ops > 0, "Unexpected empty libCEED composite operator!");
207 : CeedOperatorField field;
208 : PalaceCeedCall(ceed, CeedOperatorGetFieldByName(sub_ops[0], "u_1", &field));
209 : PalaceCeedCall(ceed, CeedOperatorFieldGetVector(field, &F_gf_vec));
210 : PalaceCeedCall(ceed, CeedOperatorGetFieldByName(sub_ops[0], "u_2", &field));
211 : PalaceCeedCall(ceed, CeedOperatorFieldGetVector(field, &G_gf_vec));
212 : if constexpr (std::is_same<VecType, ComplexVector>::value)
213 : {
214 : ceed::InitCeedVector(F_gf.Real(), ceed, &F_gf_vec, false);
215 : ceed::InitCeedVector(G_gf.Real(), ceed, &G_gf_vec, false);
216 : }
217 : else
218 : {
219 : ceed::InitCeedVector(F_gf, ceed, &F_gf_vec, false);
220 : ceed::InitCeedVector(G_gf, ceed, &G_gf_vec, false);
221 : }
222 : }
223 :
224 : // Each thread writes to non-overlapping entries of the estimates vector.
225 : CeedVector estimates_vec;
226 : ceed::InitCeedVector(estimates, ceed, &estimates_vec);
227 :
228 : // Do the integration (both input vectors are passive). For the complex case, add sum of
229 : // squares of real and imaginary parts to the estimates before square root.
230 : PalaceCeedCall(ceed,
231 : CeedOperatorApplyAdd(integ_op[utils::GetThreadNum()], CEED_VECTOR_NONE,
232 : estimates_vec, CEED_REQUEST_IMMEDIATE));
233 : if constexpr (std::is_same<VecType, ComplexVector>::value)
234 : {
235 : ceed::InitCeedVector(F_gf.Imag(), ceed, &F_gf_vec, false);
236 : ceed::InitCeedVector(G_gf.Imag(), ceed, &G_gf_vec, false);
237 : PalaceCeedCall(ceed,
238 : CeedOperatorApplyAdd(integ_op[utils::GetThreadNum()], CEED_VECTOR_NONE,
239 : estimates_vec, CEED_REQUEST_IMMEDIATE));
240 : }
241 :
242 : // Cleanup.
243 : PalaceCeedCall(ceed, CeedVectorDestroy(&estimates_vec));
244 : }
245 :
246 0 : return estimates;
247 0 : }
248 :
249 : } // namespace
250 :
251 : template <typename VecType>
252 0 : GradFluxErrorEstimator<VecType>::GradFluxErrorEstimator(
253 : const MaterialOperator &mat_op, FiniteElementSpace &nd_fespace,
254 : FiniteElementSpaceHierarchy &rt_fespaces, double tol, int max_it, int print,
255 : bool use_mg)
256 0 : : nd_fespace(nd_fespace), rt_fespace(rt_fespaces.GetFinestFESpace()),
257 0 : projector(MaterialPropertyCoefficient(mat_op.GetAttributeToMaterial(),
258 : mat_op.GetPermittivityReal()),
259 : rt_fespaces, nd_fespace, tol, max_it, print, use_mg),
260 0 : integ_op(nd_fespace.GetMesh().GetNE(), nd_fespace.GetVSize()),
261 0 : E_gf(nd_fespace.GetVSize()), D(rt_fespace.GetTrueVSize()), D_gf(rt_fespace.GetVSize())
262 : {
263 0 : E_gf.UseDevice(true);
264 0 : D.UseDevice(true);
265 0 : D_gf.UseDevice(true);
266 :
267 : // Construct the libCEED operator used for integrating the element-wise error. The
268 : // discontinuous flux is ε E = ε ∇V.
269 : const auto &mesh = nd_fespace.GetMesh();
270 0 : PalacePragmaOmp(parallel if (ceed::internal::NumCeeds() > 1))
271 : {
272 : Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
273 : for (const auto &[geom, data] : mesh.GetCeedGeomFactorData(ceed))
274 : {
275 : // Only integrate over domain elements (not on the boundary).
276 : if (mfem::Geometry::Dimension[geom] < mesh.Dimension())
277 : {
278 : continue;
279 : }
280 :
281 : // Create libCEED vector wrappers for use with libCEED operators.
282 : CeedVector E_gf_vec, D_gf_vec;
283 : if constexpr (std::is_same<VecType, ComplexVector>::value)
284 : {
285 : ceed::InitCeedVector(E_gf.Real(), ceed, &E_gf_vec);
286 : ceed::InitCeedVector(D_gf.Real(), ceed, &D_gf_vec);
287 : }
288 : else
289 : {
290 : ceed::InitCeedVector(E_gf, ceed, &E_gf_vec);
291 : ceed::InitCeedVector(D_gf, ceed, &D_gf_vec);
292 : }
293 :
294 : // Construct mesh element restriction for elements of this element geometry type.
295 : CeedElemRestriction mesh_elem_restr;
296 : PalaceCeedCall(ceed, CeedElemRestrictionCreate(
297 : ceed, static_cast<CeedInt>(data.indices.size()), 1, 1,
298 : mesh.GetNE(), mesh.GetNE(), CEED_MEM_HOST, CEED_USE_POINTER,
299 : data.indices.data(), &mesh_elem_restr));
300 :
301 : // Element restriction and basis objects for inputs.
302 : CeedElemRestriction nd_restr =
303 : nd_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
304 : CeedElemRestriction rt_restr =
305 : rt_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
306 : CeedBasis nd_basis = nd_fespace.GetCeedBasis(ceed, geom);
307 : CeedBasis rt_basis = rt_fespace.GetCeedBasis(ceed, geom);
308 :
309 : // Construct coefficient for discontinuous flux, then smooth flux.
310 : auto mat_sqrtepsilon = linalg::MatrixSqrt(mat_op.GetPermittivityReal());
311 : auto mat_invsqrtepsilon = linalg::MatrixPow(mat_op.GetPermittivityReal(), -0.5);
312 : MaterialPropertyCoefficient sqrtepsilon_func(mat_op.GetAttributeToMaterial(),
313 : mat_sqrtepsilon);
314 : MaterialPropertyCoefficient invsqrtepsilon_func(mat_op.GetAttributeToMaterial(),
315 : mat_invsqrtepsilon);
316 : auto ctx =
317 : ceed::PopulateCoefficientContext(mesh.SpaceDimension(), &sqrtepsilon_func,
318 : mesh.SpaceDimension(), &invsqrtepsilon_func);
319 :
320 : // Assemble the libCEED operator. Inputs: E (for discontinuous flux), then smooth
321 : // flux.
322 : ceed::CeedQFunctionInfo info;
323 : info.assemble_q_data = false;
324 : switch (10 * mesh.SpaceDimension() + mesh.Dimension())
325 : {
326 : case 22:
327 : info.apply_qf = f_apply_hcurlhdiv_error_22;
328 : info.apply_qf_path = PalaceQFunctionRelativePath(f_apply_hcurlhdiv_error_22_loc);
329 : break;
330 : case 33:
331 : info.apply_qf = f_apply_hcurlhdiv_error_33;
332 : info.apply_qf_path = PalaceQFunctionRelativePath(f_apply_hcurlhdiv_error_33_loc);
333 : break;
334 : default:
335 : MFEM_ABORT("Invalid value of (dim, space_dim) = ("
336 : << mesh.Dimension() << ", " << mesh.SpaceDimension()
337 : << ") for GradFluxErrorEstimator!");
338 : }
339 : info.trial_ops = info.test_ops = ceed::EvalMode::Interp;
340 :
341 : CeedOperator sub_op;
342 : ceed::AssembleCeedElementErrorIntegrator(
343 : info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed, E_gf_vec,
344 : D_gf_vec, nd_restr, rt_restr, nd_basis, rt_basis, mesh_elem_restr, data.geom_data,
345 : data.geom_data_restr, &sub_op);
346 : integ_op.AddSubOperator(sub_op); // Sub-operator owned by ceed::Operator
347 :
348 : // Element restriction and passive input vectors are owned by the operator.
349 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&mesh_elem_restr));
350 : PalaceCeedCall(ceed, CeedVectorDestroy(&E_gf_vec));
351 : PalaceCeedCall(ceed, CeedVectorDestroy(&D_gf_vec));
352 : }
353 : }
354 :
355 : // Finalize the operator (call CeedOperatorCheckReady).
356 0 : integ_op.Finalize();
357 0 : }
358 :
359 : template <typename VecType>
360 0 : void GradFluxErrorEstimator<VecType>::AddErrorIndicator(const VecType &E, double Et,
361 : ErrorIndicator &indicator) const
362 : {
363 0 : auto estimates =
364 0 : ComputeErrorEstimates(E, E_gf, D, D_gf, nd_fespace, rt_fespace, projector, integ_op);
365 0 : linalg::Sqrt(estimates, (Et > 0.0) ? 0.5 / Et : 1.0); // Correct factor of 1/2 in energy
366 0 : indicator.AddIndicator(estimates);
367 0 : }
368 :
369 : template <typename VecType>
370 0 : CurlFluxErrorEstimator<VecType>::CurlFluxErrorEstimator(
371 : const MaterialOperator &mat_op, FiniteElementSpace &rt_fespace,
372 : FiniteElementSpaceHierarchy &nd_fespaces, double tol, int max_it, int print,
373 : bool use_mg)
374 0 : : rt_fespace(rt_fespace), nd_fespace(nd_fespaces.GetFinestFESpace()),
375 0 : projector(MaterialPropertyCoefficient(mat_op.GetAttributeToMaterial(),
376 : mat_op.GetInvPermeability()),
377 : nd_fespaces, rt_fespace, tol, max_it, print, use_mg),
378 0 : integ_op(nd_fespace.GetMesh().GetNE(), rt_fespace.GetVSize()),
379 0 : B_gf(rt_fespace.GetVSize()), H(nd_fespace.GetTrueVSize()), H_gf(nd_fespace.GetVSize())
380 : {
381 0 : B_gf.UseDevice(true);
382 0 : H.UseDevice(true);
383 0 : H_gf.UseDevice(true);
384 :
385 : // Construct the libCEED operator used for integrating the element-wise error. The
386 : // discontinuous flux is μ⁻¹ B ≃ μ⁻¹ ∇ × E.
387 : const auto &mesh = rt_fespace.GetMesh();
388 0 : PalacePragmaOmp(parallel if (ceed::internal::NumCeeds() > 1))
389 : {
390 : Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
391 : for (const auto &[geom, data] : mesh.GetCeedGeomFactorData(ceed))
392 : {
393 : // Only integrate over domain elements (not on the boundary).
394 : if (mfem::Geometry::Dimension[geom] < mesh.Dimension())
395 : {
396 : continue;
397 : }
398 :
399 : // Create libCEED vector wrappers for use with libCEED operators.
400 : CeedVector B_gf_vec, H_gf_vec;
401 : if constexpr (std::is_same<VecType, ComplexVector>::value)
402 : {
403 : ceed::InitCeedVector(B_gf.Real(), ceed, &B_gf_vec);
404 : ceed::InitCeedVector(H_gf.Real(), ceed, &H_gf_vec);
405 : }
406 : else
407 : {
408 : ceed::InitCeedVector(B_gf, ceed, &B_gf_vec);
409 : ceed::InitCeedVector(H_gf, ceed, &H_gf_vec);
410 : }
411 :
412 : // Construct mesh element restriction for elements of this element geometry type.
413 : CeedElemRestriction mesh_elem_restr;
414 : PalaceCeedCall(ceed, CeedElemRestrictionCreate(
415 : ceed, static_cast<CeedInt>(data.indices.size()), 1, 1,
416 : mesh.GetNE(), mesh.GetNE(), CEED_MEM_HOST, CEED_USE_POINTER,
417 : data.indices.data(), &mesh_elem_restr));
418 :
419 : // Element restriction and basis objects for inputs.
420 : CeedElemRestriction rt_restr =
421 : rt_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
422 : CeedElemRestriction nd_restr =
423 : nd_fespace.GetCeedElemRestriction(ceed, geom, data.indices);
424 : CeedBasis rt_basis = rt_fespace.GetCeedBasis(ceed, geom);
425 : CeedBasis nd_basis = nd_fespace.GetCeedBasis(ceed, geom);
426 :
427 : // Construct coefficient for discontinuous flux, then smooth flux.
428 : auto mat_invsqrtmu = linalg::MatrixSqrt(mat_op.GetInvPermeability());
429 : auto mat_sqrtmu = linalg::MatrixPow(mat_op.GetInvPermeability(), -0.5);
430 : MaterialPropertyCoefficient invsqrtmu_func(mat_op.GetAttributeToMaterial(),
431 : mat_invsqrtmu);
432 : MaterialPropertyCoefficient sqrtmu_func(mat_op.GetAttributeToMaterial(), mat_sqrtmu);
433 : auto ctx = ceed::PopulateCoefficientContext(mesh.SpaceDimension(), &invsqrtmu_func,
434 : mesh.SpaceDimension(), &sqrtmu_func);
435 :
436 : // Assemble the libCEED operator. Inputs: B (for discontinuous flux), then smooth
437 : // flux. Currently only supports 3D, since curl in 2D requires special treatment.
438 : ceed::CeedQFunctionInfo info;
439 : info.assemble_q_data = false;
440 : switch (10 * mesh.SpaceDimension() + mesh.Dimension())
441 : {
442 : case 33:
443 : info.apply_qf = f_apply_hdivhcurl_error_33;
444 : info.apply_qf_path = PalaceQFunctionRelativePath(f_apply_hdivhcurl_error_33_loc);
445 : break;
446 : default:
447 : MFEM_ABORT("Invalid value of (dim, space_dim) = ("
448 : << mesh.Dimension() << ", " << mesh.SpaceDimension()
449 : << ") for CurlFluxErrorEstimator!");
450 : }
451 : info.trial_ops = info.test_ops = ceed::EvalMode::Interp;
452 :
453 : CeedOperator sub_op;
454 : ceed::AssembleCeedElementErrorIntegrator(
455 : info, (void *)ctx.data(), ctx.size() * sizeof(CeedIntScalar), ceed, B_gf_vec,
456 : H_gf_vec, rt_restr, nd_restr, rt_basis, nd_basis, mesh_elem_restr, data.geom_data,
457 : data.geom_data_restr, &sub_op);
458 : integ_op.AddSubOperator(sub_op); // Sub-operator owned by ceed::Operator
459 :
460 : // Element restriction and passive input vectors are owned by the operator.
461 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&mesh_elem_restr));
462 : PalaceCeedCall(ceed, CeedVectorDestroy(&B_gf_vec));
463 : PalaceCeedCall(ceed, CeedVectorDestroy(&H_gf_vec));
464 : }
465 : }
466 :
467 : // Finalize the operator (call CeedOperatorCheckReady).
468 0 : integ_op.Finalize();
469 0 : }
470 :
471 : template <typename VecType>
472 0 : void CurlFluxErrorEstimator<VecType>::AddErrorIndicator(const VecType &B, double Et,
473 : ErrorIndicator &indicator) const
474 : {
475 0 : auto estimates =
476 0 : ComputeErrorEstimates(B, B_gf, H, H_gf, rt_fespace, nd_fespace, projector, integ_op);
477 0 : linalg::Sqrt(estimates, (Et > 0.0) ? 0.5 / Et : 1.0); // Correct factor of 1/2 in energy
478 0 : indicator.AddIndicator(estimates);
479 0 : }
480 :
481 : template <typename VecType>
482 0 : TimeDependentFluxErrorEstimator<VecType>::TimeDependentFluxErrorEstimator(
483 : const MaterialOperator &mat_op, FiniteElementSpaceHierarchy &nd_fespaces,
484 : FiniteElementSpaceHierarchy &rt_fespaces, double tol, int max_it, int print,
485 : bool use_mg)
486 0 : : grad_estimator(mat_op, nd_fespaces.GetFinestFESpace(), rt_fespaces, tol, max_it, print,
487 : use_mg),
488 0 : curl_estimator(mat_op, rt_fespaces.GetFinestFESpace(), nd_fespaces, tol, max_it, print,
489 : use_mg)
490 : {
491 0 : }
492 :
493 : template <typename VecType>
494 0 : void TimeDependentFluxErrorEstimator<VecType>::AddErrorIndicator(
495 : const VecType &E, const VecType &B, double Et, ErrorIndicator &indicator) const
496 : {
497 0 : auto grad_estimates =
498 0 : ComputeErrorEstimates(E, grad_estimator.E_gf, grad_estimator.D, grad_estimator.D_gf,
499 : grad_estimator.nd_fespace, grad_estimator.rt_fespace,
500 0 : grad_estimator.projector, grad_estimator.integ_op);
501 0 : auto curl_estimates =
502 0 : ComputeErrorEstimates(B, curl_estimator.B_gf, curl_estimator.H, curl_estimator.H_gf,
503 : curl_estimator.rt_fespace, curl_estimator.nd_fespace,
504 0 : curl_estimator.projector, curl_estimator.integ_op);
505 0 : grad_estimates += curl_estimates; // Sum of squares
506 0 : linalg::Sqrt(grad_estimates,
507 : (Et > 0.0) ? 0.5 / Et : 1.0); // Correct factor of 1/2 in energy
508 0 : indicator.AddIndicator(grad_estimates);
509 0 : }
510 :
511 : template class FluxProjector<Vector>;
512 : template class FluxProjector<ComplexVector>;
513 : template class GradFluxErrorEstimator<Vector>;
514 : template class GradFluxErrorEstimator<ComplexVector>;
515 : template class CurlFluxErrorEstimator<Vector>;
516 : template class CurlFluxErrorEstimator<ComplexVector>;
517 : template class TimeDependentFluxErrorEstimator<Vector>;
518 : template class TimeDependentFluxErrorEstimator<ComplexVector>;
519 :
520 : } // namespace palace
|
