Skip to content

Commit c29a353

Browse files
hughcarshughcars
committed
Add optimized libCEED output functionals
- Add libCEED-backed surface/output functional infrastructure and qfunctions. - Route supported surface, port, and mode-overlap reductions through the new path. - Cover the functional and port semantics with focused regression tests.
1 parent f84b43d commit c29a353

34 files changed

Lines changed: 8124 additions & 48 deletions

‎palace/fem/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,13 +9,17 @@ target_sources(${LIB_TARGET_NAME}
99
PRIVATE
1010
${CMAKE_CURRENT_SOURCE_DIR}/bilinearform.cpp
1111
${CMAKE_CURRENT_SOURCE_DIR}/coefficient.cpp
12+
${CMAKE_CURRENT_SOURCE_DIR}/ceed_group_operator.cpp
1213
${CMAKE_CURRENT_SOURCE_DIR}/errorindicator.cpp
14+
${CMAKE_CURRENT_SOURCE_DIR}/facenbrexchange.cpp
1315
${CMAKE_CURRENT_SOURCE_DIR}/fespace.cpp
1416
${CMAKE_CURRENT_SOURCE_DIR}/gridfunction.cpp
1517
${CMAKE_CURRENT_SOURCE_DIR}/integrator.cpp
1618
${CMAKE_CURRENT_SOURCE_DIR}/interpolator.cpp
1719
${CMAKE_CURRENT_SOURCE_DIR}/lumpedelement.cpp
1820
${CMAKE_CURRENT_SOURCE_DIR}/mesh.cpp
21+
${CMAKE_CURRENT_SOURCE_DIR}/output_functionals.cpp
22+
${CMAKE_CURRENT_SOURCE_DIR}/postprocessing_backend.cpp
1923
${CMAKE_CURRENT_SOURCE_DIR}/integ/curlcurl.cpp
2024
${CMAKE_CURRENT_SOURCE_DIR}/integ/curlcurlmass.cpp
2125
${CMAKE_CURRENT_SOURCE_DIR}/integ/diffusion.cpp
@@ -30,6 +34,7 @@ target_sources(${LIB_TARGET_NAME}
3034
${CMAKE_CURRENT_SOURCE_DIR}/libceed/basis.cpp
3135
${CMAKE_CURRENT_SOURCE_DIR}/libceed/ceed.cpp
3236
${CMAKE_CURRENT_SOURCE_DIR}/libceed/coefficient.cpp
37+
${CMAKE_CURRENT_SOURCE_DIR}/libceed/functional.cpp
3338
${CMAKE_CURRENT_SOURCE_DIR}/libceed/integrator.cpp
3439
${CMAKE_CURRENT_SOURCE_DIR}/libceed/operator.cpp
3540
${CMAKE_CURRENT_SOURCE_DIR}/libceed/restriction.cpp

‎palace/fem/ceed_group_operator.cpp‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2+
// SPDX-License-Identifier: Apache-2.0
3+
4+
#include "fem/ceed_group_operator.hpp"
5+
6+
#include <mfem.hpp>
7+
8+
namespace palace
9+
{
10+
11+
namespace fem
12+
{
13+
14+
void DestroyGroupOperators(std::vector<CeedGroupOperator> &groups)
15+
{
16+
for (auto &group : groups)
17+
{
18+
for (auto &[field_vec, source] : group.field_vec_sources)
19+
{
20+
(void)source;
21+
PalaceCeedCall(group.ceed, CeedVectorDestroy(&field_vec));
22+
}
23+
group.field_vec_sources.clear();
24+
if (group.out_vec)
25+
{
26+
PalaceCeedCall(group.ceed, CeedVectorDestroy(&group.out_vec));
27+
group.out_size = 0;
28+
}
29+
if (group.ctx)
30+
{
31+
PalaceCeedCall(group.ceed, CeedQFunctionContextDestroy(&group.ctx));
32+
}
33+
if (group.op)
34+
{
35+
PalaceCeedCall(group.ceed, CeedOperatorDestroy(&group.op));
36+
}
37+
group.field_sources.clear();
38+
}
39+
groups.clear();
40+
}
41+
42+
void ApplyAddGroupOperators(const std::vector<CeedGroupOperator> &groups,
43+
const std::array<const Vector *, 4> &srcs, const Vector &out,
44+
const Vector *imported)
45+
{
46+
for (const auto &group : groups)
47+
{
48+
if (group.field_vec_sources.size() != group.field_sources.size())
49+
{
50+
group.field_vec_sources.clear();
51+
group.field_vec_sources.reserve(group.field_sources.size());
52+
for (const auto &[name, source] : group.field_sources)
53+
{
54+
CeedOperatorField field;
55+
CeedVector field_vec;
56+
PalaceCeedCall(group.ceed,
57+
CeedOperatorGetFieldByName(group.op, name.c_str(), &field));
58+
PalaceCeedCall(group.ceed, CeedOperatorFieldGetVector(field, &field_vec));
59+
group.field_vec_sources.emplace_back(field_vec, source);
60+
}
61+
}
62+
for (auto &[field_vec, source] : group.field_vec_sources)
63+
{
64+
// Source index 4 selects an optional imported vector, used by surface reductions
65+
// and boundary point fields for face-neighbor field values. The operator's
66+
// restriction slices and transposes the shared vector to the per-element layout.
67+
const Vector *sv = (source < 4) ? srcs[source] : imported;
68+
MFEM_ASSERT(sv, "Missing source vector for libCEED field input!");
69+
ceed::InitCeedVector(*sv, group.ceed, &field_vec, false);
70+
}
71+
CeedMemType out_mem;
72+
PalaceCeedCall(group.ceed, CeedGetPreferredMemType(group.ceed, &out_mem));
73+
if (!mfem::Device::Allows(mfem::Backend::DEVICE_MASK) && out_mem == CEED_MEM_DEVICE)
74+
{
75+
out_mem = CEED_MEM_HOST;
76+
}
77+
auto *out_data = const_cast<Vector &>(out).ReadWrite(out_mem == CEED_MEM_DEVICE);
78+
const CeedSize out_size = out.Size();
79+
if (!group.out_vec || group.out_size != out_size)
80+
{
81+
if (group.out_vec)
82+
{
83+
PalaceCeedCall(group.ceed, CeedVectorDestroy(&group.out_vec));
84+
}
85+
PalaceCeedCall(group.ceed, CeedVectorCreate(group.ceed, out_size, &group.out_vec));
86+
group.out_size = out_size;
87+
}
88+
PalaceCeedCall(group.ceed,
89+
CeedVectorSetArray(group.out_vec, out_mem, CEED_USE_POINTER, out_data));
90+
PalaceCeedCall(group.ceed, CeedOperatorApplyAdd(group.op, CEED_VECTOR_NONE,
91+
group.out_vec, CEED_REQUEST_IMMEDIATE));
92+
PalaceCeedCall(group.ceed, CeedVectorTakeArray(group.out_vec, out_mem, nullptr));
93+
}
94+
}
95+
96+
} // namespace fem
97+
98+
} // namespace palace

‎palace/fem/ceed_group_operator.hpp‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
// Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2+
// SPDX-License-Identifier: Apache-2.0
3+
4+
#ifndef PALACE_FEM_CEED_GROUP_OPERATOR_HPP
5+
#define PALACE_FEM_CEED_GROUP_OPERATOR_HPP
6+
7+
#include <array>
8+
#include <string>
9+
#include <utility>
10+
#include <vector>
11+
#include "fem/libceed/ceed.hpp"
12+
#include "linalg/vector.hpp"
13+
14+
namespace palace
15+
{
16+
17+
namespace fem
18+
{
19+
20+
// An assembled libCEED operator over one group of elements, with the named passive field
21+
// inputs re-pointed at caller data (by source vector index) on each evaluation.
22+
struct CeedGroupOperator
23+
{
24+
Ceed ceed = nullptr;
25+
CeedOperator op = nullptr;
26+
std::vector<std::pair<std::string, int>> field_sources;
27+
// Optional retained QFunction context handle for in-place runtime updates (e.g.
28+
// far-field frequency) without reassembly; nullptr if the operator has no context or
29+
// the context is not updated. Owned by the group (destroyed with it).
30+
CeedQFunctionContext ctx = nullptr;
31+
// Cached passive field vectors for field_sources, populated on first apply to avoid
32+
// repeated string lookups in libCEED during ParaView point-field output.
33+
mutable std::vector<std::pair<CeedVector, int>> field_vec_sources;
34+
// Reusable output vector wrapper. The pointed-to MFEM Vector data is supplied at
35+
// apply time, but the libCEED vector object itself can be retained across repeated
36+
// postprocessing evaluations instead of being created/destroyed for every group apply.
37+
mutable CeedVector out_vec = nullptr;
38+
mutable CeedSize out_size = 0;
39+
};
40+
41+
// Re-point the passive field inputs of each group operator at the given source vectors
42+
// and accumulate into the output vector with CeedOperatorApplyAdd. A field source index
43+
// of 4 (out of the srcs range) selects the optional imported vector instead, used to
44+
// feed face-neighbor (ghost) field values exchanged for two-sided interior boundaries on
45+
// parallel interfaces.
46+
void ApplyAddGroupOperators(const std::vector<CeedGroupOperator> &groups,
47+
const std::array<const Vector *, 4> &srcs, const Vector &out,
48+
const Vector *imported = nullptr);
49+
50+
// Destroy the libCEED references owned by a group-operator vector and clear it. The Ceed
51+
// context itself is borrowed and is not destroyed here.
52+
void DestroyGroupOperators(std::vector<CeedGroupOperator> &groups);
53+
54+
} // namespace fem
55+
56+
} // namespace palace
57+
58+
#endif // PALACE_FEM_CEED_GROUP_OPERATOR_HPP

‎palace/fem/coefficient.hpp‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -562,18 +562,20 @@ template <>
562562
inline double EnergyDensityCoefficient<EnergyDensityType::MAGNETIC>::GetLocalEnergyDensity(
563563
mfem::ElementTransformation &T) const
564564
{
565-
if (T.GetSpaceDim() == 2)
565+
if (T.GetSpaceDim() == 2 && U.Real().VectorDim() == 1)
566566
{
567-
// In 2D, B is the scalar out-of-plane component Bz (L2). The magnetic energy density
568-
// is 1/2 μ⁻¹_zz |Bz|², matching the integrated domain energy.
567+
// In 2D, B is scalar-valued and represents the out-of-plane flux density B_z.
568+
// Magnetic energy density is therefore 1/2 * mu^{-1}_{zz} * |B_z|^2, matching
569+
// the integrated domain energy.
570+
const double invmu_zz = mat_op.GetInvPermeabilityZZ(T.Attribute);
569571
double Bz = U.Real().GetValue(T, T.GetIntPoint());
570-
double dot = Bz * Bz;
572+
double dot = invmu_zz * Bz * Bz;
571573
if (U.HasImag())
572574
{
573575
Bz = U.Imag().GetValue(T, T.GetIntPoint());
574-
dot += Bz * Bz;
576+
dot += invmu_zz * Bz * Bz;
575577
}
576-
return 0.5 * mat_op.GetInvPermeabilityZZ(T.Attribute) * dot * scaling;
578+
return 0.5 * dot * scaling;
577579
}
578580
double V_data[3];
579581
mfem::Vector V(V_data, T.GetSpaceDim());

0 commit comments

Comments
 (0)