Compare commits
1
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Julian Andrej
|
fafaef82b1 |
+1
-1
@@ -214,7 +214,7 @@ miniapps/adjoint/adjoint_advection_diffusion
|
||||
|
||||
miniapps/dfem/dfem-minimal-surface
|
||||
miniapps/dfem/dfem-minimal-surface-output
|
||||
miniapps/dfem/dfem-hyperelasticity
|
||||
miniapps/dfem/dfem-hyperelasticity_energy
|
||||
|
||||
miniapps/electromagnetics/volta
|
||||
miniapps/electromagnetics/tesla
|
||||
|
||||
@@ -16,10 +16,8 @@
|
||||
#include "kernels.hpp"
|
||||
#include "util.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
namespace ker = mfem::kernels::internal;
|
||||
|
||||
@@ -594,10 +592,10 @@ class DerivativeAssemble
|
||||
const std::array<DofToQuadMap, n_outputs> output_dtq_maps;
|
||||
const std::array<bool, n_inputs> input_is_dependent;
|
||||
const size_t trial_field_uf;
|
||||
/// Column space of every row block. GetDerivative differentiates w.r.t. one
|
||||
/// field, so there is exactly one trial space. Null when that field is not
|
||||
/// an FE space, in which case nothing can be assembled.
|
||||
const size_t test_field_uf;
|
||||
const ParFiniteElementSpace *test_fes;
|
||||
const ParFiniteElementSpace *trial_fes;
|
||||
const int test_vdim;
|
||||
/// Per-output row geometry of the quadrature point cache. DerivativeSetup
|
||||
/// lays that cache out over every output FieldOperator, so reading it needs
|
||||
/// all of them.
|
||||
@@ -605,19 +603,7 @@ class DerivativeAssemble
|
||||
const std::array<int, n_outputs> out_op_dim;
|
||||
const std::array<int, n_outputs> out_offsets;
|
||||
const int output_size_on_qp;
|
||||
/// Output field ids, in order of first appearance among @a outputs. Each one
|
||||
/// is a row block of the derivative: see @ref compute_group_field_ids.
|
||||
const std::vector<int> group_field_ids;
|
||||
/// Output FieldOperator -> row block it contributes to.
|
||||
const std::array<int, n_outputs> out_group;
|
||||
/// Row block -> position in ctx.outfds (-1 if the field is not an output of
|
||||
/// the operator, which should not happen).
|
||||
const std::vector<int> group_outfd_idx;
|
||||
/// Row block -> test space, null for spaces without a basis.
|
||||
const std::vector<const ParFiniteElementSpace *> group_fes;
|
||||
const std::vector<int> group_test_vdim;
|
||||
const std::vector<int> group_num_test_dof;
|
||||
const std::vector<bool> group_assemblable;
|
||||
const int num_test_dof;
|
||||
const int trial_vdim;
|
||||
const int trial_op_dim;
|
||||
const int num_trial_dof;
|
||||
@@ -625,36 +611,7 @@ class DerivativeAssemble
|
||||
const int num_trial_dof_1d;
|
||||
const int total_trial_op_dim;
|
||||
mutable Vector inputs_trial_op_dim;
|
||||
/// One element matrix bank per row block. Outputs on the same test field
|
||||
/// share a bank and accumulate into it; outputs on different test fields are
|
||||
/// different row blocks and need banks of their own, because their element
|
||||
/// matrices have different shapes and are filled through different
|
||||
/// ElementRestrictions.
|
||||
mutable std::vector<Vector> group_Ae_mem;
|
||||
|
||||
/// @brief Distinct output field ids, in order of first appearance.
|
||||
///
|
||||
/// Outputs sharing a field id form one row block: they are summed into a
|
||||
/// single element matrix and produce a single assembled matrix. This is the
|
||||
/// multi-output case, e.g. Outputs<Value<U>, Gradient<U>> giving mass plus
|
||||
/// diffusion. Outputs on different field ids are separate row blocks of
|
||||
///
|
||||
/// dR/dU = [ dR_U/dU ; dR_Y/dU ]
|
||||
///
|
||||
/// and are assembled into one matrix each.
|
||||
static std::vector<int> compute_group_field_ids(const outputs_t &outs)
|
||||
{
|
||||
std::vector<int> ids;
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
const int fid = get<o>(outs).GetFieldId();
|
||||
if (std::find(ids.begin(), ids.end(), fid) == ids.end())
|
||||
{
|
||||
ids.push_back(fid);
|
||||
}
|
||||
});
|
||||
return ids;
|
||||
}
|
||||
mutable Vector Ae_mem;
|
||||
|
||||
public:
|
||||
DerivativeAssemble() = delete;
|
||||
@@ -701,18 +658,29 @@ public:
|
||||
ctx_in.ir)),
|
||||
input_is_dependent(compute_input_is_dependent(inputs, derivative_id)),
|
||||
trial_field_uf(find_union_field_index(ctx_in, derivative_id)),
|
||||
trial_fes(
|
||||
[&]() -> const ParFiniteElementSpace *
|
||||
test_field_uf(
|
||||
find_union_field_index(ctx_in, get<0>(outputs).GetFieldId())),
|
||||
test_fes(
|
||||
[&]
|
||||
{
|
||||
const auto *fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[test_field_uf].data);
|
||||
MFEM_ASSERT(fes != nullptr && *fes != nullptr,
|
||||
"LocalQFBackend: test space is not a ParFiniteElementSpace");
|
||||
return *fes;
|
||||
}()),
|
||||
trial_fes(
|
||||
[&]
|
||||
{
|
||||
// Not every field is an FE space: a derivative w.r.t. a ParameterSpace or
|
||||
// a QuadratureSpace has no basis to assemble columns against. That is
|
||||
// only an error if assembly is actually requested, so keep it null here
|
||||
// and report it in operator() instead of aborting registration.
|
||||
if (trial_field_uf >= ctx_in.unionfds.size()) { return nullptr; }
|
||||
const auto *fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[trial_field_uf].data);
|
||||
return fes ? *fes : nullptr;
|
||||
MFEM_ASSERT(fes != nullptr && *fes != nullptr,
|
||||
"LocalQFBackend: trial space is not a ParFiniteElementSpace");
|
||||
return *fes;
|
||||
}()),
|
||||
// All outputs are attached to the same test field, so vdim is common to
|
||||
// them; only the operator dimension differs, and that lives in out_op_dim.
|
||||
test_vdim(get<0>(outputs).vdim),
|
||||
out_vdim(get_vdim(outputs)),
|
||||
out_op_dim(compute_out_op_dim(outputs)),
|
||||
out_offsets(compute_out_offsets(out_vdim, out_op_dim)),
|
||||
@@ -723,96 +691,7 @@ public:
|
||||
for_constexpr<n_outputs>([&](auto o) { s += get<o>(outputs).size_on_qp; });
|
||||
return s;
|
||||
}()),
|
||||
group_field_ids(compute_group_field_ids(outputs_in)),
|
||||
out_group(
|
||||
[&]
|
||||
{
|
||||
std::array<int, n_outputs> g {};
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
const int fid = get<o>(outputs_in).GetFieldId();
|
||||
const auto &ids = group_field_ids;
|
||||
g[o] = static_cast<int>(std::find(ids.begin(), ids.end(), fid)
|
||||
- ids.begin());
|
||||
});
|
||||
return g;
|
||||
}()),
|
||||
group_outfd_idx(
|
||||
[&]
|
||||
{
|
||||
std::vector<int> idx(group_field_ids.size(), -1);
|
||||
for (size_t g = 0; g < group_field_ids.size(); g++)
|
||||
{
|
||||
for (size_t f = 0; f < ctx_in.outfds.size(); f++)
|
||||
{
|
||||
if (static_cast<int>(ctx_in.outfds[f].id) == group_field_ids[g])
|
||||
{
|
||||
idx[g] = static_cast<int>(f);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
return idx;
|
||||
}()),
|
||||
group_fes(
|
||||
[&]
|
||||
{
|
||||
std::vector<const ParFiniteElementSpace *> v(group_field_ids.size(),
|
||||
nullptr);
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
if (group_outfd_idx[g] < 0) { continue; }
|
||||
const auto *fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.outfds[group_outfd_idx[g]].data);
|
||||
v[g] = fes ? *fes : nullptr;
|
||||
}
|
||||
return v;
|
||||
}()),
|
||||
group_test_vdim(
|
||||
[&]
|
||||
{
|
||||
// Outputs in a group share a field, hence a vdim; only the operator
|
||||
// dimension differs between them, and that lives in out_op_dim.
|
||||
std::vector<int> v(group_field_ids.size(), 0);
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
v[out_group[o]] = get<o>(outputs_in).vdim;
|
||||
});
|
||||
return v;
|
||||
}()),
|
||||
group_num_test_dof(
|
||||
[&]
|
||||
{
|
||||
std::vector<int> v(group_field_ids.size(), 0);
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
if (group_fes[g] == nullptr) { continue; }
|
||||
v[g] = group_fes[g]->GetFE(0)->GetDof();
|
||||
}
|
||||
return v;
|
||||
}()),
|
||||
group_assemblable(
|
||||
[&]
|
||||
{
|
||||
// A row block can only be assembled when both its trial and test fields are
|
||||
// ParFiniteElementSpaces. Quadrature and parameter spaces have no element
|
||||
// basis or ElementRestriction through which to form a SparseMatrix.
|
||||
//
|
||||
// Identity outputs are also excluded for now as it has no supported mapping from
|
||||
// its pointwise quadrature rows into that FE row space.
|
||||
std::vector<bool> v(group_field_ids.size(), false);
|
||||
if (trial_fes == nullptr) { return v; }
|
||||
for (size_t g = 0; g < v.size(); g++) { v[g] = group_fes[g] != nullptr; }
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
using output_fop_t = std::decay_t<decltype(get<o>(outputs_in))>;
|
||||
if constexpr (is_identity_fop_v<output_fop_t>)
|
||||
{
|
||||
v[out_group[o]] = false;
|
||||
}
|
||||
});
|
||||
return v;
|
||||
}()),
|
||||
num_test_dof(test_fes->GetFE(0)->GetDof()),
|
||||
trial_vdim(compute_trial_vdim(inputs, derivative_id)), trial_op_dim(
|
||||
[&]
|
||||
{
|
||||
@@ -826,7 +705,7 @@ public:
|
||||
});
|
||||
return top;
|
||||
}()),
|
||||
num_trial_dof(trial_fes ? trial_fes->GetFE(0)->GetDof() : 0),
|
||||
num_trial_dof(trial_fes->GetFE(0)->GetDof()),
|
||||
dim(ctx_in.mesh.Dimension()), ne(ctx_in.nentities),
|
||||
nq(ctx_in.ir.GetNPoints()), q1d(tensor_1d_size(nq, dim)),
|
||||
num_trial_dof_1d(tensor_1d_size(num_trial_dof, dim)), total_trial_op_dim(
|
||||
@@ -837,12 +716,14 @@ public:
|
||||
return compute_total_trial_op_dim(
|
||||
inputs, input_is_dependent, in_qp_sizes);
|
||||
}()),
|
||||
inputs_trial_op_dim(), group_Ae_mem()
|
||||
inputs_trial_op_dim(), Ae_mem()
|
||||
{
|
||||
MFEM_ASSERT(ctx.unionfds.size() == nfields,
|
||||
"LocalQFBackend: unionfds size mismatch");
|
||||
MFEM_ASSERT(trial_field_uf != SIZE_MAX,
|
||||
"DerivativeAssemble: trial field not found in unionfds");
|
||||
MFEM_ASSERT(test_field_uf != SIZE_MAX,
|
||||
"DerivativeAssemble: test field not found in unionfds");
|
||||
|
||||
MFEM_ASSERT(trial_vdim > 0,
|
||||
"LocalQFBackend: could not determine trial vdim");
|
||||
@@ -860,35 +741,35 @@ public:
|
||||
: 0;
|
||||
});
|
||||
|
||||
group_Ae_mem.resize(group_field_ids.size());
|
||||
for (size_t g = 0; g < group_Ae_mem.size(); g++)
|
||||
{
|
||||
if (!group_assemblable[g]) { continue; }
|
||||
const int elem_mat_size = group_num_test_dof[g] * group_test_vdim[g] *
|
||||
num_trial_dof * trial_vdim;
|
||||
group_Ae_mem[g].SetSize(elem_mat_size * ne,
|
||||
Device::GetDeviceMemoryType());
|
||||
group_Ae_mem[g].UseDevice(true);
|
||||
group_Ae_mem[g] = 0.0;
|
||||
}
|
||||
const int elem_mat_size =
|
||||
num_test_dof * test_vdim * num_trial_dof * trial_vdim;
|
||||
Ae_mem.SetSize(elem_mat_size * ne, Device::GetDeviceMemoryType());
|
||||
Ae_mem.UseDevice(true);
|
||||
Ae_mem = 0.0;
|
||||
}
|
||||
|
||||
/// @brief Assemble one SparseMatrix per assemblable output field.
|
||||
///
|
||||
/// @a A is indexed by position in ctx.outfds, so A[f] receives the row block
|
||||
/// belonging to output field f. Slots this integrator writes nothing
|
||||
/// assemblable to -- quadrature or parameter spaces, or fields it does not
|
||||
/// touch at all -- are left alone, which lets several integrators fill
|
||||
/// different row blocks of the same vector. Ownership of the matrices passes
|
||||
/// to the caller.
|
||||
void operator()(std::vector<SparseMatrix *> &A) const
|
||||
void operator()(SparseMatrix *&A) const
|
||||
{
|
||||
if (ctx.attr.Size() == 0) { return; }
|
||||
// Every output is contracted into one element matrix Ae, sized from the
|
||||
// test space of get<0>(outputs), and filled through a single test
|
||||
// ElementRestriction.
|
||||
//
|
||||
// WIP:
|
||||
// This takes care of single-field, multiple-outputs case.
|
||||
// For a multiple fields case, outputs on a second field would need a second
|
||||
// matrix -- the derivative then eould be a block column with one row block per
|
||||
// test space.
|
||||
//
|
||||
// For now we just add a check that all outputs are attached to the same test field, and abort if not.
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
MFEM_VERIFY(get<o>(outputs).GetFieldId() == get<0>(outputs).GetFieldId(),
|
||||
"DerivativeAssemble: every output FieldOperator has to be "
|
||||
"attached to the same test field; assembling outputs that "
|
||||
"span several fields is not supported");
|
||||
});
|
||||
|
||||
const bool any_assemblable =
|
||||
std::find(group_assemblable.begin(), group_assemblable.end(), true) !=
|
||||
group_assemblable.end();
|
||||
if (!any_assemblable) { return; }
|
||||
if (ctx.attr.Size() == 0) { return; }
|
||||
|
||||
if (!(use_sum_factorization && (dim == 2 || dim == 3)))
|
||||
{
|
||||
@@ -896,78 +777,43 @@ public:
|
||||
"for tensor-product 2D/3D elements only");
|
||||
}
|
||||
|
||||
MFEM_VERIFY(trial_fes != nullptr,
|
||||
"DerivativeAssemble: the differentiated field is not a "
|
||||
"ParFiniteElementSpace, so the columns of the derivative "
|
||||
"have no basis to be assembled against");
|
||||
DerivativeAssembleHO::Run(dim,
|
||||
q1d,
|
||||
ctx,
|
||||
qp_cache,
|
||||
Ae_mem,
|
||||
inputs,
|
||||
outputs,
|
||||
input_dtq_maps,
|
||||
output_dtq_maps,
|
||||
out_vdim,
|
||||
out_op_dim,
|
||||
out_offsets,
|
||||
output_size_on_qp,
|
||||
inputs_trial_op_dim,
|
||||
test_vdim,
|
||||
num_test_dof,
|
||||
num_trial_dof,
|
||||
num_trial_dof_1d,
|
||||
trial_vdim,
|
||||
total_trial_op_dim,
|
||||
nq,
|
||||
ne,
|
||||
q1d,
|
||||
dim);
|
||||
|
||||
A = new SparseMatrix;
|
||||
A->OverrideSize(test_fes->GetVSize(), trial_fes->GetVSize());
|
||||
|
||||
const auto *test_restr = dynamic_cast<const ElementRestriction *>(
|
||||
test_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC));
|
||||
const auto *trial_restr = dynamic_cast<const ElementRestriction *>(
|
||||
trial_fes->GetElementRestriction(
|
||||
ElementDofOrdering::LEXICOGRAPHIC));
|
||||
MFEM_VERIFY(trial_restr != nullptr,
|
||||
trial_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC));
|
||||
MFEM_VERIFY(test_restr != nullptr && trial_restr != nullptr,
|
||||
"DerivativeAssemble SparseMatrix assembly requires "
|
||||
"H1/conforming ElementRestriction spaces");
|
||||
|
||||
if (A.size() < ctx.outfds.size())
|
||||
{
|
||||
A.resize(ctx.outfds.size(), nullptr);
|
||||
}
|
||||
|
||||
for (size_t g = 0; g < group_Ae_mem.size(); g++)
|
||||
{
|
||||
if (!group_assemblable[g]) { continue; }
|
||||
|
||||
// The kernel accumulates into Ae, so the bank has to start clean on
|
||||
// every call: otherwise assembling twice (the next Newton step, say)
|
||||
// would double every entry.
|
||||
group_Ae_mem[g] = 0.0;
|
||||
|
||||
DerivativeAssembleHO::Run(dim,
|
||||
q1d,
|
||||
ctx,
|
||||
qp_cache,
|
||||
group_Ae_mem[g],
|
||||
inputs,
|
||||
outputs,
|
||||
input_dtq_maps,
|
||||
output_dtq_maps,
|
||||
out_vdim,
|
||||
out_op_dim,
|
||||
out_offsets,
|
||||
out_group,
|
||||
static_cast<int>(g),
|
||||
output_size_on_qp,
|
||||
inputs_trial_op_dim,
|
||||
group_test_vdim[g],
|
||||
group_num_test_dof[g],
|
||||
num_trial_dof,
|
||||
num_trial_dof_1d,
|
||||
trial_vdim,
|
||||
total_trial_op_dim,
|
||||
nq,
|
||||
ne,
|
||||
q1d,
|
||||
dim);
|
||||
|
||||
const ParFiniteElementSpace *test_fes = group_fes[g];
|
||||
const auto *test_restr = dynamic_cast<const ElementRestriction *>(
|
||||
test_fes->GetElementRestriction(
|
||||
ElementDofOrdering::LEXICOGRAPHIC));
|
||||
MFEM_VERIFY(test_restr != nullptr,
|
||||
"DerivativeAssemble SparseMatrix assembly requires "
|
||||
"H1/conforming ElementRestriction spaces");
|
||||
|
||||
const int f = group_outfd_idx[g];
|
||||
MFEM_VERIFY(A[f] == nullptr,
|
||||
"DerivativeAssemble: output field already carries an "
|
||||
"assembled matrix; two integrators contributing to the "
|
||||
"same row block cannot be assembled into one matrix");
|
||||
|
||||
auto *M = new SparseMatrix;
|
||||
M->OverrideSize(test_fes->GetVSize(), trial_fes->GetVSize());
|
||||
test_restr->FillSparseMatrix(group_Ae_mem[g], *M, *trial_restr);
|
||||
A[f] = M;
|
||||
}
|
||||
test_restr->FillSparseMatrix(Ae_mem, *A, *trial_restr);
|
||||
}
|
||||
|
||||
template<typename backend_t = LocalQFHOBackend<3>, int T_Q1D = 0>
|
||||
@@ -982,8 +828,6 @@ public:
|
||||
const std::array<int, n_outputs> &out_vdim,
|
||||
const std::array<int, n_outputs> &out_op_dim,
|
||||
const std::array<int, n_outputs> &out_offsets,
|
||||
const std::array<int, n_outputs> &out_group,
|
||||
const int group,
|
||||
const int output_size_on_qp,
|
||||
const Vector &inputs_trial_op_dim,
|
||||
const int test_vdim,
|
||||
@@ -1009,7 +853,6 @@ public:
|
||||
"DerivativeAssemble: nq exceeds backend quadrature capacity");
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
if (out_group[o] != group) { return; }
|
||||
MFEM_VERIFY(out_op_dim[o] <= DIM,
|
||||
"DerivativeAssemble: test_op_dim exceeds spatial DIM");
|
||||
});
|
||||
@@ -1053,25 +896,14 @@ public:
|
||||
|
||||
// Each output contributes its own rows of the cache, contracted
|
||||
// against its own test basis operation; map_quadrature_data_to_fields
|
||||
// accumulates, so the element matrix is the sum over the outputs that
|
||||
// belong to this row block. Outputs on other test fields are skipped
|
||||
// here and picked up by their own launch, which has its own Ae.
|
||||
//
|
||||
// A launch only ever happens for an assemblable row block, and such a
|
||||
// block contains no Identity output, so the is_identity_fop_v test
|
||||
// below never rejects anything at run time. It is there to keep
|
||||
// assemble_element_mat_sumfact from being instantiated for a fop it
|
||||
// cannot handle.
|
||||
// accumulates, so the element matrix is the sum over outputs for the
|
||||
// same field.
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
using output_fop_t = std::decay_t<decltype(get<o>(outputs))>;
|
||||
|
||||
if constexpr (!is_identity_fop_v<output_fop_t>)
|
||||
{
|
||||
// Uniform across the thread block, so returning early here
|
||||
// cannot desynchronise the barriers below.
|
||||
if (out_group[o] != group) { return; }
|
||||
|
||||
// The outputs share fhat_storage, so one has to be done with it
|
||||
// before the next zeroes it.
|
||||
MFEM_SYNC_THREAD;
|
||||
|
||||
@@ -15,18 +15,12 @@
|
||||
#include "kernels.hpp"
|
||||
#include "util.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <vector>
|
||||
|
||||
namespace mfem::future::LocalQFImpl
|
||||
{
|
||||
|
||||
// Assemble the diagonal of one row block of a cached Jacobian (tensor 2D/3D).
|
||||
//
|
||||
// The derivative is a block column, one row block per output field. Only a
|
||||
// block whose test space is the trial space is square, and only a square block
|
||||
// has a diagonal at all, so the row block is chosen per call and checked.
|
||||
// Assemble diagonal of cached Jacobian (square trial == test, tensor 2D/3D)
|
||||
|
||||
template<int derivative_id,
|
||||
typename qfunc_t,
|
||||
@@ -53,47 +47,22 @@ class DerivativeAssembleDiagonal
|
||||
const std::array<DofToQuadMap, n_outputs> output_dtq_maps;
|
||||
const std::array<bool, n_inputs> input_is_dependent;
|
||||
const size_t trial_field_uf;
|
||||
/// Column space of every row block; null when the differentiated field is
|
||||
/// not an FE space, in which case no block has a diagonal.
|
||||
const ParFiniteElementSpace *trial_fes;
|
||||
const size_t test_field_uf;
|
||||
const bool is_square;
|
||||
const int test_vdim;
|
||||
const std::array<int, n_outputs> out_vdim;
|
||||
const std::array<int, n_outputs> out_op_dim;
|
||||
const std::array<int, n_outputs> out_offsets;
|
||||
const int output_size_on_qp;
|
||||
/// Row blocks: the distinct output field ids, in order of first appearance.
|
||||
/// Outputs sharing a field id are summed into one diagonal, which is how
|
||||
/// Value<U> + Gradient<U> becomes mass plus diffusion.
|
||||
const std::vector<int> group_field_ids;
|
||||
const std::array<int, n_outputs> out_group;
|
||||
const std::vector<const ParFiniteElementSpace *> group_fes;
|
||||
const std::vector<int> group_test_vdim;
|
||||
const std::vector<int> group_num_test_dof;
|
||||
const std::vector<int> group_num_test_dof_1d;
|
||||
/// Whether a row block has a diagonal: its test space has to be an FE space
|
||||
/// and has to *be* the trial space, and no output on it may be an Identity.
|
||||
const std::vector<bool> group_has_diagonal;
|
||||
const int num_test_dof;
|
||||
const int num_test_dof_1d;
|
||||
const int trial_vdim;
|
||||
const int total_trial_op_dim;
|
||||
const int num_trial_dof_1d;
|
||||
const int residual_size_on_qp;
|
||||
const int dim, ne, nq, q1d;
|
||||
const std::array<int, n_inputs> inputs_trial_op_dim;
|
||||
mutable std::vector<Vector> group_Ye_mem;
|
||||
|
||||
/// Distinct output field ids, in order of first appearance.
|
||||
static std::vector<int> compute_group_field_ids(const outputs_t &outs)
|
||||
{
|
||||
std::vector<int> ids;
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
const int fid = get<o>(outs).GetFieldId();
|
||||
if (std::find(ids.begin(), ids.end(), fid) == ids.end())
|
||||
{
|
||||
ids.push_back(fid);
|
||||
}
|
||||
});
|
||||
return ids;
|
||||
}
|
||||
mutable Vector Ye_mem;
|
||||
|
||||
public:
|
||||
DerivativeAssembleDiagonal() = delete;
|
||||
@@ -140,14 +109,19 @@ public:
|
||||
ctx_in.ir)),
|
||||
input_is_dependent(compute_input_is_dependent(inputs, derivative_id)),
|
||||
trial_field_uf(find_union_field_index(ctx_in, derivative_id)),
|
||||
trial_fes(
|
||||
[&]() -> const ParFiniteElementSpace *
|
||||
test_field_uf(
|
||||
find_union_field_index(ctx_in, get<0>(outputs).GetFieldId())),
|
||||
is_square(
|
||||
[&]
|
||||
{
|
||||
if (trial_field_uf >= ctx_in.unionfds.size()) { return nullptr; }
|
||||
const auto *fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
const auto *test_fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[test_field_uf].data);
|
||||
const auto *trial_fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[trial_field_uf].data);
|
||||
return fes ? *fes : nullptr;
|
||||
return test_fes && trial_fes && *test_fes && *trial_fes &&
|
||||
(*test_fes == *trial_fes);
|
||||
}()),
|
||||
test_vdim(get<0>(outputs).vdim),
|
||||
out_vdim(get_vdim(outputs_in)),
|
||||
out_op_dim(compute_out_op_dim(outputs_in)),
|
||||
out_offsets(compute_out_offsets(out_vdim, out_op_dim)),
|
||||
@@ -158,96 +132,16 @@ public:
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{ s += get<o>(outputs_in).size_on_qp; });
|
||||
return s;
|
||||
}()),
|
||||
group_field_ids(compute_group_field_ids(outputs_in)),
|
||||
out_group(
|
||||
}()), num_test_dof(
|
||||
[&]
|
||||
{
|
||||
std::array<int, n_outputs> g {};
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
const int fid = get<o>(outputs_in).GetFieldId();
|
||||
const auto &ids = group_field_ids;
|
||||
g[o] = static_cast<int>(std::find(ids.begin(), ids.end(), fid)
|
||||
- ids.begin());
|
||||
});
|
||||
return g;
|
||||
}()),
|
||||
group_fes(
|
||||
[&]
|
||||
{
|
||||
std::vector<const ParFiniteElementSpace *> v(group_field_ids.size(),
|
||||
nullptr);
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
const size_t uf = find_union_field_index(ctx_in, group_field_ids[g]);
|
||||
if (uf >= ctx_in.unionfds.size()) { continue; }
|
||||
const auto *fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[uf].data);
|
||||
v[g] = fes ? *fes : nullptr;
|
||||
}
|
||||
return v;
|
||||
}()),
|
||||
group_test_vdim(
|
||||
[&]
|
||||
{
|
||||
std::vector<int> v(group_field_ids.size(), 0);
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
v[out_group[o]] = get<o>(outputs_in).vdim;
|
||||
});
|
||||
return v;
|
||||
}()),
|
||||
group_num_test_dof(
|
||||
[&]
|
||||
{
|
||||
std::vector<int> v(group_field_ids.size(), 0);
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
if (group_fes[g] == nullptr) { continue; }
|
||||
v[g] = group_fes[g]->GetFE(0)->GetDof();
|
||||
}
|
||||
return v;
|
||||
}()),
|
||||
group_num_test_dof_1d(
|
||||
[&]
|
||||
{
|
||||
std::vector<int> v(group_field_ids.size(), 0);
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
if (group_num_test_dof[g] > 0)
|
||||
{
|
||||
v[g] = tensor_1d_size(group_num_test_dof[g],
|
||||
ctx_in.mesh.Dimension());
|
||||
}
|
||||
}
|
||||
return v;
|
||||
}()),
|
||||
group_has_diagonal(
|
||||
[&]
|
||||
{
|
||||
// A diagonal needs row space == column space, so only a row block on the
|
||||
// trial space qualifies. Squareness alone cannot pick a block when
|
||||
// several output fields share that space, which is why the caller names
|
||||
// the row. Identity outputs are quadrature point data and are excluded
|
||||
// for the same reason as in DerivativeAssemble: they cannot be
|
||||
// contracted, and every output on a field lands in the same block.
|
||||
std::vector<bool> v(group_field_ids.size(), false);
|
||||
if (trial_fes == nullptr) { return v; }
|
||||
for (size_t g = 0; g < v.size(); g++)
|
||||
{
|
||||
v[g] = (group_fes[g] != nullptr) && (group_fes[g] == trial_fes);
|
||||
}
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
using output_fop_t = std::decay_t<decltype(get<o>(outputs_in))>;
|
||||
if constexpr (is_identity_fop_v<output_fop_t>)
|
||||
{
|
||||
v[out_group[o]] = false;
|
||||
}
|
||||
});
|
||||
return v;
|
||||
const auto *test_fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[test_field_uf].data);
|
||||
MFEM_ASSERT(test_fes != nullptr && *test_fes != nullptr,
|
||||
"LocalQFBackend: test space is not a ParFiniteElementSpace");
|
||||
return (*test_fes)->GetFE(0)->GetDof();
|
||||
}()),
|
||||
num_test_dof_1d(tensor_1d_size(num_test_dof, ctx_in.mesh.Dimension())),
|
||||
trial_vdim(compute_trial_vdim(inputs, derivative_id)), total_trial_op_dim(
|
||||
[&]
|
||||
{
|
||||
@@ -257,9 +151,15 @@ public:
|
||||
inputs, input_is_dependent, input_size_on_qp);
|
||||
}()),
|
||||
num_trial_dof_1d(
|
||||
trial_fes ? tensor_1d_size(trial_fes->GetFE(0)->GetDof(),
|
||||
ctx_in.mesh.Dimension())
|
||||
: 0),
|
||||
[&]
|
||||
{
|
||||
const auto *trial_fes = std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx_in.unionfds[trial_field_uf].data);
|
||||
MFEM_ASSERT(trial_fes != nullptr && *trial_fes != nullptr,
|
||||
"LocalQFBackend: trial space is not a ParFiniteElementSpace");
|
||||
const int num_trial_dof = (*trial_fes)->GetFE(0)->GetDof();
|
||||
return tensor_1d_size(num_trial_dof, ctx_in.mesh.Dimension());
|
||||
}()),
|
||||
residual_size_on_qp(output_size_on_qp * trial_vdim * total_trial_op_dim),
|
||||
dim(ctx_in.mesh.Dimension()), ne(ctx_in.nentities),
|
||||
nq(ctx_in.ir.GetNPoints()), q1d(tensor_1d_size(nq, dim)),
|
||||
@@ -275,62 +175,53 @@ public:
|
||||
});
|
||||
return itod;
|
||||
}()),
|
||||
group_Ye_mem()
|
||||
Ye_mem()
|
||||
{
|
||||
MFEM_ASSERT(ctx.unionfds.size() == nfields,
|
||||
"LocalQFBackend: unionfds size mismatch");
|
||||
MFEM_ASSERT(
|
||||
trial_field_uf != SIZE_MAX,
|
||||
"DerivativeAssembleDiagonal: trial field not found in unionfds");
|
||||
MFEM_ASSERT(
|
||||
test_field_uf != SIZE_MAX,
|
||||
"DerivativeAssembleDiagonal: test field not found in unionfds");
|
||||
MFEM_ASSERT(trial_vdim > 0,
|
||||
"LocalQFBackend: could not determine trial vdim");
|
||||
MFEM_ASSERT(total_trial_op_dim > 0,
|
||||
"LocalQFBackend: no dependent inputs found");
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
MFEM_ASSERT(out_vdim[o] == group_test_vdim[out_group[o]],
|
||||
"DerivativeAssembleDiagonal: outputs on one field must "
|
||||
"share its vdim");
|
||||
MFEM_CONTRACT_VAR(o);
|
||||
MFEM_ASSERT(out_vdim[o] == test_vdim,
|
||||
"DerivativeAssembleDiagonal: all outputs must share the "
|
||||
"test field vdim");
|
||||
});
|
||||
|
||||
group_Ye_mem.resize(group_field_ids.size());
|
||||
for (size_t g = 0; g < group_Ye_mem.size(); g++)
|
||||
if (is_square)
|
||||
{
|
||||
if (!group_has_diagonal[g]) { continue; }
|
||||
group_Ye_mem[g].SetSize(group_num_test_dof[g] * group_test_vdim[g] *
|
||||
ne);
|
||||
group_Ye_mem[g].UseDevice(true);
|
||||
Ye_mem.SetSize(num_test_dof * test_vdim * ne);
|
||||
Ye_mem.UseDevice(true);
|
||||
}
|
||||
}
|
||||
|
||||
/// Index of the row block for output field @a field_id, or -1.
|
||||
int FindGroup(int field_id) const
|
||||
{
|
||||
const auto &ids = group_field_ids;
|
||||
const auto it = std::find(ids.begin(), ids.end(), field_id);
|
||||
return (it == ids.end()) ? -1 : static_cast<int>(it - ids.begin());
|
||||
}
|
||||
|
||||
template<typename Backend>
|
||||
void run_kernels(const int g) const
|
||||
void run_kernels() const
|
||||
{
|
||||
Backend::Run(dim,
|
||||
q1d,
|
||||
ctx,
|
||||
qp_cache,
|
||||
group_Ye_mem[g],
|
||||
Ye_mem,
|
||||
inputs,
|
||||
outputs,
|
||||
output_dtq_maps,
|
||||
input_dtq_maps,
|
||||
out_group,
|
||||
g,
|
||||
group_test_vdim[g],
|
||||
test_vdim,
|
||||
out_op_dim,
|
||||
out_offsets,
|
||||
output_size_on_qp,
|
||||
group_num_test_dof[g],
|
||||
group_num_test_dof_1d[g],
|
||||
num_test_dof,
|
||||
num_test_dof_1d,
|
||||
trial_vdim,
|
||||
total_trial_op_dim,
|
||||
residual_size_on_qp,
|
||||
@@ -341,14 +232,9 @@ public:
|
||||
dim);
|
||||
}
|
||||
|
||||
/// Add this integrator's contribution to the diagonal of the row block of
|
||||
/// output field @a out_field_id. Adds nothing if the integrator writes no
|
||||
/// square, basis-backed block for that field; the caller is responsible for
|
||||
/// rejecting a row that no integrator can serve.
|
||||
void operator()(const int out_field_id, Vector &diag_e) const
|
||||
void operator()(Vector &diag_e) const
|
||||
{
|
||||
const int g = FindGroup(out_field_id);
|
||||
if (g < 0 || !group_has_diagonal[g]) { return; }
|
||||
if (!is_square) { return; }
|
||||
if (ctx.attr.Size() == 0) { return; }
|
||||
|
||||
if (!(use_sum_factorization && (dim == 2 || dim == 3)))
|
||||
@@ -356,28 +242,27 @@ public:
|
||||
MFEM_ABORT("DerivativeAssembleDiagonal optimized path is implemented "
|
||||
"for tensor-product 2D/3D elements only");
|
||||
}
|
||||
MFEM_VERIFY(group_num_test_dof_1d[g] == num_trial_dof_1d,
|
||||
MFEM_VERIFY(num_test_dof_1d == num_trial_dof_1d,
|
||||
"DerivativeAssembleDiagonal requires matching tensor dofs");
|
||||
const auto &limits = DeviceDofQuadLimits::Get();
|
||||
MFEM_VERIFY(group_num_test_dof_1d[g] <= limits.MAX_D1D, "");
|
||||
MFEM_VERIFY(q1d <= limits.MAX_Q1D, "");
|
||||
MFEM_VERIFY(num_test_dof_1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
|
||||
MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
|
||||
|
||||
group_Ye_mem[g] = 0.0;
|
||||
Ye_mem = 0.0;
|
||||
|
||||
if (q1d <= LocalQFLOBackendMQ1())
|
||||
{
|
||||
run_kernels<DerivativeAssembleDiagonalLO>(g);
|
||||
run_kernels<DerivativeAssembleDiagonalLO>();
|
||||
}
|
||||
else if (q1d <= LocalQFHOBackendMQ1())
|
||||
{
|
||||
run_kernels<DerivativeAssembleDiagonalHO>(g);
|
||||
run_kernels<DerivativeAssembleDiagonalHO>();
|
||||
}
|
||||
else
|
||||
{
|
||||
MFEM_ABORT("Unsupported quadrature order for LocalQF backend");
|
||||
}
|
||||
|
||||
diag_e += group_Ye_mem[g];
|
||||
diag_e += Ye_mem;
|
||||
}
|
||||
|
||||
template<typename backend_t = LocalQFLOBackend<3>, int T_Q1D = 0>
|
||||
@@ -389,8 +274,6 @@ public:
|
||||
const outputs_t &outputs,
|
||||
const std::array<DofToQuadMap, n_outputs> &output_dtq_maps,
|
||||
const std::array<DofToQuadMap, n_inputs> &input_dtq_maps,
|
||||
const std::array<int, n_outputs> &out_group,
|
||||
const int row_group,
|
||||
const int test_vdim,
|
||||
const std::array<int, n_outputs> &out_op_dim,
|
||||
const std::array<int, n_outputs> &out_offsets,
|
||||
@@ -456,92 +339,85 @@ public:
|
||||
}
|
||||
MFEM_SYNC_THREAD;
|
||||
|
||||
// Accumulate every output belonging to the requested row block.
|
||||
// This sums multiple contributions, such as Value<U> +
|
||||
// Gradient<U>, while skipping outputs on the other row blocks.
|
||||
// The row is a run time choice, so unlike the field id it cannot
|
||||
// gate the instantiation; is_identity_fop_v still does, since
|
||||
// eval_test has no meaning for quadrature point data.
|
||||
// Accumulate every (output o, test op k, dependent input s,
|
||||
// trial op m) block of the cached Jacobian into the diagonal via
|
||||
// the backend driver.
|
||||
for_constexpr<n_outputs>([&](auto o)
|
||||
{
|
||||
using test_fop_t = std::decay_t<decltype(get<o>(outputs))>;
|
||||
if constexpr (!is_identity_fop_v<test_fop_t>)
|
||||
const auto &out_dtq = output_dtq_maps[o];
|
||||
const int test_op_dim = out_op_dim[static_cast<int>(o)];
|
||||
|
||||
// Test-basis factor along a spatial axis
|
||||
const auto eval_test =
|
||||
[&](const int k, const int axis, const int q, const int d)
|
||||
{
|
||||
if (out_group[static_cast<int>(o)] != row_group) { return; }
|
||||
const auto &out_dtq = output_dtq_maps[o];
|
||||
const int test_op_dim = out_op_dim[static_cast<int>(o)];
|
||||
|
||||
// Test-basis factor along a spatial axis
|
||||
const auto eval_test =
|
||||
[&](const int k, const int axis, const int q, const int d)
|
||||
const auto &B = out_dtq.B;
|
||||
const auto &G = out_dtq.G;
|
||||
if constexpr (is_value_fop<test_fop_t>::value)
|
||||
{
|
||||
const auto &B = out_dtq.B;
|
||||
const auto &G = out_dtq.G;
|
||||
if constexpr (is_value_fop<test_fop_t>::value)
|
||||
{
|
||||
return (k == 0) ? B(q, 0, d) : 0.0;
|
||||
}
|
||||
else if constexpr (is_gradient_fop<test_fop_t>::value)
|
||||
{
|
||||
return (k == axis) ? G(q, 0, d) : B(q, 0, d);
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
for (int k = 0; k < test_op_dim; k++)
|
||||
{
|
||||
const int row = out_offsets[static_cast<int>(o)] +
|
||||
vd * test_op_dim + k;
|
||||
int m_offset = 0;
|
||||
for_constexpr<n_inputs>([&](auto s)
|
||||
{
|
||||
using fop_t = std::decay_t<decltype(get<s>(inputs))>;
|
||||
const int trial_op_dim =
|
||||
inputs_trial_op_dim[static_cast<int>(s)];
|
||||
if (trial_op_dim == 0) { return; }
|
||||
|
||||
const auto &in_dtq = input_dtq_maps[s];
|
||||
const auto eval_input =
|
||||
[&](const int m, const int axis, const int q,
|
||||
const int d)
|
||||
{
|
||||
if constexpr (is_value_fop<fop_t>::value)
|
||||
{
|
||||
return (m == 0) ? in_dtq.B(q, 0, d) : 0.0;
|
||||
}
|
||||
else if constexpr (is_gradient_fop<fop_t>::value)
|
||||
{
|
||||
return (m == axis) ? in_dtq.G(q, 0, d)
|
||||
: in_dtq.B(q, 0, d);
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
for (int m = 0; m < trial_op_dim; m++)
|
||||
{
|
||||
const int col = m_offset + m;
|
||||
backend_t::DiagContract(
|
||||
s_diag,
|
||||
num_test_dof_1d,
|
||||
q1d,
|
||||
nz_dof,
|
||||
[&](int axis, int q, int d)
|
||||
{ return eval_test(k, axis, q, d); },
|
||||
[&](int axis, int q, int d)
|
||||
{ return eval_input(m, axis, q, d); },
|
||||
[&](int q) { return qpdc(q, col, vd, row); },
|
||||
[&](int dx, int dy, int dz, real_t u)
|
||||
{ Y(dx, dy, dz) += u; });
|
||||
}
|
||||
m_offset += trial_op_dim;
|
||||
});
|
||||
return (k == 0) ? B(q, 0, d) : 0.0;
|
||||
}
|
||||
else if constexpr (is_gradient_fop<test_fop_t>::value)
|
||||
{
|
||||
return (k == axis) ? G(q, 0, d) : B(q, 0, d);
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
for (int k = 0; k < test_op_dim; k++)
|
||||
{
|
||||
const int row =
|
||||
out_offsets[static_cast<int>(o)] + vd * test_op_dim + k;
|
||||
int m_offset = 0;
|
||||
for_constexpr<n_inputs>([&](auto s)
|
||||
{
|
||||
using fop_t = std::decay_t<decltype(get<s>(inputs))>;
|
||||
const int trial_op_dim =
|
||||
inputs_trial_op_dim[static_cast<int>(s)];
|
||||
if (trial_op_dim == 0) { return; }
|
||||
|
||||
const auto &in_dtq = input_dtq_maps[s];
|
||||
const auto eval_input =
|
||||
[&](const int m, const int axis, const int q,
|
||||
const int d)
|
||||
{
|
||||
if constexpr (is_value_fop<fop_t>::value)
|
||||
{
|
||||
return (m == 0) ? in_dtq.B(q, 0, d) : 0.0;
|
||||
}
|
||||
else if constexpr (is_gradient_fop<fop_t>::value)
|
||||
{
|
||||
return (m == axis) ? in_dtq.G(q, 0, d)
|
||||
: in_dtq.B(q, 0, d);
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
for (int m = 0; m < trial_op_dim; m++)
|
||||
{
|
||||
const int col = m_offset + m;
|
||||
backend_t::DiagContract(
|
||||
s_diag,
|
||||
num_test_dof_1d,
|
||||
q1d,
|
||||
nz_dof,
|
||||
[&](int axis, int q, int d)
|
||||
{ return eval_test(k, axis, q, d); },
|
||||
[&](int axis, int q, int d)
|
||||
{ return eval_input(m, axis, q, d); },
|
||||
[&](int q) { return qpdc(q, col, vd, row); },
|
||||
[&](int dx, int dy, int dz, real_t u)
|
||||
{ Y(dx, dy, dz) += u; });
|
||||
}
|
||||
m_offset += trial_op_dim;
|
||||
});
|
||||
}
|
||||
});
|
||||
}
|
||||
|
||||
@@ -384,8 +384,6 @@ public:
|
||||
for_constexpr<n_outputs>([&](auto oc)
|
||||
{
|
||||
constexpr size_t o = oc.value, ao = n_inputs + o;
|
||||
using out_fop_t =
|
||||
std::decay_t<tuple_element_t<o, outputs_t>>;
|
||||
const auto &tangent = get<ao>(shadow_args);
|
||||
const int tv = out_vdim[o], to = out_op_dim[o];
|
||||
for (int i = 0; i < tv; i++)
|
||||
@@ -396,26 +394,8 @@ public:
|
||||
const int cache_idx =
|
||||
row * trial_vdim * total_trial_op_dim +
|
||||
j * total_trial_op_dim + col_m;
|
||||
// An Identity output is flat quadrature
|
||||
// point data: its FieldOperator vdim counts
|
||||
// components, not rows of the q-function
|
||||
// argument's shape. The two index form
|
||||
// assumes vdim == extents[0] and would run
|
||||
// off the end of, say, a
|
||||
// tensor<real_t, DIM, DIM> bound to a vdim
|
||||
// DIM*DIM space, so read it flat with the
|
||||
// same column major packing that
|
||||
// identity_qp_write_value writes.
|
||||
if constexpr (is_identity_fop_v<out_fop_t>)
|
||||
{
|
||||
cache_tensor(q, cache_idx, e) =
|
||||
qf_flat_value(tangent, i * to + k);
|
||||
}
|
||||
else
|
||||
{
|
||||
cache_tensor(q, cache_idx, e) =
|
||||
qf_value_at(tangent, i, k);
|
||||
}
|
||||
cache_tensor(q, cache_idx, e) =
|
||||
qf_value_at(tangent, i, k);
|
||||
}
|
||||
}
|
||||
});
|
||||
|
||||
+71
-209
@@ -109,29 +109,19 @@ struct derivative_action_t
|
||||
std::function<void *()> qfunc_shadow;
|
||||
};
|
||||
|
||||
/// @brief Type alias for a function that assembles the SparseMatrix row blocks
|
||||
/// of a derivative operator
|
||||
///
|
||||
/// The derivative of a residual with several output (test) fields is a block
|
||||
/// column, one row block per output field, so the callback fills a vector
|
||||
/// indexed by position in the operator's output FieldDescriptors rather than a
|
||||
/// single matrix.
|
||||
/// @brief Type alias for a function that assembles the SparseMatrix of a
|
||||
/// derivative operator
|
||||
using assemble_derivative_sparsematrix_callback_t =
|
||||
std::function<void(std::vector<SparseMatrix *> &)>;
|
||||
std::function<void(SparseMatrix *&)>;
|
||||
|
||||
/// @brief Type alias for a function that assembles the HypreParMatrix row
|
||||
/// blocks of a derivative operator
|
||||
///
|
||||
/// @see assemble_derivative_sparsematrix_callback_t for the indexing.
|
||||
/// @brief Type alias for a function that assembles the HypreParMatrix of a
|
||||
/// derivative operator
|
||||
using assemble_derivative_hypreparmatrix_callback_t =
|
||||
std::function<void(std::vector<HypreParMatrix *> &)>;
|
||||
std::function<void(HypreParMatrix *&)>;
|
||||
|
||||
/// @brief Type alias for a function that assembles the diagonal of one row
|
||||
/// block of a derivative operator into an E-vector
|
||||
///
|
||||
/// The first argument names the output (test) field whose row block is wanted;
|
||||
/// only a square block has a diagonal at all.
|
||||
using assemble_diagonal_callback_t = std::function<void(int, Vector &)>;
|
||||
/// @brief Type alias for a function that assembles the diagonal of a derivative
|
||||
/// operator into an E-vector
|
||||
using assemble_diagonal_callback_t = std::function<void(Vector &)>;
|
||||
|
||||
/// @brief Type alias for a function that applies the appropriate restriction to
|
||||
/// the solution and parameters
|
||||
@@ -258,73 +248,60 @@ MakeDerivativeHypreParMatrixAssemble(
|
||||
std::vector<assemble_derivative_sparsematrix_callback_t> &sparse_callbacks,
|
||||
const IntegratorContext &ctx)
|
||||
{
|
||||
using blocks_t = std::vector<HypreParMatrix *>;
|
||||
return [derivative_idx, &sparse_callbacks, ctx](blocks_t &A)
|
||||
return [derivative_idx, &sparse_callbacks, ctx](HypreParMatrix *&A)
|
||||
{
|
||||
MFEM_VERIFY(ctx.outfds.size() == 1,
|
||||
"HypreParMatrix assembly requires a single output field");
|
||||
|
||||
const size_t trial_field_idx = FindIdx(derivative_idx, ctx.unionfds);
|
||||
MFEM_VERIFY(trial_field_idx != SIZE_MAX,
|
||||
"derivative field not found for HypreParMatrix assembly");
|
||||
|
||||
const auto *test_fes_ptr =
|
||||
std::get_if<const ParFiniteElementSpace *>(&ctx.outfds[0].data);
|
||||
const auto *trial_fes_ptr =
|
||||
std::get_if<const ParFiniteElementSpace *>(
|
||||
&ctx.unionfds[trial_field_idx].data);
|
||||
MFEM_VERIFY(test_fes_ptr && *test_fes_ptr,
|
||||
"HypreParMatrix assembly requires a ParFiniteElementSpace "
|
||||
"output field");
|
||||
MFEM_VERIFY(trial_fes_ptr && *trial_fes_ptr,
|
||||
"HypreParMatrix assembly requires a ParFiniteElementSpace "
|
||||
"derivative field");
|
||||
const ParFiniteElementSpace *trial_fes = *trial_fes_ptr;
|
||||
|
||||
// Every integrator fills the row blocks it contributes to; the local
|
||||
// blocks are then RAP'd one by one, each with its own test space.
|
||||
std::vector<SparseMatrix *> spmat(ctx.outfds.size(), nullptr);
|
||||
const ParFiniteElementSpace *test_fes = *test_fes_ptr;
|
||||
const ParFiniteElementSpace *trial_fes = *trial_fes_ptr;
|
||||
MFEM_VERIFY(test_fes->GetComm() == trial_fes->GetComm(),
|
||||
"test and trial spaces must use the same MPI communicator");
|
||||
|
||||
SparseMatrix *spmat = nullptr;
|
||||
for (const auto &f : sparse_callbacks)
|
||||
{
|
||||
f(spmat);
|
||||
}
|
||||
|
||||
bool any = false;
|
||||
for (const auto *m : spmat) { any = any || (m != nullptr); }
|
||||
MFEM_VERIFY(any,
|
||||
MFEM_VERIFY(spmat != nullptr,
|
||||
"internal error: sparse derivative assembly returned NULL");
|
||||
MFEM_VERIFY(spmat->Finalized(),
|
||||
"local derivative matrix must be finalized");
|
||||
|
||||
if (A.size() < ctx.outfds.size())
|
||||
if (test_fes == trial_fes)
|
||||
{
|
||||
A.resize(ctx.outfds.size(), nullptr);
|
||||
HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
|
||||
test_fes->GetDofOffsets(), spmat);
|
||||
A = RAP(&dA, test_fes->Dof_TrueDof_Matrix());
|
||||
}
|
||||
else
|
||||
{
|
||||
HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
|
||||
trial_fes->GlobalVSize(),
|
||||
test_fes->GetDofOffsets(),
|
||||
trial_fes->GetDofOffsets(), spmat);
|
||||
A = RAP(test_fes->Dof_TrueDof_Matrix(), &dA,
|
||||
trial_fes->Dof_TrueDof_Matrix());
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < spmat.size(); i++)
|
||||
{
|
||||
if (spmat[i] == nullptr) { continue; }
|
||||
MFEM_VERIFY(spmat[i]->Finalized(),
|
||||
"local derivative matrix must be finalized");
|
||||
|
||||
const auto *test_fes_ptr =
|
||||
std::get_if<const ParFiniteElementSpace *>(&ctx.outfds[i].data);
|
||||
MFEM_VERIFY(test_fes_ptr && *test_fes_ptr,
|
||||
"HypreParMatrix assembly requires a ParFiniteElementSpace "
|
||||
"output field");
|
||||
const ParFiniteElementSpace *test_fes = *test_fes_ptr;
|
||||
MFEM_VERIFY(test_fes->GetComm() == trial_fes->GetComm(),
|
||||
"test and trial spaces must use the same "
|
||||
"MPI communicator");
|
||||
|
||||
if (test_fes == trial_fes)
|
||||
{
|
||||
HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
|
||||
test_fes->GetDofOffsets(), spmat[i]);
|
||||
A[i] = RAP(&dA, test_fes->Dof_TrueDof_Matrix());
|
||||
}
|
||||
else
|
||||
{
|
||||
HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
|
||||
trial_fes->GlobalVSize(),
|
||||
test_fes->GetDofOffsets(),
|
||||
trial_fes->GetDofOffsets(), spmat[i]);
|
||||
A[i] = RAP(test_fes->Dof_TrueDof_Matrix(), &dA,
|
||||
trial_fes->Dof_TrueDof_Matrix());
|
||||
}
|
||||
|
||||
delete spmat[i];
|
||||
}
|
||||
delete spmat;
|
||||
};
|
||||
}
|
||||
|
||||
@@ -586,64 +563,17 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief Assemble the row blocks of the derivative operator into
|
||||
/// SparseMatrices.
|
||||
///
|
||||
/// A q-function may write to several output (test) fields. Differentiating
|
||||
/// w.r.t. one field then gives a block column with one row block per output
|
||||
/// field,
|
||||
///
|
||||
/// dR/dU = [ dR_U/dU ; dR_Y/dU ],
|
||||
///
|
||||
/// all sharing the trial space of the differentiated field. Each row block
|
||||
/// is an ordinary test x trial matrix and is assembled separately; the
|
||||
/// stacked matrix, if wanted, is one HypreParMatrixFromBlocks call away.
|
||||
///
|
||||
/// @param A Resized to the number of output fields and indexed the same way,
|
||||
/// so A[f] is the row block of output field f. Fields whose rows cannot be
|
||||
/// materialised -- quadrature and parameter spaces, which have no basis to
|
||||
/// contract against -- are left as nullptr. Ownership passes to the caller.
|
||||
void Assemble(std::vector<SparseMatrix *> &A)
|
||||
{
|
||||
MFEM_VERIFY(!assemble_derivative_sparsematrix_callbacks.empty(),
|
||||
"derivative can't be assembled into a SparseMatrix");
|
||||
EnsureQpCache();
|
||||
|
||||
A.assign(outfds.size(), nullptr);
|
||||
for (const auto &f : assemble_derivative_sparsematrix_callbacks)
|
||||
{
|
||||
f(A);
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief Assemble the derivative operator into a SparseMatrix.
|
||||
///
|
||||
/// Convenience overload for the common case of a derivative with a single
|
||||
/// assemblable row block. Use Assemble(std::vector<SparseMatrix *> &) when
|
||||
/// the q-function writes to more than one test field.
|
||||
///
|
||||
/// @param A The SparseMatrix to assemble the derivative operator into. Can
|
||||
/// be an uninitialized object.
|
||||
void Assemble(SparseMatrix *&A)
|
||||
{
|
||||
std::vector<SparseMatrix *> blocks;
|
||||
Assemble(blocks);
|
||||
A = SingleBlock(blocks);
|
||||
}
|
||||
|
||||
/// @brief Assemble the row blocks of the derivative operator into
|
||||
/// HypreParMatrices.
|
||||
///
|
||||
/// @see Assemble(std::vector<SparseMatrix *> &) for the indexing and the
|
||||
/// block structure.
|
||||
void Assemble(std::vector<HypreParMatrix *> &A)
|
||||
{
|
||||
MFEM_VERIFY(!assemble_derivative_hypreparmatrix_callbacks.empty(),
|
||||
"derivative can't be assembled into a HypreParMatrix");
|
||||
MFEM_ASSERT(!assemble_derivative_sparsematrix_callbacks.empty(),
|
||||
"derivative can't be assembled into a SparseMatrix");
|
||||
EnsureQpCache();
|
||||
|
||||
A.assign(outfds.size(), nullptr);
|
||||
for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
|
||||
for (const auto &f : assemble_derivative_sparsematrix_callbacks)
|
||||
{
|
||||
f(A);
|
||||
}
|
||||
@@ -651,17 +581,18 @@ public:
|
||||
|
||||
/// @brief Assemble the derivative operator into a HypreParMatrix.
|
||||
///
|
||||
/// Convenience overload for the common case of a derivative with a single
|
||||
/// assemblable row block. Use Assemble(std::vector<HypreParMatrix *> &) when
|
||||
/// the q-function writes to more than one test field.
|
||||
///
|
||||
/// @param A The HypreParMatrix to assemble the derivative operator into. Can
|
||||
/// be an uninitialized object.
|
||||
void Assemble(HypreParMatrix *&A)
|
||||
{
|
||||
std::vector<HypreParMatrix *> blocks;
|
||||
Assemble(blocks);
|
||||
A = SingleBlock(blocks);
|
||||
MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
|
||||
"derivative can't be assembled into a HypreParMatrix");
|
||||
EnsureQpCache();
|
||||
|
||||
for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
|
||||
{
|
||||
f(A);
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief Assemble the derivative of a functional into a Vector.
|
||||
@@ -704,89 +635,35 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief Assemble the diagonal of one row block of the derivative
|
||||
/// operator into a T-vector.
|
||||
/// @brief Assemble the diagonal of the derivative operator into a T-vector.
|
||||
///
|
||||
/// The derivative is a block column, one row block per output field, all
|
||||
/// sharing the trial space of the differentiated field. Only a block whose
|
||||
/// test space *is* that trial space is square, and only a square block has a
|
||||
/// diagonal. Squareness alone cannot pick the block when several output
|
||||
/// fields live on the trial space, so the row is named rather than inferred.
|
||||
/// Output field ids may differ from the differentiated field's id while
|
||||
/// sharing its ParFiniteElementSpace, in which case several row blocks are
|
||||
/// square at once. Naming the output field id resolves the ambiguity.
|
||||
///
|
||||
/// @param out_field_id The output field whose row block is wanted.
|
||||
/// @param diag The vector to receive the diagonal, resized to that field's
|
||||
/// T-dof size.
|
||||
void AssembleDiagonal(size_t out_field_id, Vector &diag) const
|
||||
/// @param diag The vector to receive the diagonal (must be T-dof sized).
|
||||
void AssembleDiagonal(Vector &diag) const override
|
||||
{
|
||||
MFEM_VERIFY(!assemble_diagonal_callbacks.empty(),
|
||||
MFEM_ASSERT(!assemble_diagonal_callbacks.empty(),
|
||||
"derivative can't assemble diagonal");
|
||||
EnsureQpCache();
|
||||
|
||||
const size_t diagonal_idx = FindIdx(out_field_id, outfds);
|
||||
MFEM_VERIFY(diagonal_idx != SIZE_MAX,
|
||||
"AssembleDiagonal: field " << out_field_id << " is not an "
|
||||
"output field of this operator");
|
||||
MFEM_ASSERT(outfds.size() == 1,
|
||||
"AssembleDiagonal currently requires a single output field");
|
||||
|
||||
const auto *test_pf =
|
||||
std::get_if<const ParFiniteElementSpace *>(&outfds[diagonal_idx].data);
|
||||
std::get_if<const ParFiniteElementSpace *>(&outfds[0].data);
|
||||
MFEM_VERIFY(test_pf && *test_pf,
|
||||
"AssembleDiagonal: test field must be a ParFiniteElementSpace");
|
||||
const auto *trial_pf =
|
||||
std::get_if<const ParFiniteElementSpace *>(&direction.data);
|
||||
MFEM_VERIFY(trial_pf && *trial_pf,
|
||||
"AssembleDiagonal: the differentiated field must be a "
|
||||
"ParFiniteElementSpace");
|
||||
MFEM_VERIFY(*test_pf == *trial_pf,
|
||||
"AssembleDiagonal: the requested row block is not square and "
|
||||
"so has no diagonal; its test space is not the trial space "
|
||||
"of the differentiated field");
|
||||
|
||||
prepare_residual(outfds, out_rcache, daction_e);
|
||||
for (auto *v : daction_e) { *v = 0.0; }
|
||||
|
||||
for (const auto &f : assemble_diagonal_callbacks)
|
||||
{
|
||||
f(static_cast<int>(out_field_id), *daction_e[diagonal_idx]);
|
||||
f(*daction_e[0]);
|
||||
}
|
||||
|
||||
restriction_transpose(outfds, out_rcache, daction_e, daction_l);
|
||||
prolongation_transpose(outfds[diagonal_idx],
|
||||
*daction_l[diagonal_idx], diag);
|
||||
}
|
||||
|
||||
/// @brief Assemble the diagonal of the derivative operator into a T-vector.
|
||||
///
|
||||
/// Uses the row block of the differentiated field, i.e. the usual diagonal
|
||||
/// block dR_U/dU for GetDerivative(U). Call
|
||||
/// AssembleDiagonal(size_t, Vector &) for any other row block.
|
||||
///
|
||||
/// @param diag The vector to receive the diagonal.
|
||||
void AssembleDiagonal(Vector &diag) const override
|
||||
{
|
||||
AssembleDiagonal(direction.id, diag);
|
||||
prolongation_transpose(outfds[0], *daction_l[0], diag);
|
||||
}
|
||||
|
||||
private:
|
||||
/// The single non-null row block of @a blocks, for the scalar Assemble
|
||||
/// overloads.
|
||||
template <typename mat_t>
|
||||
static mat_t *SingleBlock(const std::vector<mat_t *> &blocks)
|
||||
{
|
||||
mat_t *single = nullptr;
|
||||
int count = 0;
|
||||
for (auto *m : blocks)
|
||||
{
|
||||
if (m != nullptr) { single = m; count++; }
|
||||
}
|
||||
MFEM_VERIFY(count == 1,
|
||||
"the derivative has " << count << " assemblable row blocks, "
|
||||
"not one; assemble it into a std::vector instead");
|
||||
return single;
|
||||
}
|
||||
|
||||
/// Derivative action callbacks. Depending on the requested derivatives in
|
||||
/// DifferentiableOperator the callbacks represent certain combinations of
|
||||
/// actions of derivatives of the forward operator.
|
||||
@@ -1602,23 +1479,16 @@ void DifferentiableOperator::AddIntegrator(
|
||||
constexpr size_t derivative_idx = decltype(derivative_id)::value;
|
||||
using callback_outputs_t = std::decay_t<decltype(outputs)>;
|
||||
|
||||
|
||||
// NOTE: before disable_assemble was looking for any identity outputs,
|
||||
// But this would silently disable assembly for other valid outputs across
|
||||
// different fields.
|
||||
bool any_assemblable_output = false;
|
||||
bool disable_assemble = false;
|
||||
for_constexpr([&](auto j)
|
||||
{
|
||||
using output_fop_t =
|
||||
std::decay_t<tuple_element_t<j, callback_outputs_t>>;
|
||||
if constexpr (!is_identity_fop_v<output_fop_t>)
|
||||
using output_fop_t = tuple_element_t<j, callback_outputs_t>;
|
||||
if constexpr (is_identity_fop_v<std::decay_t<output_fop_t>>)
|
||||
{
|
||||
any_assemblable_output = true;
|
||||
disable_assemble = true;
|
||||
}
|
||||
}, std::make_index_sequence<tuple_size<callback_outputs_t>::value> {});
|
||||
|
||||
const bool disable_assemble = !any_assemblable_output;
|
||||
|
||||
// Setup the qp cache for the derivative
|
||||
setup_callbacks[callback_key].push_back(
|
||||
MakeDerivativeSetupCallback(
|
||||
@@ -1644,21 +1514,13 @@ void DifferentiableOperator::AddIntegrator(
|
||||
backend_t::template MakeDerivativeAssemble<derivative_idx>(
|
||||
callback_ctx, qf, inputs, outputs, callback_qp_cache));
|
||||
|
||||
// Assemble the derivative into a HypreParMatrix. This one runs
|
||||
// every sparse callback registered under the key, so it is
|
||||
// registered once and not once per integrator.
|
||||
if (assemble_hypreparmatrix_callbacks[callback_key].empty())
|
||||
{
|
||||
assemble_hypreparmatrix_callbacks[callback_key].push_back(
|
||||
MakeDerivativeHypreParMatrixAssemble(
|
||||
derivative_idx,
|
||||
assemble_sparsematrix_callbacks[callback_key],
|
||||
callback_ctx));
|
||||
}
|
||||
}
|
||||
// Assemble the derivative into a HypreParMatrix
|
||||
assemble_hypreparmatrix_callbacks[callback_key].push_back(
|
||||
MakeDerivativeHypreParMatrixAssemble(
|
||||
derivative_idx,
|
||||
assemble_sparsematrix_callbacks[callback_key],
|
||||
callback_ctx));
|
||||
|
||||
if (!disable_assemble)
|
||||
{
|
||||
// Assemble the diagonal of the derivative into an L-vector
|
||||
assemble_diagonal_cbs[callback_key].push_back(
|
||||
backend_t::template MakeDerivativeAssembleDiagonal<derivative_idx>(
|
||||
|
||||
@@ -14,20 +14,8 @@ if (MFEM_USE_MPI)
|
||||
MAIN dfem-minimal-surface.cpp
|
||||
LIBRARIES mfem)
|
||||
|
||||
add_mfem_miniapp(dfem-hyperelasticity
|
||||
MAIN dfem-hyperelasticity.cpp
|
||||
LIBRARIES mfem)
|
||||
|
||||
add_mfem_miniapp(dfem-nonlinear-poisson
|
||||
MAIN dfem-nonlinear-poisson.cpp
|
||||
LIBRARIES mfem)
|
||||
|
||||
add_mfem_miniapp(dfem-hyperbolic-heat
|
||||
MAIN dfem-hyperbolic-heat.cpp
|
||||
LIBRARIES mfem)
|
||||
|
||||
add_mfem_miniapp(dfem-reaction-diffusion
|
||||
MAIN dfem-reaction-diffusion.cpp
|
||||
add_mfem_miniapp(dfem-hyperelasticity-energy
|
||||
MAIN dfem-hyperelasticity_energy.cpp
|
||||
LIBRARIES mfem)
|
||||
|
||||
# Add parallel tests.
|
||||
@@ -37,35 +25,9 @@ if (MFEM_USE_MPI)
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-minimal-surface> -der 0 -o 1 -r 2 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
|
||||
add_test(NAME dfem-hyperelasticity=${MFEM_MPI_NP}
|
||||
add_test(NAME dfem-hyperelasticity-energy=${MFEM_MPI_NP}
|
||||
COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperelasticity> -o 1 -rs 0 -no-vis
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperelasticity-energy> -o 1 -rs 0 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
|
||||
add_test(NAME dfem-nonlinear-poisson=${MFEM_MPI_NP}
|
||||
COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-nonlinear-poisson>
|
||||
-k grad -f energy -o 1 -r 2 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
|
||||
add_test(NAME dfem-nonlinear-poisson-residual=${MFEM_MPI_NP}
|
||||
COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-nonlinear-poisson>
|
||||
-k sol -f residual -o 1 -r 2 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
|
||||
add_test(NAME dfem-reaction-diffusion=${MFEM_MPI_NP}
|
||||
COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-reaction-diffusion>
|
||||
-b 1 -o 1 -r 2 -pc 2 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
|
||||
foreach(AL partial legacy full)
|
||||
add_test(NAME dfem-hyperbolic-heat-${AL}=${MFEM_MPI_NP}
|
||||
COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
|
||||
${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperbolic-heat>
|
||||
-al ${AL} -o 1 -r 0 -tf 0.1 -dt 0.05 -no-vis
|
||||
${MPIEXEC_POSTFLAGS})
|
||||
endforeach()
|
||||
endif()
|
||||
endif()
|
||||
|
||||
@@ -1,839 +0,0 @@
|
||||
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
|
||||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
//
|
||||
// This file is part of the MFEM library. For more information and source code
|
||||
// availability visit https://mfem.org.
|
||||
//
|
||||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
// CONTRIBUTING.md for details.
|
||||
//
|
||||
// --------------------------------------------
|
||||
// Nonlinear Hyperbolic Heat Equation with dFEM
|
||||
// --------------------------------------------
|
||||
//
|
||||
// Compile with: make dfem-hyperbolic-heat
|
||||
//
|
||||
// Sample runs: mpirun -np 4 dfem-hyperbolic-heat
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -b 0
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -tau 0.02
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -q 5 -rc 4
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -r 1 -tf 0.6
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -dt 0.2 -tf 2
|
||||
// mpirun -np 4 dfem-hyperbolic-heat -al legacy
|
||||
//
|
||||
// Description: This miniapp solves the hyperbolic (Maxwell-Cattaneo) heat
|
||||
// equation with a nonlinear conductivity,
|
||||
//
|
||||
// tau d^2T/dt^2 + rho c dT/dt = div( k(T) grad T ) + Q,
|
||||
//
|
||||
// on the unit cube (set the compile-time constant dim to 2 for
|
||||
// the unit square), starting from a Gaussian hot spot at rest.
|
||||
// The relaxation term tau d^2T/dt^2 turns the parabolic heat
|
||||
// equation into a damped wave equation to comply with the
|
||||
// finite speed sqrt(k/tau) of heat propagation.
|
||||
// Small tau recovers diffusive behavior (Fourier's law),
|
||||
// large tau a wave-like behavior.
|
||||
//
|
||||
// The conductivity k(T) = k0 (1 + beta T^2) is the only
|
||||
// nonlinearity. It is written as a dFEM q-function, so both the
|
||||
// diffusion residual
|
||||
//
|
||||
// F(T) = int k(T) grad T . grad w dx
|
||||
//
|
||||
// and its tangent J_F(T) come from the same pointwise
|
||||
// description; nothing here differentiates k(T) by hand.
|
||||
//
|
||||
// We recommend viewing examples 16 and 23, and dfem miniapps
|
||||
// before viewing this example.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "../../fem/dfem/doperator.hpp"
|
||||
#include "../../fem/dfem/backends/local_qf/prelude.hpp"
|
||||
|
||||
#include <memory>
|
||||
|
||||
using namespace mfem;
|
||||
|
||||
// This example code demonstrates the use of new features in MFEM that are in
|
||||
// development but exposed through the mfem::future namespace. All features
|
||||
// under this namespace might change their interface or behavior in upcoming
|
||||
// releases until they have stabilized.
|
||||
using namespace mfem::future;
|
||||
using mfem::future::tensor;
|
||||
|
||||
#ifdef MFEM_USE_ENZYME
|
||||
using dscalar_t = real_t;
|
||||
#else
|
||||
using mfem::future::dual;
|
||||
using dscalar_t = dual<real_t, real_t>;
|
||||
#endif
|
||||
|
||||
// Space dimension. The q-function and the diffusion operator are templated on
|
||||
// it, so it is fixed at compile time; set it to 2 to run the same problem on
|
||||
// the unit square.
|
||||
constexpr int dim = 3;
|
||||
|
||||
// Field IDs used by the dFEM integrator.
|
||||
static constexpr int Temperature = 0;
|
||||
static constexpr int Coords = 1;
|
||||
|
||||
enum class PreconditionerType { None, Diagonal };
|
||||
|
||||
/// Map an AssemblyLevel name onto the MFEM assembly level used for the mass operators.
|
||||
/// Only partial assembly keeps the whole solve matrix free; the others are here
|
||||
/// so the cost of that choice can be measured against the alternatives.
|
||||
static AssemblyLevel ParseAssemblyLevel(const char *name)
|
||||
{
|
||||
const std::string s(name);
|
||||
if (s == "partial") { return AssemblyLevel::PARTIAL; }
|
||||
if (s == "legacy") { return AssemblyLevel::LEGACY; }
|
||||
if (s == "full") { return AssemblyLevel::FULL; }
|
||||
MFEM_ABORT("Unknown assembly level '" << name << "'. Available levels: "
|
||||
"partial, legacy, full.");
|
||||
return AssemblyLevel::LEGACY;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Pointwise description of the nonlinear diffusion term.
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
/// Quadrature point kernel for F(T) = int k(T) grad T . grad w dx with the
|
||||
/// temperature dependent conductivity k(T) = k0 (1 + beta T^2).
|
||||
template <int DIM>
|
||||
struct NonlinearDiffusionQFunction
|
||||
{
|
||||
real_t k0 = 1.0;
|
||||
real_t beta = 1.0;
|
||||
|
||||
MFEM_HOST_DEVICE inline
|
||||
auto operator()(const dscalar_t &T,
|
||||
const tensor<dscalar_t, DIM> &dTdxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
tensor<dscalar_t, DIM> &dwdx) const
|
||||
{
|
||||
const auto invJ = inv(J);
|
||||
const auto dTdx = dTdxi * invJ;
|
||||
const dscalar_t k = k0 * (1.0_r + beta * T * T);
|
||||
dwdx = k * dTdx * transpose(invJ) * det(J) * w;
|
||||
}
|
||||
};
|
||||
|
||||
/// dFEM wrapper for the nonlinear diffusion residual F(T) and its tangent
|
||||
/// J_F(T).
|
||||
///
|
||||
/// Neither F nor J_F applies any essential-dof treatment. That belongs to the
|
||||
/// operator that assembles the full residual, since the boundary condition
|
||||
/// constrains the acceleration, not the diffusion term by itself.
|
||||
template <int DIM>
|
||||
class NonlinearDiffusionOperator
|
||||
{
|
||||
public:
|
||||
NonlinearDiffusionOperator(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
real_t k0, real_t beta)
|
||||
{
|
||||
auto &mesh_nodes =
|
||||
*static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
|
||||
mesh_nodes_fes = mesh_nodes.ParFESpace();
|
||||
mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
|
||||
|
||||
const std::vector<FieldDescriptor> inputs =
|
||||
{
|
||||
{Temperature, &fes}, {Coords, mesh_nodes_fes}
|
||||
};
|
||||
// The output operator is the gradient of the test function basis, which
|
||||
// completes the weak form B^T D(B T, B x, w).
|
||||
const std::vector<FieldDescriptor> outputs = {{Temperature, &fes}};
|
||||
|
||||
dop = std::make_shared<DifferentiableOperator>(
|
||||
inputs, outputs, *fes.GetParMesh());
|
||||
|
||||
Array<int> all_domain_attr;
|
||||
if (fes.GetMesh()->attributes.Size() > 0)
|
||||
{
|
||||
all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
|
||||
all_domain_attr = 1;
|
||||
}
|
||||
|
||||
NonlinearDiffusionQFunction<DIM> qf;
|
||||
qf.k0 = k0;
|
||||
qf.beta = beta;
|
||||
|
||||
// Requesting the derivative with respect to Temperature is what makes
|
||||
// J_F available later through GetDerivative.
|
||||
auto derivatives = std::integer_sequence<size_t, Temperature> {};
|
||||
dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
Inputs < Value<Temperature>, Gradient<Temperature>,
|
||||
Gradient<Coords>, Weight > {},
|
||||
Outputs<Gradient<Temperature>> {},
|
||||
ir, all_domain_attr, derivatives);
|
||||
}
|
||||
|
||||
/// F(T), the nonlinear diffusion residual.
|
||||
void Mult(const Vector &T, Vector &F) const
|
||||
{
|
||||
MultiVector X{T, mesh_nodes_tdofs};
|
||||
MultiVector Y{F};
|
||||
dop->Mult(X, Y);
|
||||
}
|
||||
|
||||
/// J_F(T), the tangent of F linearized at @a T. The returned operator
|
||||
/// captures @a T, so it has to be rebuilt whenever the linearization point
|
||||
/// moves.
|
||||
std::shared_ptr<DerivativeOperator> GetGradient(const Vector &T) const
|
||||
{
|
||||
MultiVector X{T, mesh_nodes_tdofs};
|
||||
return dop->GetDerivative(Temperature, X);
|
||||
}
|
||||
|
||||
private:
|
||||
ParFiniteElementSpace *mesh_nodes_fes = nullptr;
|
||||
Vector mesh_nodes_tdofs;
|
||||
std::shared_ptr<DifferentiableOperator> dop;
|
||||
};
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Newton residual for one implicit stage.
|
||||
// ----------------------------------------------------------------------------
|
||||
///
|
||||
/// Residual of the semi-discrete system, written in the acceleration unknown
|
||||
/// a = d^2T/dt^2 that MFEM's second order ODE solvers solve for. With the stage
|
||||
/// quantities requested by the integrator,
|
||||
///
|
||||
/// T_stage = T + fac0 a, v_stage = dT/dt + fac1 a,
|
||||
///
|
||||
/// the weak form M_tau a + M_rhoc v + F(T) - Q = 0 becomes
|
||||
///
|
||||
/// R(a) = M_tau a + M_rhoc v_stage + F(T_stage) - Q.
|
||||
///
|
||||
class ResidualOperator : public Operator
|
||||
{
|
||||
public:
|
||||
/// Newton tangent
|
||||
///
|
||||
/// J_R(a) = M_tau + fac1 M_rhoc + fac0 J_F(T + fac0 a).
|
||||
///
|
||||
/// None of the three terms is a matrix: the mass operators are partially
|
||||
/// assembled and J_F is a dFEM derivative action, so this class only ever
|
||||
/// composes actions and diagonals.
|
||||
class JacobianOperator : public Operator
|
||||
{
|
||||
public:
|
||||
JacobianOperator(const ResidualOperator &oper, const Vector &T_stage) :
|
||||
Operator(oper.Height()), oper(oper), z(oper.Height()), w(oper.Height())
|
||||
{
|
||||
diffusion_tangent = oper.diffusion.GetGradient(T_stage);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override
|
||||
{
|
||||
// Essential directions are removed before the apply and restored as
|
||||
// identity rows afterwards.
|
||||
z = x;
|
||||
z.SetSubVector(oper.ess_tdofs, 0.0);
|
||||
|
||||
oper.M_tau.Mult(z, y);
|
||||
oper.M_rhoc.AddMult(z, y, oper.fac1);
|
||||
|
||||
MultiVector W{w};
|
||||
diffusion_tangent->Mult(z, W);
|
||||
y.Add(oper.fac0, w);
|
||||
|
||||
auto d_y = y.ReadWrite();
|
||||
const auto d_x = x.Read();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_y[d_dofs[i]] = d_x[d_dofs[i]];
|
||||
});
|
||||
}
|
||||
|
||||
void AssembleDiagonal(Vector &diag) const override
|
||||
{
|
||||
oper.M_tau.AssembleDiagonal(diag);
|
||||
|
||||
oper.M_rhoc.AssembleDiagonal(w);
|
||||
diag.Add(oper.fac1, w);
|
||||
|
||||
diffusion_tangent->AssembleDiagonal(w);
|
||||
diag.Add(oper.fac0, w);
|
||||
|
||||
auto d_diag = diag.ReadWrite();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_diag[d_dofs[i]] = 1.0;
|
||||
});
|
||||
}
|
||||
|
||||
private:
|
||||
const ResidualOperator &oper;
|
||||
mutable Vector z, w;
|
||||
std::shared_ptr<DerivativeOperator> diffusion_tangent;
|
||||
};
|
||||
|
||||
/// @a M_tau and @a M_rhoc are the mass operators as returned by
|
||||
/// ParBilinearForm::FormSystemMatrix, i.e. constrained true-dof operators.
|
||||
ResidualOperator(const Operator &M_tau, const Operator &M_rhoc,
|
||||
const NonlinearDiffusionOperator<dim> &diffusion,
|
||||
const Array<int> &ess_tdofs) :
|
||||
Operator(M_tau.Height()),
|
||||
M_tau(M_tau), M_rhoc(M_rhoc), diffusion(diffusion), ess_tdofs(ess_tdofs),
|
||||
T_stage(height), v_stage(height), z(height) { }
|
||||
|
||||
/// Describe the stage the time integrator is asking about. @a T, @a dTdt and
|
||||
/// @a Q are referenced, not copied, so they have to outlive the Newton solve.
|
||||
void SetParameters(real_t fac0_, real_t fac1_, const Vector *T_,
|
||||
const Vector *dTdt_, const Vector *Q_)
|
||||
{
|
||||
fac0 = fac0_;
|
||||
fac1 = fac1_;
|
||||
T = T_;
|
||||
dTdt = dTdt_;
|
||||
Q = Q_;
|
||||
}
|
||||
|
||||
void Mult(const Vector &a, Vector &R) const override
|
||||
{
|
||||
MFEM_ASSERT(T && dTdt && Q, "call SetParameters() first");
|
||||
|
||||
add(*T, fac0, a, T_stage);
|
||||
add(*dTdt, fac1, a, v_stage);
|
||||
|
||||
M_tau.Mult(a, R);
|
||||
M_rhoc.AddMult(v_stage, R);
|
||||
|
||||
diffusion.Mult(T_stage, z);
|
||||
R += z;
|
||||
|
||||
R -= *Q;
|
||||
|
||||
R.SetSubVector(ess_tdofs, 0.0);
|
||||
}
|
||||
|
||||
Operator &GetGradient(const Vector &a) const override
|
||||
{
|
||||
// The tangent has to be linearized at the stage temperature, not at the
|
||||
// temperature of the previous step.
|
||||
add(*T, fac0, a, T_stage);
|
||||
jacobian = std::make_shared<JacobianOperator>(*this, T_stage);
|
||||
return *jacobian;
|
||||
}
|
||||
|
||||
private:
|
||||
const Operator &M_tau;
|
||||
const Operator &M_rhoc;
|
||||
const NonlinearDiffusionOperator<dim> &diffusion;
|
||||
const Array<int> &ess_tdofs;
|
||||
|
||||
real_t fac0 = 0.0, fac1 = 0.0;
|
||||
const Vector *T = nullptr;
|
||||
const Vector *dTdt = nullptr;
|
||||
const Vector *Q = nullptr;
|
||||
|
||||
mutable Vector T_stage, v_stage, z;
|
||||
mutable std::shared_ptr<JacobianOperator> jacobian;
|
||||
};
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// The second order time dependent operator.
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
/// Semi-discrete hyperbolic heat equation
|
||||
///
|
||||
/// M_tau a + M_rhoc v + F(T) - Q = 0,
|
||||
///
|
||||
class HyperbolicHeatOperator : public SecondOrderTimeDependentOperator
|
||||
{
|
||||
public:
|
||||
HyperbolicHeatOperator(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
const Array<int> &ess_bdr,
|
||||
real_t tau, real_t rho_c, real_t k0, real_t beta,
|
||||
real_t source, PreconditionerType pc_type,
|
||||
AssemblyLevel assembly);
|
||||
|
||||
/// Compute a = M_tau^{-1} ( Q - M_rhoc v - F(T) ).
|
||||
using SecondOrderTimeDependentOperator::Mult;
|
||||
void Mult(const Vector &T, const Vector &dTdt,
|
||||
Vector &d2Tdt2) const override;
|
||||
|
||||
/// Solve R(a) = M_tau a + M_rhoc (dTdt + fac1 a) + F(T + fac0 a) - Q = 0 for
|
||||
/// the unknown acceleration a = d2Tdt2.
|
||||
using SecondOrderTimeDependentOperator::ImplicitSolve;
|
||||
void ImplicitSolve(const real_t fac0, const real_t fac1,
|
||||
const Vector &T, const Vector &dTdt,
|
||||
Vector &d2Tdt2) override;
|
||||
|
||||
const Array<int> &GetEssentialTrueDofs() const { return ess_tdof_list; }
|
||||
|
||||
/// Newton iterations taken by the last ImplicitSolve.
|
||||
int GetNewtonIterations() const { return newton_solver.GetNumIterations(); }
|
||||
|
||||
/// Compute the physical L2 norm directly from a true-dof vector.
|
||||
real_t ComputeL2Norm(const Vector &x) const;
|
||||
|
||||
~HyperbolicHeatOperator() override;
|
||||
|
||||
private:
|
||||
/// Set up one mass operator weighted by the constant coefficient @a w at the
|
||||
/// requested @a assembly level, and hand back its constrained true-dof form
|
||||
/// in @a M.
|
||||
void SetupMass(ParBilinearForm &m, ConstantCoefficient &w,
|
||||
const IntegrationRule &ir, AssemblyLevel assembly,
|
||||
OperatorHandle &M);
|
||||
|
||||
ParFiniteElementSpace &fes;
|
||||
Array<int> ess_tdof_list; // empty for pure Neumann boundary conditions
|
||||
|
||||
// Mass operators, weighted by tau and by rho*c. The forms have to outlive
|
||||
|
||||
// the handles: under partial assembly the operator in the handle refers to
|
||||
// the quadrature point data held by the form.
|
||||
ConstantCoefficient tau_coeff, rhoc_coeff;
|
||||
real_t tau;
|
||||
ParBilinearForm m_tau, m_rhoc;
|
||||
OperatorHandle M_tau, M_rhoc;
|
||||
|
||||
std::unique_ptr<NonlinearDiffusionOperator<dim>> diffusion;
|
||||
Vector Q;
|
||||
|
||||
CGSolver M_solver;
|
||||
OperatorJacobiSmoother M_prec;
|
||||
|
||||
std::unique_ptr<ResidualOperator> residual_op;
|
||||
NewtonSolver newton_solver;
|
||||
std::unique_ptr<IterativeSolver> krylov;
|
||||
std::unique_ptr<Solver> krylov_prec;
|
||||
|
||||
mutable Vector cached_accel;
|
||||
mutable Vector rhs, z;
|
||||
};
|
||||
|
||||
static std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
|
||||
{
|
||||
switch (type)
|
||||
{
|
||||
case PreconditionerType::None: return nullptr;
|
||||
case PreconditionerType::Diagonal:
|
||||
return std::make_unique<OperatorJacobiSmoother>();
|
||||
default:
|
||||
MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
void HyperbolicHeatOperator::SetupMass(ParBilinearForm &m,
|
||||
ConstantCoefficient &w,
|
||||
const IntegrationRule &ir,
|
||||
AssemblyLevel assembly,
|
||||
OperatorHandle &M)
|
||||
{
|
||||
m.SetAssemblyLevel(assembly);
|
||||
m.AddDomainIntegrator(new MassIntegrator(w, &ir));
|
||||
m.Assemble();
|
||||
m.FormSystemMatrix(ess_tdof_list, M);
|
||||
}
|
||||
|
||||
HyperbolicHeatOperator::HyperbolicHeatOperator(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
const Array<int> &ess_bdr,
|
||||
real_t tau, real_t rho_c,
|
||||
real_t k0, real_t beta,
|
||||
real_t source,
|
||||
PreconditionerType pc_type,
|
||||
AssemblyLevel assembly)
|
||||
: SecondOrderTimeDependentOperator(fes.GetTrueVSize(), (real_t) 0.0),
|
||||
fes(fes),
|
||||
tau_coeff(tau), rhoc_coeff(rho_c),
|
||||
tau(tau),
|
||||
m_tau(&fes), m_rhoc(&fes),
|
||||
M_solver(fes.GetComm()),
|
||||
newton_solver(fes.GetComm()),
|
||||
cached_accel(height), rhs(height), z(height)
|
||||
{
|
||||
fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
|
||||
|
||||
/// Setup the mass operators for the acceleration and velocity terms.
|
||||
SetupMass(m_tau, tau_coeff, ir, assembly, M_tau);
|
||||
SetupMass(m_rhoc, rhoc_coeff, ir, assembly, M_rhoc);
|
||||
|
||||
/// Create the nonlinear diffusion operator as a dFEM operator
|
||||
diffusion = std::make_unique<NonlinearDiffusionOperator<dim>>(fes, ir, k0,
|
||||
beta);
|
||||
|
||||
/// The source term is a constant coefficient, so we can assemble it once and for all.
|
||||
ConstantCoefficient source_coeff(source);
|
||||
ParLinearForm b(&fes);
|
||||
b.AddDomainIntegrator(new DomainLFIntegrator(source_coeff, &ir));
|
||||
b.Assemble();
|
||||
Q.SetSize(height);
|
||||
b.ParallelAssemble(Q);
|
||||
|
||||
M_solver.iterative_mode = false;
|
||||
M_solver.SetRelTol(1e-10);
|
||||
M_solver.SetAbsTol(0.0);
|
||||
M_solver.SetMaxIter(200);
|
||||
M_solver.SetPrintLevel(0);
|
||||
M_solver.SetPreconditioner(M_prec);
|
||||
M_solver.SetOperator(*M_tau);
|
||||
|
||||
/// Set up the nonlinear residual operator used by the Newton solver.
|
||||
residual_op = std::make_unique<ResidualOperator>(*M_tau, *M_rhoc, *diffusion,
|
||||
ess_tdof_list);
|
||||
|
||||
/// Set up the Newton solver and Krylov solver for linearized systems.
|
||||
if (beta == 0.0)
|
||||
{
|
||||
krylov = std::make_unique<CGSolver>(fes.GetComm());
|
||||
}
|
||||
else
|
||||
{
|
||||
krylov = std::make_unique<GMRESSolver>(fes.GetComm());
|
||||
}
|
||||
krylov->SetRelTol(1e-10);
|
||||
krylov->SetAbsTol(0.0);
|
||||
krylov->SetMaxIter(500);
|
||||
krylov->SetPrintLevel(0);
|
||||
|
||||
krylov_prec = MakePreconditioner(pc_type);
|
||||
if (krylov_prec) { krylov->SetPreconditioner(*krylov_prec); }
|
||||
|
||||
newton_solver.SetOperator(*residual_op);
|
||||
newton_solver.SetSolver(*krylov);
|
||||
#ifdef MFEM_USE_SINGLE
|
||||
newton_solver.SetRelTol(1e-5);
|
||||
newton_solver.SetAbsTol(1e-8);
|
||||
#else
|
||||
newton_solver.SetRelTol(1e-8);
|
||||
newton_solver.SetAbsTol(1e-12);
|
||||
#endif
|
||||
newton_solver.SetMaxIter(25);
|
||||
newton_solver.SetPrintLevel(0);
|
||||
|
||||
cached_accel = 0.0;
|
||||
}
|
||||
|
||||
real_t HyperbolicHeatOperator::ComputeL2Norm(const Vector &x) const
|
||||
{
|
||||
M_tau->Mult(x, z);
|
||||
return sqrt(InnerProduct(fes.GetComm(), x, z) / tau);
|
||||
}
|
||||
|
||||
void HyperbolicHeatOperator::Mult(const Vector &T, const Vector &dTdt,
|
||||
Vector &d2Tdt2) const
|
||||
{
|
||||
diffusion->Mult(T, rhs);
|
||||
M_rhoc->Mult(dTdt, z);
|
||||
rhs += z;
|
||||
rhs.Neg();
|
||||
rhs += Q;
|
||||
|
||||
rhs.SetSubVector(ess_tdof_list, 0.0);
|
||||
M_solver.Mult(rhs, d2Tdt2);
|
||||
d2Tdt2.SetSubVector(ess_tdof_list, 0.0);
|
||||
|
||||
cached_accel = d2Tdt2;
|
||||
}
|
||||
|
||||
void HyperbolicHeatOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
|
||||
const Vector &T, const Vector &dTdt,
|
||||
Vector &d2Tdt2)
|
||||
{
|
||||
residual_op->SetParameters(fac0, fac1, &T, &dTdt, &Q);
|
||||
|
||||
d2Tdt2 = cached_accel;
|
||||
d2Tdt2.SetSubVector(ess_tdof_list, 0.0);
|
||||
|
||||
Vector zero;
|
||||
newton_solver.Mult(zero, d2Tdt2);
|
||||
MFEM_VERIFY(newton_solver.GetConverged(), "Nonlinear solve failed.");
|
||||
|
||||
cached_accel = d2Tdt2;
|
||||
}
|
||||
|
||||
HyperbolicHeatOperator::~HyperbolicHeatOperator() = default;
|
||||
|
||||
static real_t InitialTemperature(const Vector &x)
|
||||
{
|
||||
real_t r2 = 0.0;
|
||||
for (int d = 0; d < x.Size(); d++)
|
||||
{
|
||||
const real_t s = x(d) - 0.5;
|
||||
r2 += s * s;
|
||||
}
|
||||
return exp(-30.0 * r2);
|
||||
}
|
||||
|
||||
static real_t InitialRate(const Vector &) { return 0.0; }
|
||||
|
||||
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI and HYPRE
|
||||
Mpi::Init(argc, argv);
|
||||
Hypre::Init();
|
||||
|
||||
// 2. Parse command-line options
|
||||
int order = 2;
|
||||
int refinements = 0;
|
||||
const char *device_config = "cpu";
|
||||
int ode_solver_type = 11;
|
||||
real_t t_final = 0.4;
|
||||
real_t dt = 0.02;
|
||||
real_t tau = 1.0;
|
||||
real_t rho_c = 1.0;
|
||||
real_t k0 = 1.0;
|
||||
real_t beta = 1.0;
|
||||
real_t source = 0.0;
|
||||
int prec_type = static_cast<int>(PreconditionerType::Diagonal);
|
||||
const char *assembly_name = "partial";
|
||||
bool dirichlet = true;
|
||||
bool visualization = false;
|
||||
bool paraview = false;
|
||||
int vis_steps = 5;
|
||||
int visport = 19916;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&device_config, "-d", "--device",
|
||||
"Device configuration string, see Device::Configure().");
|
||||
// FEM
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
"Finite element order (polynomial degree).");
|
||||
args.AddOption(&refinements, "-r", "--refinements",
|
||||
"Number of uniform refinements.");
|
||||
args.AddOption(&prec_type, "-pc", "--preconditioner",
|
||||
"Preconditioner for the Newton tangent: 0 = none, "
|
||||
"1 = diagonal/Jacobi.");
|
||||
args.AddOption(&assembly_name, "-al", "--assembly-level",
|
||||
"Assembly level of the two mass operators: partial, legacy, or full.");
|
||||
// Time integration
|
||||
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
|
||||
SecondOrderODESolver::Types.c_str());
|
||||
args.AddOption(&t_final, "-tf", "--t-final",
|
||||
"Final time; start time is 0.");
|
||||
args.AddOption(&dt, "-dt", "--time-step", "Time step.");
|
||||
// Parameters/Physics
|
||||
args.AddOption(&tau, "-tau", "--tau",
|
||||
"Thermal relaxation coefficient of d^2T/dt^2. The wave speed "
|
||||
"is sqrt(k/tau); tau -> 0 recovers the parabolic limit.");
|
||||
args.AddOption(&rho_c, "-rc", "--rho-c",
|
||||
"Volumetric heat capacity rho*c, the damping of the wave.");
|
||||
args.AddOption(&k0, "-k", "--conductivity",
|
||||
"Conductivity k0 in k(T) = k0 (1 + beta T^2).");
|
||||
args.AddOption(&beta, "-b", "--beta",
|
||||
"Nonlinearity beta in k(T) = k0 (1 + beta T^2). "
|
||||
"beta = 0 gives a linear, symmetric problem.");
|
||||
args.AddOption(&source, "-q", "--source", "Constant volumetric source Q.");
|
||||
args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu", "--neumann",
|
||||
"BC switch: T = 0 or zero flux on the whole boundary.");
|
||||
// Visualization
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.AddOption(¶view, "-pv", "--paraview", "-no-pv", "--no-paraview",
|
||||
"Enable or disable ParaView DataCollection output.");
|
||||
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
|
||||
"Visualize every n-th timestep.");
|
||||
args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
|
||||
args.ParseCheck();
|
||||
|
||||
MFEM_VERIFY(tau > 0.0, "-tau must be positive: the hyperbolic term is what "
|
||||
"makes this a second order problem in time.");
|
||||
MFEM_VERIFY(rho_c >= 0.0, "-rc must be non-negative.");
|
||||
MFEM_VERIFY(k0 > 0.0, "-k must be positive.");
|
||||
MFEM_VERIFY(beta >= 0.0, "-b must be non-negative to keep k(T) positive.");
|
||||
|
||||
const AssemblyLevel assembly = ParseAssemblyLevel(assembly_name);
|
||||
|
||||
// 3. Enable hardware devices such as GPUs, and programming models such as
|
||||
// CUDA
|
||||
Device device(device_config);
|
||||
if (Mpi::Root()) { device.Print(); }
|
||||
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "k(T) = " << k0 << " (1 + " << beta << " T^2)"
|
||||
<< ", wave speed sqrt(k0/tau) = " << sqrt(k0 / tau)
|
||||
<< ", relaxation time tau/(rho c) = "
|
||||
<< (rho_c > 0.0 ? tau / rho_c
|
||||
: std::numeric_limits<real_t>::infinity())
|
||||
<< ", tangent: "
|
||||
<< (beta == 0.0 ? "symmetric (CG)" : "non-symmetric (GMRES)")
|
||||
<< ", mass assembly: " << assembly_name << std::endl;
|
||||
}
|
||||
|
||||
// 4. Define the second order ODE solver used for time integration
|
||||
std::unique_ptr<SecondOrderODESolver> ode_solver(
|
||||
SecondOrderODESolver::Select(ode_solver_type));
|
||||
|
||||
// 5. Create the unit square/cube mesh and refine it
|
||||
Mesh mesh = (dim == 2)
|
||||
? Mesh::MakeCartesian2D(8, 8, Element::QUADRILATERAL)
|
||||
: Mesh::MakeCartesian3D(8, 8, 8, Element::HEXAHEDRON);
|
||||
for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
|
||||
mesh.SetCurvature(order);
|
||||
|
||||
// 6. Define a parallel mesh. dFEM reads the mesh coordinates as a field, so
|
||||
// the nodes have to exist.
|
||||
ParMesh pmesh(MPI_COMM_WORLD, mesh);
|
||||
mesh.Clear();
|
||||
pmesh.EnsureNodes();
|
||||
|
||||
// 7. Define a scalar H1 space for the temperature. GlobalTrueVSize is
|
||||
// collective, so every rank has to reach it.
|
||||
H1_FECollection fec(order, dim);
|
||||
ParFiniteElementSpace fes(&pmesh, &fec, 1);
|
||||
const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "Number of temperature unknowns: " << global_dofs
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
// 8. Select quadrature for the nonlinear weak term
|
||||
// k(T) grad(T).grad(w). For k(T) = k0 (1 + beta T^2), its
|
||||
// polynomial degree on affine elements is up to 4*p - 2.
|
||||
const int integration_order =
|
||||
(beta == 0.0) ? 2 * order : 4 * order - 2;
|
||||
const IntegrationRule &ir =
|
||||
IntRules.Get(pmesh.GetTypicalElementGeometry(), integration_order);
|
||||
|
||||
// 9. Set the initial conditions. We assume a gaussian hot spot at rest.
|
||||
ParGridFunction T_gf(&fes), dTdt_gf(&fes);
|
||||
|
||||
FunctionCoefficient T_0(InitialTemperature);
|
||||
T_gf.ProjectCoefficient(T_0);
|
||||
Vector T;
|
||||
T_gf.GetTrueDofs(T);
|
||||
|
||||
FunctionCoefficient dTdt_0(InitialRate);
|
||||
dTdt_gf.ProjectCoefficient(dTdt_0);
|
||||
Vector dTdt;
|
||||
dTdt_gf.GetTrueDofs(dTdt);
|
||||
|
||||
// 10. Initialize the hyperbolic heat operator
|
||||
Array<int> ess_bdr;
|
||||
if (pmesh.bdr_attributes.Size())
|
||||
{
|
||||
ess_bdr.SetSize(pmesh.bdr_attributes.Max());
|
||||
ess_bdr = dirichlet ? 1 : 0;
|
||||
}
|
||||
|
||||
HyperbolicHeatOperator oper(fes, ir, ess_bdr, tau, rho_c, k0, beta, source,
|
||||
static_cast<PreconditionerType>(prec_type), assembly);
|
||||
|
||||
// The projected initial data does not necessarily satisfy T = 0 exactly on
|
||||
// the boundary, and the operator keeps the acceleration at zero there, so the
|
||||
// constraint is imposed on the initial state instead.
|
||||
T.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
|
||||
dTdt.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
|
||||
T_gf.SetFromTrueDofs(T);
|
||||
dTdt_gf.SetFromTrueDofs(dTdt);
|
||||
|
||||
// 11. Set up visualization
|
||||
socketstream sout;
|
||||
if (visualization)
|
||||
{
|
||||
sout.open("localhost", visport);
|
||||
if (!sout)
|
||||
{
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "Unable to connect to GLVis server at localhost:"
|
||||
<< visport << "\nGLVis visualization disabled.\n";
|
||||
}
|
||||
visualization = false;
|
||||
}
|
||||
else
|
||||
{
|
||||
sout.precision(8);
|
||||
sout << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
|
||||
<< "\n";
|
||||
sout << "solution\n" << pmesh << T_gf;
|
||||
sout << "window_title 'Temperature'\n" << std::flush;
|
||||
}
|
||||
}
|
||||
|
||||
// Create a ParaView data collection
|
||||
ParaViewDataCollection pd("dfem-hyperbolic-heat-output", &pmesh);
|
||||
if (paraview)
|
||||
{
|
||||
pd.RegisterField("temperature", &T_gf);
|
||||
pd.RegisterField("rate", &dTdt_gf);
|
||||
pd.SetDataFormat(VTKFormat::BINARY);
|
||||
if (order > 1)
|
||||
{
|
||||
pd.SetHighOrderOutput(true);
|
||||
pd.SetLevelsOfDetail(order);
|
||||
}
|
||||
pd.SetCycle(0);
|
||||
pd.SetTime(0.0);
|
||||
pd.Save();
|
||||
}
|
||||
|
||||
{
|
||||
const real_t T_norm = oper.ComputeL2Norm(T);
|
||||
const real_t dTdt_norm = oper.ComputeL2Norm(dTdt);
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "step t ||T||_L2 ||dT/dt||_L2 newton\n";
|
||||
mfem::out << " 0 " << std::setw(8) << 0.0
|
||||
<< " " << std::setw(11) << T_norm
|
||||
<< " " << std::setw(13) << dTdt_norm
|
||||
<< " " << std::setw(7) << "-" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
// 12. Perform time integration
|
||||
ode_solver->Init(oper);
|
||||
real_t t = 0.0;
|
||||
|
||||
bool last_step = false;
|
||||
for (int ti = 1; !last_step; ti++)
|
||||
{
|
||||
if (t + dt >= t_final - dt / 2) { last_step = true; }
|
||||
|
||||
ode_solver->Step(T, dTdt, t, dt);
|
||||
|
||||
const real_t T_norm = oper.ComputeL2Norm(T);
|
||||
const real_t dTdt_norm = oper.ComputeL2Norm(dTdt);
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << std::setw(4) << ti << " " << std::setw(8) << t
|
||||
<< " " << std::setw(11) << T_norm
|
||||
<< " " << std::setw(13) << dTdt_norm
|
||||
<< " " << std::setw(7) << oper.GetNewtonIterations()
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
if ((last_step || (ti % vis_steps) == 0) && (visualization || paraview))
|
||||
{
|
||||
T_gf.SetFromTrueDofs(T);
|
||||
dTdt_gf.SetFromTrueDofs(dTdt);
|
||||
|
||||
if (visualization)
|
||||
{
|
||||
sout << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
|
||||
<< "\n";
|
||||
sout << "solution\n" << pmesh << T_gf << std::flush;
|
||||
}
|
||||
|
||||
if (paraview)
|
||||
{
|
||||
pd.SetCycle(ti);
|
||||
pd.SetTime(t);
|
||||
pd.Save();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
+110
-98
@@ -18,7 +18,7 @@
|
||||
// Sample runs: mpirun -np 4 dfem-hyperelasticity -o 1 -rs 0 -no-vis
|
||||
// mpirun -np 4 dfem-hyperelasticity -mat linear-elastic -no-vis
|
||||
// mpirun -np 4 dfem-hyperelasticity -mat mooney-rivlin -no-vis
|
||||
// mpirun -np 4 dfem-hyperelasticity -o 2 -rs 1 -pc 2 -no-vis
|
||||
// mpirun -np 4 dfem-hyperelasticity -mat holzapfel -no-vis
|
||||
//
|
||||
// Description: This miniapp solves a quasistatic solid mechanics problem on
|
||||
// the 3D beam used by the Hooke miniapp. The material response is
|
||||
@@ -26,9 +26,6 @@
|
||||
// residual is obtained with DifferentiableOperator::GetDerivative
|
||||
// and Newton's Hessian-vector products are obtained with the new
|
||||
// DifferentiableOperator::GetSecondDerivative functionality.
|
||||
// The same second derivative can also be assembled into a
|
||||
// HypreParMatrix, which -pc 2 uses to build a BoomerAMG
|
||||
// preconditioner for the matrix-free Newton tangent.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "../../fem/dfem/doperator.hpp"
|
||||
@@ -45,6 +42,9 @@ using mfem::future::dual;
|
||||
using dscalar_t = dual<real_t, real_t>;
|
||||
#endif
|
||||
|
||||
#ifdef MFEM_USE_ENZYME
|
||||
|
||||
|
||||
constexpr int dim = 3;
|
||||
constexpr int Displacement = 0;
|
||||
constexpr int Coords = 1;
|
||||
@@ -55,13 +55,13 @@ enum class MaterialType
|
||||
NeoHookean,
|
||||
LinearElastic,
|
||||
MooneyRivlin,
|
||||
Holzapfel
|
||||
};
|
||||
|
||||
enum class PreconditionerType
|
||||
{
|
||||
None,
|
||||
Diagonal,
|
||||
AMG
|
||||
Diagonal
|
||||
};
|
||||
|
||||
MaterialType ParseMaterial(const char *material)
|
||||
@@ -79,8 +79,13 @@ MaterialType ParseMaterial(const char *material)
|
||||
{
|
||||
return MaterialType::MooneyRivlin;
|
||||
}
|
||||
if (name == "holzapfel" || name == "fiber" || name == "fiber-reinforced")
|
||||
{
|
||||
return MaterialType::Holzapfel;
|
||||
}
|
||||
MFEM_ABORT("Unknown material '" << name
|
||||
<< "'. Available materials: neo-hookean, linear-elastic, mooney-rivlin.");
|
||||
<< "'. Available materials: neo-hookean, linear-elastic, "
|
||||
<< "mooney-rivlin, holzapfel.");
|
||||
return MaterialType::NeoHookean;
|
||||
}
|
||||
|
||||
@@ -182,11 +187,61 @@ struct MooneyRivlinEnergy :
|
||||
}
|
||||
};
|
||||
|
||||
template <typename dscalar_t>
|
||||
struct HolzapfelEnergy :
|
||||
HyperelasticEnergyQFunction<HolzapfelEnergy<dscalar_t>, dscalar_t>
|
||||
{
|
||||
real_t c = 50.0;
|
||||
real_t kappa = 100.0;
|
||||
real_t k1 = 10.0;
|
||||
real_t k2 = 20.0;
|
||||
//tensor<real_t, dim> a0 = {1.0, 0.0, 0.0};
|
||||
tensor<real_t, dim> a0 = {0.7071067811865475, 0.7071067811865475, 0.0}; // Normalized fiber direction in the x-y plane at 45 degrees a = (1, 1, 0) / sqrt(2)
|
||||
tensor<real_t, dim, dim> A = make_tensor<dim, dim>(
|
||||
[&](int i, int j) { return a0(i) * a0(j); });
|
||||
|
||||
// Holzapfel-type transversely reinforced model with one fiber family a0:
|
||||
// Ψ(F) = c/2 (Ī₁ - dim) + κ/2 log(J)^2
|
||||
// + k1/(2 k2) (exp(k2 (Ī₄ - 1)^2) - 1),
|
||||
// where Ī₁ = J^(-2/3) tr(C), I₄ = A : C = a0.C.a0
|
||||
// A = a0 ⊗ a0 is the fiber direction tensor.
|
||||
// We use the full I4 invariant, to avoid auxetic behavior
|
||||
//
|
||||
// In this case we assume that the fiber direction is oriented 45 degrees in the x-y plane, i.e. a0 = (1, 1, 0)/sqrt(2).
|
||||
// For a more general case one could start from an external fiber "field", and provide it as an input
|
||||
// to the q-function, and then compute A = a0 ⊗ a0 at each quadrature point.
|
||||
// This would require a different q-function signature including the fiber direction as well.
|
||||
MFEM_HOST_DEVICE inline
|
||||
dscalar_t psi(const tensor<dscalar_t, dim, dim> &F,
|
||||
const tensor<dscalar_t, dim, dim> & /* dudx */) const
|
||||
{
|
||||
// Kinematic quantities
|
||||
const auto C = transpose(F) * F;
|
||||
const auto J = det(F);
|
||||
const auto Jm23 = pow(J, -2.0_r / 3.0_r);
|
||||
|
||||
// Strain invariants
|
||||
const auto I1_bar = Jm23 * tr(C);
|
||||
const auto I4 = ddot(C, A);
|
||||
const auto fiber_strain = I4 - 1.0_r;
|
||||
const auto log_J = log(J);
|
||||
|
||||
// Strain energy density components
|
||||
const auto psi_vol = 0.5_r * kappa * log_J * log_J;
|
||||
const auto psi_iso = 0.5_r * c * (I1_bar - real_t(dim));
|
||||
const auto psi_aniso = (k1 / (2.0_r * k2)) * (exp(k2 * fiber_strain *
|
||||
fiber_strain) - 1.0_r);
|
||||
// NOTE: in practice the anisotropic term should contribute only in tension, i.e. when I4 > 1, or fiber_strain > 0.
|
||||
// mathematically this would introduce a non-smoothness in the energy functional that needs to taken care of.
|
||||
|
||||
return psi_vol + psi_iso + psi_aniso;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
class HyperelasticOperator : public Operator
|
||||
{
|
||||
public:
|
||||
// Matrix-free Hessian-vector product used by Newton's method. This wraps the
|
||||
// functional second-derivative interface and applies the same essential-dof
|
||||
// treatment as the Hooke elasticity Jacobian operator.
|
||||
@@ -238,17 +293,6 @@ public:
|
||||
});
|
||||
}
|
||||
|
||||
// Matrix counterpart of AssembleDiagonal, for preconditioners that need
|
||||
// a real matrix. @a A can be uninitialized; it is allocated by dFEM and
|
||||
// ownership is passed to the caller. Eliminating the essential rows and
|
||||
// columns puts 1.0 on their diagonal, matching Mult above.
|
||||
void AssembleHessian(HypreParMatrix *&A) const
|
||||
{
|
||||
hessian->Assemble(A);
|
||||
auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
|
||||
delete Ae;
|
||||
}
|
||||
|
||||
private:
|
||||
const HyperelasticOperator &oper;
|
||||
Vector state;
|
||||
@@ -266,11 +310,9 @@ public:
|
||||
MaterialType material) :
|
||||
Operator(fes.GetTrueVSize()),
|
||||
fes(fes),
|
||||
ir(ir),
|
||||
qspace(*fes.GetParMesh(), ir),
|
||||
qspace_vec(qspace, 1),
|
||||
q(qspace_vec),
|
||||
material(material)
|
||||
q(qspace_vec)
|
||||
{
|
||||
auto &mesh_nodes =
|
||||
*static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
|
||||
@@ -338,6 +380,19 @@ public:
|
||||
ir, all_domain_attr, derivatives, second_derivatives);
|
||||
break;
|
||||
}
|
||||
case MaterialType::Holzapfel:
|
||||
{
|
||||
// Fiber-reinforced Holzapfel-type material with fibers aligned to
|
||||
// the beam axis. The functional registration lets dFEM derive both
|
||||
// the residual and the Hessian-vector product from the energy.
|
||||
HolzapfelEnergy<dscalar_t> energy;
|
||||
internal_energy_dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
energy,
|
||||
Inputs<Gradient<Displacement>, Gradient<Coords>, Weight> {},
|
||||
Outputs<FunctionalValue<Energy>> {},
|
||||
ir, all_domain_attr, derivatives, second_derivatives);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// The first variation of a functional is exposed as a stateless
|
||||
@@ -384,11 +439,9 @@ public:
|
||||
private:
|
||||
ParFiniteElementSpace &fes;
|
||||
ParFiniteElementSpace *mesh_nodes_fes = nullptr;
|
||||
const IntegrationRule &ir;
|
||||
QuadratureSpace qspace;
|
||||
VectorQuadratureSpace qspace_vec;
|
||||
QuadratureFunction q;
|
||||
MaterialType material;
|
||||
Vector mesh_nodes_tdofs;
|
||||
Array<int> ess_tdofs;
|
||||
Array<int> prescribed_tdofs;
|
||||
@@ -402,78 +455,21 @@ private:
|
||||
mutable std::shared_ptr<HessianOperator> hessian;
|
||||
};
|
||||
|
||||
// BoomerAMG preconditioner for the matrix-free Hessian. HypreBoomerAMG needs a
|
||||
// real HypreParMatrix, so on every SetOperator we let dFEM assemble the second
|
||||
// derivative and setup AMG.
|
||||
class HessianAMG : public Solver
|
||||
{
|
||||
public:
|
||||
HessianAMG(ParFiniteElementSpace *fes_) : fes(fes_)
|
||||
{
|
||||
amg.SetPrintLevel(0);
|
||||
}
|
||||
|
||||
void SetOperator(const Operator &op) override
|
||||
{
|
||||
const auto *H =
|
||||
dynamic_cast<const HyperelasticOperator::HessianOperator *>(&op);
|
||||
MFEM_VERIFY(H, "HessianAMG requires a HessianOperator");
|
||||
height = width = op.Height();
|
||||
|
||||
delete A;
|
||||
A = nullptr;
|
||||
H->AssembleHessian(A);
|
||||
amg.SetOperator(*A);
|
||||
// Tell BoomerAMG this is a dim-component displacement system rather than
|
||||
// a scalar one. order_bynodes MUST be true: the state space is built as
|
||||
// ParFiniteElementSpace(&pmesh, &fec, dim, Ordering::byNODES), while
|
||||
// hypre's default assumes byVDIM. Getting this flag wrong hands hypre a
|
||||
// bogus dof -> function map and converges *worse* than passing no systems
|
||||
// options at all. Must follow SetOperator: the dof map is sized from
|
||||
// height, which SetOperator establishes.
|
||||
amg.SetSystemsOptions(fes->GetVDim(), /*order_bynodes=*/true);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
|
||||
|
||||
~HessianAMG() { delete A; }
|
||||
|
||||
private:
|
||||
ParFiniteElementSpace *fes;
|
||||
HypreParMatrix *A = nullptr;
|
||||
HypreBoomerAMG amg;
|
||||
};
|
||||
|
||||
// Build the preconditioner selected by -pc. Newton hands the current tangent to
|
||||
// the Krylov solver on every iteration, and IterativeSolver::SetOperator
|
||||
// forwards it to the preconditioner, so the AMG hierarchy is rebuilt from the
|
||||
// freshly assembled Hessian at each Newton step.
|
||||
std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type,
|
||||
ParFiniteElementSpace *fes)
|
||||
{
|
||||
switch (type)
|
||||
{
|
||||
case PreconditionerType::None:
|
||||
return nullptr;
|
||||
case PreconditionerType::Diagonal:
|
||||
return std::make_unique<OperatorJacobiSmoother>();
|
||||
case PreconditionerType::AMG:
|
||||
return std::make_unique<HessianAMG>(fes);
|
||||
default:
|
||||
MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
#endif // MFEM_USE_ENZYME
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
Mpi::Init(argc, argv);
|
||||
const int num_procs = Mpi::WorldSize();
|
||||
const int myid = Mpi::WorldRank();
|
||||
Hypre::Init();
|
||||
|
||||
|
||||
#ifndef MFEM_USE_ENZYME
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "This miniapp requires MFEM_USE_ENZYME=YES because it uses "
|
||||
<< "dFEM functional second derivatives.\n";
|
||||
}
|
||||
return 0;
|
||||
#else
|
||||
int order = 1;
|
||||
const char *device_config = "cpu";
|
||||
int serial_refinement_levels = 0;
|
||||
@@ -483,6 +479,7 @@ int main(int argc, char *argv[])
|
||||
bool visualization = true;
|
||||
bool paraview = false;
|
||||
int visport = 19916;
|
||||
const char *outfolder = "./Output";
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
@@ -492,9 +489,9 @@ int main(int argc, char *argv[])
|
||||
args.AddOption(&serial_refinement_levels, "-rs", "--ref-serial",
|
||||
"Number of uniform refinements on the serial mesh.");
|
||||
args.AddOption(&material_name, "-mat", "--material",
|
||||
"Material: neo-hookean, linear-elastic, mooney-rivlin.");
|
||||
"Material: neo-hookean, linear-elastic, mooney-rivlin, or holzapfel.");
|
||||
args.AddOption(&prec_type, "-pc", "--preconditioner",
|
||||
"Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
|
||||
"Preconditioner: 0=none, 1=diagonal.");
|
||||
args.AddOption(&cg_tol, "-tol", "--cg-tol",
|
||||
"Relative tolerance for the CG solver.");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
@@ -504,6 +501,8 @@ int main(int argc, char *argv[])
|
||||
"--no-paraview",
|
||||
"Enable or disable ParaView DataCollection output.");
|
||||
args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
|
||||
args.AddOption(&outfolder, "-of", "--output-folder",
|
||||
"Output folder for ParaView DataCollection files.");
|
||||
args.ParseCheck();
|
||||
|
||||
const MaterialType material = ParseMaterial(material_name);
|
||||
@@ -570,9 +569,18 @@ int main(int argc, char *argv[])
|
||||
cg.SetMaxIter(10000);
|
||||
cg.SetPrintLevel(2);
|
||||
|
||||
std::unique_ptr<Solver> pc =
|
||||
MakePreconditioner(static_cast<PreconditionerType>(prec_type), &fes);
|
||||
if (pc) { cg.SetPreconditioner(*pc); }
|
||||
std::unique_ptr<Solver> pc;
|
||||
switch (static_cast<PreconditionerType>(prec_type))
|
||||
{
|
||||
case PreconditionerType::None:
|
||||
break;
|
||||
case PreconditionerType::Diagonal:
|
||||
pc = std::make_unique<OperatorJacobiSmoother>();
|
||||
cg.SetPreconditioner(*pc);
|
||||
break;
|
||||
default:
|
||||
MFEM_ABORT("Unknown preconditioner type: " << prec_type);
|
||||
}
|
||||
|
||||
NewtonSolver newton(MPI_COMM_WORLD);
|
||||
newton.SetSolver(cg);
|
||||
@@ -596,15 +604,18 @@ int main(int argc, char *argv[])
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
socketstream sol_sock(vishost, visport);
|
||||
sol_sock << "parallel " << num_procs << " " << myid << "\n";
|
||||
sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
|
||||
<< "\n";
|
||||
sol_sock.precision(8);
|
||||
sol_sock << "solution\n" << pmesh << U_gf << std::flush;
|
||||
}
|
||||
|
||||
if (paraview)
|
||||
{
|
||||
// Create a ParaView data collection
|
||||
ParaViewDataCollection pd("dfem-hyperelasticity-output", &pmesh);
|
||||
// Create a ParaView data collection. Save() creates the output directory
|
||||
// tree under the prefix path itself, with the ranks synchronized.
|
||||
ParaViewDataCollection pd("dfem-hyperelasticity", &pmesh);
|
||||
pd.SetPrefixPath(outfolder);
|
||||
pd.RegisterField("displacement", &U_gf);
|
||||
pd.SetDataFormat(VTKFormat::BINARY);
|
||||
if (order > 1)
|
||||
@@ -618,4 +629,5 @@ int main(int argc, char *argv[])
|
||||
}
|
||||
|
||||
return 0;
|
||||
#endif // MFEM_USE_ENZYME
|
||||
}
|
||||
@@ -1,622 +0,0 @@
|
||||
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
|
||||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
//
|
||||
// This file is part of the MFEM library. For more information and source code
|
||||
// availability visit https://mfem.org.
|
||||
//
|
||||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
// CONTRIBUTING.md for details.
|
||||
//
|
||||
// ----------------------------------------------
|
||||
// Nonlinear Poisson (Diffusion) Problem with dFEM
|
||||
// ----------------------------------------------
|
||||
//
|
||||
// Compile with: make dfem-nonlinear-poisson
|
||||
//
|
||||
// Sample runs: mpirun -np 4 dfem-nonlinear-poisson -k grad -f energy -no-vis
|
||||
// mpirun -np 4 dfem-nonlinear-poisson -k grad -f residual -no-vis
|
||||
// mpirun -np 4 dfem-nonlinear-poisson -k sol -f residual -no-vis
|
||||
// mpirun -np 4 dfem-nonlinear-poisson -o 2 -r 4 -pc 2 -no-vis
|
||||
//
|
||||
// Description: This miniapp solves the nonlinear diffusion problem
|
||||
//
|
||||
// -div( kappa(u, grad u) grad u ) = f in Omega
|
||||
// u = 0 on dOmega
|
||||
//
|
||||
// on the unit square (set the compile-time constant dim to 3 for
|
||||
// the unit cube). It shows the two ways a nonlinear
|
||||
// problem can be given to dFEM, and the rule that decides which
|
||||
// one is available: the coefficient (-k) determines whether an
|
||||
// energy exists, and the formulation (-f) picks how it is
|
||||
// written.
|
||||
//
|
||||
// -k grad : kappa = 1 + |grad u|^2 depends only on the gradient,
|
||||
// so the operator is the first variation of
|
||||
// Pi(u) = int psi(|grad u|^2) dx with 2 psi' = kappa.
|
||||
// The tangent is symmetric positive definite, and both
|
||||
// formulations apply.
|
||||
//
|
||||
// -k sol : kappa = 1 + u^2 depends on the solution, so the
|
||||
// operator is NOT the gradient of any functional: the
|
||||
// tangent gains int kappa'(u) du (grad u . grad v) dx,
|
||||
// which is not symmetric. Only -f residual applies, and
|
||||
// the linear solves need GMRES rather than CG.
|
||||
//
|
||||
// -f energy : the q-function returns the energy density psi as
|
||||
// a FunctionalValue. The residual is GetDerivative,
|
||||
// the Newton tangent is GetSecondDerivative.
|
||||
//
|
||||
// -f residual : the q-function returns the pointwise residual, as
|
||||
// in the minimal surface miniapp, and the tangent
|
||||
// is GetDerivative.
|
||||
//
|
||||
// With -k grad both formulations describe the same problem and
|
||||
// produce the same answer. Convergence of both against a
|
||||
// manufactured solution is verified in
|
||||
// tests/unit/dfem/test_nonlinear_poisson.cpp.
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "../../fem/dfem/doperator.hpp"
|
||||
#include "../../fem/dfem/backends/local_qf/prelude.hpp"
|
||||
|
||||
using namespace mfem;
|
||||
|
||||
// This example code demonstrates the use of new features in MFEM that are in
|
||||
// development but exposed through the mfem::future namespace. All features
|
||||
// under this namespace might change their interface or behavior in upcoming
|
||||
// releases until they have stabilized.
|
||||
using namespace mfem::future;
|
||||
using mfem::future::tensor;
|
||||
|
||||
#ifdef MFEM_USE_ENZYME
|
||||
using dscalar_t = real_t;
|
||||
#else
|
||||
using mfem::future::dual;
|
||||
using dscalar_t = dual<real_t, real_t>;
|
||||
#endif
|
||||
|
||||
// Space dimension. The q-functions and the operator are templated on it, so it
|
||||
// is fixed at compile time; set it to 3 to run the same problem on the unit
|
||||
// cube.
|
||||
constexpr int dim = 2;
|
||||
|
||||
// Field IDs shared by every integrator in this miniapp.
|
||||
static constexpr int Solution = 0;
|
||||
static constexpr int Coords = 1;
|
||||
static constexpr int Energy = 2;
|
||||
|
||||
// The coefficient determines whether an energy formulation exists at all.
|
||||
enum class KappaType
|
||||
{
|
||||
GradientDependent, // kappa = 1 + |grad u|^2, symmetric, energy available
|
||||
SolutionDependent // kappa = 1 + u^2, non-symmetric, residual only
|
||||
};
|
||||
|
||||
enum class FormType
|
||||
{
|
||||
Energy, // q-function returns psi, tangent = GetSecondDerivative
|
||||
Residual // q-function returns the residual, tangent = GetDerivative
|
||||
};
|
||||
|
||||
enum class PreconditionerType { None, Diagonal, AMG };
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Pointwise description of the problem: the coefficients and one q-function per
|
||||
// formulation.
|
||||
// ----------------------------------------------------------------------------
|
||||
template <int DIM>
|
||||
struct NonlinearPoisson
|
||||
{
|
||||
// kappa = 1 + |grad u|^2, as a function of s = |grad u|^2.
|
||||
template <typename T>
|
||||
MFEM_HOST_DEVICE static inline T Kappa_grad(const T &s)
|
||||
{
|
||||
return 1.0_r + s;
|
||||
}
|
||||
|
||||
// kappa = 1 + u^2.
|
||||
template <typename T>
|
||||
MFEM_HOST_DEVICE static inline T Kappa_sol(const T &u)
|
||||
{
|
||||
return 1.0_r + u * u;
|
||||
}
|
||||
|
||||
// Energy density with 2 Psi_grad'(s) = Kappa_grad(s), so that the first
|
||||
// variation of int psi dx is int kappa grad u . grad v dx.
|
||||
template <typename T>
|
||||
MFEM_HOST_DEVICE static inline T Psi_grad(const T &s)
|
||||
{
|
||||
return 0.5_r * s + 0.25_r * s * s;
|
||||
}
|
||||
|
||||
// The q-function supplies only the energy density; dFEM forms both the first
|
||||
// variation (the residual) and the second variation (the Newton tangent)
|
||||
// from it. Only Kappa_grad has a potential, so only it can be written here.
|
||||
struct EnergyBased
|
||||
{
|
||||
MFEM_HOST_DEVICE inline
|
||||
auto operator()(const tensor<dscalar_t, DIM> &dudxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
dscalar_t &energy) const
|
||||
{
|
||||
const auto dudx = dudxi * inv(J);
|
||||
// det(J) * w is the quadrature form of dx in Pi(u) = int psi dx.
|
||||
energy = Psi_grad(sqnorm(dudx)) * det(J) * w;
|
||||
}
|
||||
};
|
||||
|
||||
// KAPPA selects the coefficient family. Both variants take the same inputs,
|
||||
// so they register identically; the solution value is interpolated either
|
||||
// way even though only Kappa_sol consumes it.
|
||||
template <KappaType KAPPA>
|
||||
struct ResidualBased
|
||||
{
|
||||
MFEM_HOST_DEVICE inline
|
||||
auto operator()(const dscalar_t &u,
|
||||
const tensor<dscalar_t, DIM> &dudxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
tensor<dscalar_t, DIM> &dvdx) const
|
||||
{
|
||||
const auto invJ = inv(J);
|
||||
const auto dudx = dudxi * invJ;
|
||||
dscalar_t kappa;
|
||||
if constexpr (KAPPA == KappaType::GradientDependent)
|
||||
{
|
||||
kappa = Kappa_grad(sqnorm(dudx));
|
||||
}
|
||||
else
|
||||
{
|
||||
kappa = Kappa_sol(u);
|
||||
}
|
||||
dvdx = kappa * dudx * transpose(invJ) * det(J) * w;
|
||||
}
|
||||
};
|
||||
};
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Depending on @a form, the DifferentiableOperator below holds either the
|
||||
// energy Pi(u), whose first derivative is the residual, or the residual R(u)
|
||||
// directly. Both expose a tangent through the same DerivativeOperator
|
||||
// interface, so the Newton plumbing is shared.
|
||||
// ----------------------------------------------------------------------------
|
||||
template <int DIM>
|
||||
class NonlinearPoissonOperator : public Operator
|
||||
{
|
||||
public:
|
||||
// Matrix-free Newton tangent. Essential directions are removed before the
|
||||
// apply and restored as identity rows afterwards.
|
||||
class JacobianOperator : public Operator
|
||||
{
|
||||
public:
|
||||
JacobianOperator(const NonlinearPoissonOperator &oper,
|
||||
const Vector &state) :
|
||||
Operator(oper.Height()), oper(oper), z(oper.Height())
|
||||
{
|
||||
MultiVector X{state, oper.mesh_nodes_tdofs};
|
||||
// In the energy form the residual is already the first derivative of
|
||||
// the energy, so its tangent is the second derivative. In the residual
|
||||
// form the residual is differentiated directly.
|
||||
tangent = (oper.form == FormType::Energy)
|
||||
? oper.dop->GetSecondDerivative(Solution, X)
|
||||
: oper.dop->GetDerivative(Solution, X);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override
|
||||
{
|
||||
z = x;
|
||||
z.SetSubVector(oper.ess_tdofs, 0.0);
|
||||
|
||||
MultiVector Y{y};
|
||||
tangent->Mult(z, Y);
|
||||
|
||||
auto d_y = y.ReadWrite();
|
||||
const auto d_x = x.Read();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_y[d_dofs[i]] = d_x[d_dofs[i]];
|
||||
});
|
||||
}
|
||||
|
||||
void AssembleDiagonal(Vector &diag) const override
|
||||
{
|
||||
tangent->AssembleDiagonal(diag);
|
||||
auto d_diag = diag.ReadWrite();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_diag[d_dofs[i]] = 1.0;
|
||||
});
|
||||
}
|
||||
|
||||
// Matrix counterpart of AssembleDiagonal, for preconditioners that need a
|
||||
// real matrix. @a A can be uninitialized; it is allocated by dFEM and
|
||||
// ownership is passed to the caller. Eliminating the essential rows and
|
||||
// columns puts 1.0 on their diagonal, matching Mult above.
|
||||
void AssembleJacobian(HypreParMatrix *&A) const
|
||||
{
|
||||
tangent->Assemble(A);
|
||||
auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
|
||||
delete Ae;
|
||||
}
|
||||
|
||||
private:
|
||||
const NonlinearPoissonOperator &oper;
|
||||
mutable Vector z;
|
||||
std::shared_ptr<DerivativeOperator> tangent;
|
||||
};
|
||||
|
||||
NonlinearPoissonOperator(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
FormType form, KappaType kappa) :
|
||||
Operator(fes.GetTrueVSize()),
|
||||
fes(fes),
|
||||
qspace(*fes.GetParMesh(), ir),
|
||||
qspace_vec(qspace, 1),
|
||||
form(form)
|
||||
{
|
||||
MFEM_VERIFY(form == FormType::Residual ||
|
||||
kappa == KappaType::GradientDependent,
|
||||
"the energy formulation requires -k grad: a solution "
|
||||
"dependent coefficient is not the gradient of a functional");
|
||||
|
||||
auto &mesh_nodes =
|
||||
*static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
|
||||
mesh_nodes_fes = mesh_nodes.ParFESpace();
|
||||
mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
|
||||
|
||||
const std::vector<FieldDescriptor> inputs =
|
||||
{
|
||||
{Solution, &fes}, {Coords, mesh_nodes_fes}
|
||||
};
|
||||
|
||||
Array<int> all_domain_attr;
|
||||
if (fes.GetMesh()->attributes.Size() > 0)
|
||||
{
|
||||
all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
|
||||
all_domain_attr = 1;
|
||||
}
|
||||
|
||||
auto derivatives = std::integer_sequence<size_t, Solution> {};
|
||||
|
||||
if (form == FormType::Energy)
|
||||
{
|
||||
// A single scalar per quadrature point, summed by dFEM into the value
|
||||
// of the functional.
|
||||
const std::vector<FieldDescriptor> outputs = {{Energy, &qspace_vec}};
|
||||
dop = std::make_shared<DifferentiableOperator>(
|
||||
inputs, outputs, *fes.GetParMesh());
|
||||
|
||||
// Requesting the second derivative here is what makes the Newton
|
||||
// tangent available later through GetSecondDerivative.
|
||||
auto second_derivatives = SecondDerivatives<Pairs::All> {};
|
||||
typename NonlinearPoisson<DIM>::EnergyBased qf;
|
||||
dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
Inputs<Gradient<Solution>, Gradient<Coords>, Weight> {},
|
||||
Outputs<FunctionalValue<Energy>> {},
|
||||
ir, all_domain_attr, derivatives, second_derivatives);
|
||||
|
||||
// The first variation of a functional is stateless: the current
|
||||
// solution is passed to gradient->Mult(X, Y) on every residual
|
||||
// evaluation, so this wrapper is built once and reused.
|
||||
gradient = dop->GetDerivative(Solution);
|
||||
}
|
||||
else
|
||||
{
|
||||
// The output operator is the gradient of the test function basis,
|
||||
// completing the diffusion-like weak form B^T D(B u, B x, w).
|
||||
const std::vector<FieldDescriptor> outputs = {{Solution, &fes}};
|
||||
dop = std::make_shared<DifferentiableOperator>(
|
||||
inputs, outputs, *fes.GetParMesh());
|
||||
|
||||
const auto in = Inputs < Value<Solution>, Gradient<Solution>,
|
||||
Gradient<Coords>, Weight > {};
|
||||
const auto out = Outputs<Gradient<Solution>> {};
|
||||
|
||||
if (kappa == KappaType::GradientDependent)
|
||||
{
|
||||
typename NonlinearPoisson<DIM>::template
|
||||
ResidualBased<KappaType::GradientDependent> qf;
|
||||
dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
qf, in, out, ir, all_domain_attr, derivatives);
|
||||
}
|
||||
else
|
||||
{
|
||||
typename NonlinearPoisson<DIM>::template
|
||||
ResidualBased<KappaType::SolutionDependent> qf;
|
||||
dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
qf, in, out, ir, all_domain_attr, derivatives);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void SetEssentialAttributes(const Array<int> &ess_bdr)
|
||||
{
|
||||
fes.GetEssentialTrueDofs(ess_bdr, ess_tdofs);
|
||||
}
|
||||
|
||||
// The load stays outside the dFEM functional: the full potential is
|
||||
// Pi(u) - (b, u), whose first variation is the residual below.
|
||||
void SetLoad(const Vector &b) { load = b; }
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override
|
||||
{
|
||||
MultiVector X{x, mesh_nodes_tdofs};
|
||||
MultiVector Y{y};
|
||||
// For a functional, GetDerivative returns the gradient action directly,
|
||||
// including the pointwise reverse seed for the summed energy.
|
||||
if (form == FormType::Energy) { gradient->Mult(X, Y); }
|
||||
else { dop->Mult(X, Y); }
|
||||
|
||||
y -= load;
|
||||
y.SetSubVector(ess_tdofs, 0.0);
|
||||
}
|
||||
|
||||
Operator& GetGradient(const Vector &x) const override
|
||||
{
|
||||
jacobian = std::make_shared<JacobianOperator>(*this, x);
|
||||
return *jacobian;
|
||||
}
|
||||
|
||||
private:
|
||||
ParFiniteElementSpace &fes;
|
||||
ParFiniteElementSpace *mesh_nodes_fes = nullptr;
|
||||
QuadratureSpace qspace;
|
||||
VectorQuadratureSpace qspace_vec;
|
||||
FormType form;
|
||||
Vector mesh_nodes_tdofs, load;
|
||||
Array<int> ess_tdofs;
|
||||
|
||||
std::shared_ptr<DifferentiableOperator> dop;
|
||||
std::shared_ptr<DerivativeOperator> gradient;
|
||||
mutable std::shared_ptr<JacobianOperator> jacobian;
|
||||
};
|
||||
|
||||
// BoomerAMG for the matrix-free tangent. HypreBoomerAMG needs a real
|
||||
// HypreParMatrix, so on every SetOperator we let dFEM assemble the tangent and
|
||||
// set up AMG. This works for both formulations, since GetDerivative and
|
||||
// GetSecondDerivative both return a DerivativeOperator that can assemble itself.
|
||||
template <int DIM>
|
||||
class JacobianAMG : public Solver
|
||||
{
|
||||
using JacobianOperator =
|
||||
typename NonlinearPoissonOperator<DIM>::JacobianOperator;
|
||||
|
||||
public:
|
||||
JacobianAMG() { amg.SetPrintLevel(0); }
|
||||
|
||||
void SetOperator(const Operator &op) override
|
||||
{
|
||||
const auto *Jop = dynamic_cast<const JacobianOperator *>(&op);
|
||||
MFEM_VERIFY(Jop, "JacobianAMG requires a JacobianOperator");
|
||||
height = width = op.Height();
|
||||
|
||||
delete A;
|
||||
A = nullptr;
|
||||
Jop->AssembleJacobian(A);
|
||||
amg.SetOperator(*A);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
|
||||
|
||||
~JacobianAMG() { delete A; }
|
||||
|
||||
private:
|
||||
HypreParMatrix *A = nullptr;
|
||||
HypreBoomerAMG amg;
|
||||
};
|
||||
|
||||
template <int DIM>
|
||||
std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
|
||||
{
|
||||
switch (type)
|
||||
{
|
||||
case PreconditionerType::None: return nullptr;
|
||||
case PreconditionerType::Diagonal: return
|
||||
std::make_unique<OperatorJacobiSmoother>();
|
||||
case PreconditionerType::AMG: return std::make_unique<JacobianAMG<DIM>>();
|
||||
default:
|
||||
MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI and HYPRE
|
||||
Mpi::Init(argc, argv);
|
||||
Hypre::Init();
|
||||
|
||||
// 2. Parse command-line options
|
||||
int order = 1;
|
||||
int refinements = 3;
|
||||
const char *device_config = "cpu";
|
||||
const char *kappa_name = "grad";
|
||||
const char *form_name = "energy";
|
||||
int prec_type = static_cast<int>(PreconditionerType::AMG);
|
||||
real_t krylov_tol = 1e-8;
|
||||
real_t source = 1.0;
|
||||
bool visualization = true;
|
||||
bool paraview = false;
|
||||
int visport = 19916;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
"Finite element order (polynomial degree).");
|
||||
args.AddOption(&refinements, "-r", "--refinements",
|
||||
"Number of uniform refinements.");
|
||||
args.AddOption(&device_config, "-d", "--device",
|
||||
"Device configuration string, see Device::Configure().");
|
||||
args.AddOption(&kappa_name, "-k", "--kappa",
|
||||
"Diffusion coefficient: grad = 1 + |grad u|^2, "
|
||||
"sol = 1 + u^2.");
|
||||
args.AddOption(&form_name, "-f", "--form",
|
||||
"dFEM formulation: energy (requires -k grad), residual.");
|
||||
args.AddOption(&prec_type, "-pc", "--preconditioner",
|
||||
"Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
|
||||
args.AddOption(&krylov_tol, "-tol", "--krylov-tol",
|
||||
"Relative tolerance for the linear solver.");
|
||||
args.AddOption(&source, "-s", "--source", "Constant source term.");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.AddOption(¶view, "-pv", "--paraview", "-no-pv", "--no-paraview",
|
||||
"Enable or disable ParaView DataCollection output.");
|
||||
args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
|
||||
args.ParseCheck();
|
||||
|
||||
const std::string kappa_str(kappa_name), form_str(form_name);
|
||||
MFEM_VERIFY(kappa_str == "grad" || kappa_str == "sol",
|
||||
"-k must be grad or sol");
|
||||
MFEM_VERIFY(form_str == "energy" || form_str == "residual",
|
||||
"-f must be energy or residual");
|
||||
|
||||
const KappaType kappa = (kappa_str == "grad")
|
||||
? KappaType::GradientDependent
|
||||
: KappaType::SolutionDependent;
|
||||
const FormType form = (form_str == "energy") ? FormType::Energy
|
||||
: FormType::Residual;
|
||||
|
||||
// 3. Enable hardware devices such as GPUs, and programming models such as
|
||||
// CUDA
|
||||
Device device(device_config);
|
||||
if (Mpi::Root()) { device.Print(); }
|
||||
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "coefficient: "
|
||||
<< (kappa == KappaType::GradientDependent ? "1 + |grad u|^2"
|
||||
: "1 + u^2")
|
||||
<< ", formulation: " << form_str
|
||||
<< ", tangent: "
|
||||
<< (kappa == KappaType::GradientDependent
|
||||
? "symmetric (CG)" : "non-symmetric (GMRES)")
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
// 4. Create the unit square/cube mesh and refine it
|
||||
Mesh mesh = (dim == 2)
|
||||
? Mesh::MakeCartesian2D(4, 4, Element::QUADRILATERAL)
|
||||
: Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
|
||||
for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
|
||||
mesh.SetCurvature(order);
|
||||
|
||||
// 5. Define a parallel mesh
|
||||
ParMesh pmesh(MPI_COMM_WORLD, mesh);
|
||||
mesh.Clear();
|
||||
pmesh.EnsureNodes();
|
||||
|
||||
// 6. Define a scalar H1 space on the parallel mesh. GlobalTrueVSize is
|
||||
// collective, so every rank has to reach it.
|
||||
H1_FECollection fec(order, dim);
|
||||
ParFiniteElementSpace fes(&pmesh, &fec, 1);
|
||||
const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
|
||||
if (Mpi::Root()) { mfem::out << "#dofs: " << global_dofs << std::endl; }
|
||||
|
||||
// 7. Set up the integration rule. It has to resolve the nonlinear
|
||||
// coefficient, not just the bilinear part, hence the extra order
|
||||
// compared to a linear problem.
|
||||
const IntegrationRule &ir =
|
||||
IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * order + 2);
|
||||
|
||||
// 8. Create the nonlinear operator for the chosen coefficient and
|
||||
// formulation
|
||||
NonlinearPoissonOperator<dim> poisson_op(fes, ir, form, kappa);
|
||||
|
||||
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
poisson_op.SetEssentialAttributes(ess_bdr);
|
||||
|
||||
// 9. Assemble the constant source term. The boundary data is homogeneous.
|
||||
ConstantCoefficient source_coeff(source);
|
||||
ParLinearForm b(&fes);
|
||||
b.AddDomainIntegrator(new DomainLFIntegrator(source_coeff, &ir));
|
||||
b.Assemble();
|
||||
|
||||
Vector B(fes.GetTrueVSize());
|
||||
b.ParallelAssemble(B);
|
||||
poisson_op.SetLoad(B);
|
||||
|
||||
ParGridFunction u_gf(&fes);
|
||||
u_gf = 0.0;
|
||||
|
||||
Vector U(fes.GetTrueVSize());
|
||||
u_gf.GetTrueDofs(U);
|
||||
|
||||
// 10. Set up the linear solver used within Newton's method. The
|
||||
// solution-dependent coefficient gives a non-symmetric tangent, so CG
|
||||
// does not apply there.
|
||||
std::unique_ptr<IterativeSolver> krylov;
|
||||
if (kappa == KappaType::GradientDependent)
|
||||
{
|
||||
krylov = std::make_unique<CGSolver>(MPI_COMM_WORLD);
|
||||
}
|
||||
else
|
||||
{
|
||||
krylov = std::make_unique<GMRESSolver>(MPI_COMM_WORLD);
|
||||
}
|
||||
krylov->SetAbsTol(0.0);
|
||||
krylov->SetRelTol(krylov_tol);
|
||||
krylov->SetMaxIter(2000);
|
||||
krylov->SetPrintLevel(0);
|
||||
|
||||
std::unique_ptr<Solver> pc =
|
||||
MakePreconditioner<dim>(static_cast<PreconditionerType>(prec_type));
|
||||
if (pc) { krylov->SetPreconditioner(*pc); }
|
||||
|
||||
// 11. Set up the nonlinear solver (Newton)
|
||||
NewtonSolver newton(MPI_COMM_WORLD);
|
||||
newton.SetOperator(poisson_op);
|
||||
newton.SetSolver(*krylov);
|
||||
newton.SetAbsTol(0.0);
|
||||
#ifdef MFEM_USE_SINGLE
|
||||
newton.SetRelTol(1e-6);
|
||||
#else
|
||||
newton.SetRelTol(1e-10);
|
||||
#endif
|
||||
newton.SetMaxIter(30);
|
||||
newton.SetPrintLevel(1);
|
||||
|
||||
// 12. Solve the nonlinear system using Newton's method
|
||||
Vector zero;
|
||||
newton.Mult(zero, U);
|
||||
MFEM_VERIFY(newton.GetConverged(), "Newton did not converge");
|
||||
|
||||
u_gf.Distribute(U);
|
||||
|
||||
// 13. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
socketstream sol_sock(vishost, visport);
|
||||
sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
|
||||
<< "\n";
|
||||
sol_sock.precision(8);
|
||||
sol_sock << "solution\n" << pmesh << u_gf << std::flush;
|
||||
}
|
||||
|
||||
// 14. Save the solution in parallel using ParaView format
|
||||
if (paraview)
|
||||
{
|
||||
ParaViewDataCollection pd("dfem-nonlinear-poisson-output", &pmesh);
|
||||
pd.RegisterField("solution", &u_gf);
|
||||
pd.SetDataFormat(VTKFormat::BINARY);
|
||||
if (order > 1)
|
||||
{
|
||||
pd.SetHighOrderOutput(true);
|
||||
pd.SetLevelsOfDetail(order);
|
||||
}
|
||||
pd.SetCycle(0);
|
||||
pd.SetTime(0.0);
|
||||
pd.Save();
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -1,592 +0,0 @@
|
||||
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
|
||||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
//
|
||||
// This file is part of the MFEM library. For more information and source code
|
||||
// availability visit https://mfem.org.
|
||||
//
|
||||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
// CONTRIBUTING.md for details.
|
||||
//
|
||||
// ------------------------------------------------
|
||||
// Nonlinear Reaction-Diffusion with dFEM (two
|
||||
// outputs from a single quadrature point function)
|
||||
// ------------------------------------------------
|
||||
//
|
||||
// Compile with: make dfem-reaction-diffusion
|
||||
//
|
||||
// Sample runs: mpirun -np 4 dfem-reaction-diffusion -b 1 -no-vis
|
||||
// mpirun -np 4 dfem-reaction-diffusion -b 1 -o 2 -r 4 -pc 2 -no-vis
|
||||
// mpirun -np 4 dfem-reaction-diffusion -b 10 -a 0 -no-vis
|
||||
// mpirun -np 4 dfem-reaction-diffusion -b 0 -check -no-vis
|
||||
//
|
||||
// Description: This miniapp solves the nonlinear reaction-diffusion problem
|
||||
//
|
||||
// -div( kappa grad u ) + alpha u + beta u^3 = f in Omega
|
||||
// u = 0 on dOmega
|
||||
//
|
||||
// on the unit square (set the compile-time constant dim to 3 for
|
||||
// the unit cube).
|
||||
//
|
||||
// This miniapp demonstrates the use of dFEM in multiple-output mode.
|
||||
// The q-function writes two outputs onto the same test field, e.g.:
|
||||
//
|
||||
// Outputs<Gradient<Solution>, Value<Solution>>
|
||||
//
|
||||
// dFEM integrates each output against its own test basis
|
||||
// operation and accumulates all of them into a single output
|
||||
// vector,
|
||||
//
|
||||
// R_i = int (kappa grad u) . grad phi_i dx
|
||||
// + int (alpha u + beta u^3 - f) phi_i dx
|
||||
//
|
||||
// so the diffusion and the reaction term share one kernel and one
|
||||
// restriction of u.
|
||||
//
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "../../fem/dfem/doperator.hpp"
|
||||
#include "../../fem/dfem/backends/local_qf/prelude.hpp"
|
||||
|
||||
#include <memory>
|
||||
|
||||
using namespace mfem;
|
||||
|
||||
// This example code demonstrates the use of new features in MFEM that are in
|
||||
// development but exposed through the mfem::future namespace. All features
|
||||
// under this namespace might change their interface or behavior in upcoming
|
||||
// releases until they have stabilized.
|
||||
using namespace mfem::future;
|
||||
using mfem::future::tensor;
|
||||
|
||||
#ifdef MFEM_USE_ENZYME
|
||||
using dscalar_t = real_t;
|
||||
#else
|
||||
using mfem::future::dual;
|
||||
using dscalar_t = dual<real_t, real_t>;
|
||||
#endif
|
||||
|
||||
// Space dimension. The q-function and the operator are templated on it, so it
|
||||
// is fixed at compile time; set it to 3 to run the same problem on the unit
|
||||
// cube.
|
||||
constexpr int dim = 2;
|
||||
|
||||
// Field IDs used by the dFEM integrator.
|
||||
static constexpr int Solution = 0;
|
||||
static constexpr int Coords = 1;
|
||||
|
||||
enum class PreconditionerType { None, Diagonal, AMG };
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Pointwise description of the problem.
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
/// Quadrature point kernel for
|
||||
///
|
||||
/// R(u; v) = int kappa grad u . grad v dx
|
||||
/// + int (alpha u + beta u^3 - f) v dx
|
||||
///
|
||||
/// The two integrands leave through two different outputs: @a dvdx is paired
|
||||
/// with Gradient<Solution> and @a v with Value<Solution>.
|
||||
/// However, both are attached to the same fieldID, so dFEM sums their contributions
|
||||
/// into one residual vector.
|
||||
template <int DIM>
|
||||
struct ReactionDiffusion
|
||||
{
|
||||
real_t kappa = 1.0;
|
||||
real_t alpha = 1.0;
|
||||
real_t beta = 1.0;
|
||||
real_t source = 1.0;
|
||||
|
||||
MFEM_HOST_DEVICE inline
|
||||
auto operator()(const dscalar_t &u,
|
||||
const tensor<dscalar_t, DIM> &dudxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
tensor<dscalar_t, DIM> &dvdx,
|
||||
dscalar_t &v) const
|
||||
{
|
||||
const auto invJ = inv(J);
|
||||
const auto dudx = dudxi * invJ;
|
||||
const real_t dxw = det(J) * w;
|
||||
|
||||
dvdx = kappa * dudx * transpose(invJ) * dxw;
|
||||
v = (alpha * u + beta * u * u * u - source) * dxw;
|
||||
}
|
||||
};
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// The residual operator. A single DifferentiableOperator holds both terms, so
|
||||
// the Newton tangent below is the derivative of the whole q-function.
|
||||
// ----------------------------------------------------------------------------
|
||||
template <int DIM>
|
||||
class ReactionDiffusionOperator : public Operator
|
||||
{
|
||||
public:
|
||||
/// Matrix-free Newton tangent K + (alpha + 3 beta u^2) M. Essential
|
||||
/// directions are removed before the apply and restored as identity rows
|
||||
/// afterwards.
|
||||
class JacobianOperator : public Operator
|
||||
{
|
||||
public:
|
||||
JacobianOperator(const ReactionDiffusionOperator &oper,
|
||||
const Vector &state) :
|
||||
Operator(oper.Height()), oper(oper), z(oper.Height())
|
||||
{
|
||||
MultiVector X{state, oper.mesh_nodes_tdofs};
|
||||
// Differentiating the q-function linearizes both outputs at once,
|
||||
// so no term has to be added by hand here.
|
||||
tangent = oper.dop->GetDerivative(Solution, X);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override
|
||||
{
|
||||
z = x;
|
||||
z.SetSubVector(oper.ess_tdofs, 0.0);
|
||||
|
||||
MultiVector Y{y};
|
||||
tangent->Mult(z, Y);
|
||||
|
||||
auto d_y = y.ReadWrite();
|
||||
const auto d_x = x.Read();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_y[d_dofs[i]] = d_x[d_dofs[i]];
|
||||
});
|
||||
}
|
||||
|
||||
void AssembleDiagonal(Vector &diag) const override
|
||||
{
|
||||
tangent->AssembleDiagonal(diag);
|
||||
auto d_diag = diag.ReadWrite();
|
||||
const auto d_dofs = oper.ess_tdofs.Read();
|
||||
mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
|
||||
{
|
||||
d_diag[d_dofs[i]] = 1.0;
|
||||
});
|
||||
}
|
||||
|
||||
// Matrix counterpart of AssembleDiagonal, for preconditioners that need a
|
||||
// real matrix. @a A can be uninitialized; it is allocated by dFEM and
|
||||
// ownership is passed to the caller. Eliminating the essential rows and
|
||||
// columns puts 1.0 on their diagonal, matching Mult above.
|
||||
void AssembleJacobian(HypreParMatrix *&A) const
|
||||
{
|
||||
tangent->Assemble(A);
|
||||
auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
|
||||
delete Ae;
|
||||
}
|
||||
|
||||
private:
|
||||
const ReactionDiffusionOperator &oper;
|
||||
mutable Vector z;
|
||||
std::shared_ptr<DerivativeOperator> tangent;
|
||||
};
|
||||
|
||||
ReactionDiffusionOperator(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
real_t kappa, real_t alpha, real_t beta,
|
||||
real_t source) :
|
||||
Operator(fes.GetTrueVSize()),
|
||||
fes(fes)
|
||||
{
|
||||
MFEM_VERIFY(kappa > 0.0, "the diffusion coefficient has to be positive");
|
||||
MFEM_VERIFY(alpha >= 0.0 && beta >= 0.0,
|
||||
"negative reaction coefficients make the tangent indefinite; "
|
||||
"this miniapp assumes -a and -b are non-negative so that the "
|
||||
"problem stays convex and CG applies");
|
||||
|
||||
auto &mesh_nodes =
|
||||
*static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
|
||||
mesh_nodes_fes = mesh_nodes.ParFESpace();
|
||||
mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
|
||||
|
||||
const std::vector<FieldDescriptor> inputs =
|
||||
{
|
||||
{Solution, &fes}, {Coords, mesh_nodes_fes}
|
||||
};
|
||||
// One output field, even though the q-function writes two outputs onto
|
||||
// it: the number of FieldOperators does not set the number of blocks.
|
||||
const std::vector<FieldDescriptor> outputs = {{Solution, &fes}};
|
||||
|
||||
Array<int> all_domain_attr;
|
||||
if (fes.GetMesh()->attributes.Size() > 0)
|
||||
{
|
||||
all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
|
||||
all_domain_attr = 1;
|
||||
}
|
||||
|
||||
dop = std::make_shared<DifferentiableOperator>(
|
||||
inputs, outputs, *fes.GetParMesh());
|
||||
|
||||
ReactionDiffusion<DIM> qf{kappa, alpha, beta, source};
|
||||
dop->AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
Inputs < Value<Solution>, Gradient<Solution>,
|
||||
Gradient<Coords>, Weight > {},
|
||||
Outputs<Gradient<Solution>, Value<Solution>> {},
|
||||
ir, all_domain_attr, Derivatives<Solution> {});
|
||||
}
|
||||
|
||||
void SetEssentialAttributes(const Array<int> &ess_bdr)
|
||||
{
|
||||
fes.GetEssentialTrueDofs(ess_bdr, ess_tdofs);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override
|
||||
{
|
||||
MultiVector X{x, mesh_nodes_tdofs};
|
||||
MultiVector Y{y};
|
||||
// Both terms and the load are inside the q-function, so there is nothing
|
||||
// to add to the residual here.
|
||||
dop->Mult(X, Y);
|
||||
y.SetSubVector(ess_tdofs, 0.0);
|
||||
}
|
||||
|
||||
Operator& GetGradient(const Vector &x) const override
|
||||
{
|
||||
jacobian = std::make_shared<JacobianOperator>(*this, x);
|
||||
return *jacobian;
|
||||
}
|
||||
|
||||
private:
|
||||
ParFiniteElementSpace &fes;
|
||||
ParFiniteElementSpace *mesh_nodes_fes = nullptr;
|
||||
Vector mesh_nodes_tdofs;
|
||||
Array<int> ess_tdofs;
|
||||
|
||||
std::shared_ptr<DifferentiableOperator> dop;
|
||||
mutable std::shared_ptr<JacobianOperator> jacobian;
|
||||
};
|
||||
|
||||
// BoomerAMG for the matrix-free tangent. HypreBoomerAMG needs a real
|
||||
// HypreParMatrix, so on every SetOperator we let dFEM assemble the tangent and
|
||||
// set up AMG. The assembled matrix carries both outputs, so the reaction term
|
||||
// is present in the preconditioner as well.
|
||||
template <int DIM>
|
||||
class JacobianAMG : public Solver
|
||||
{
|
||||
using JacobianOperator =
|
||||
typename ReactionDiffusionOperator<DIM>::JacobianOperator;
|
||||
|
||||
public:
|
||||
JacobianAMG() { amg.SetPrintLevel(0); }
|
||||
|
||||
void SetOperator(const Operator &op) override
|
||||
{
|
||||
const auto *Jop = dynamic_cast<const JacobianOperator *>(&op);
|
||||
MFEM_VERIFY(Jop, "JacobianAMG requires a JacobianOperator");
|
||||
height = width = op.Height();
|
||||
|
||||
delete A;
|
||||
A = nullptr;
|
||||
Jop->AssembleJacobian(A);
|
||||
amg.SetOperator(*A);
|
||||
}
|
||||
|
||||
void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
|
||||
|
||||
~JacobianAMG() { delete A; }
|
||||
|
||||
private:
|
||||
HypreParMatrix *A = nullptr;
|
||||
HypreBoomerAMG amg;
|
||||
};
|
||||
|
||||
template <int DIM>
|
||||
std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
|
||||
{
|
||||
switch (type)
|
||||
{
|
||||
case PreconditionerType::None: return nullptr;
|
||||
case PreconditionerType::Diagonal: return
|
||||
std::make_unique<OperatorJacobiSmoother>();
|
||||
case PreconditionerType::AMG: return std::make_unique<JacobianAMG<DIM>>();
|
||||
default:
|
||||
MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Verification of the two-output path against a classical assembly.
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
/// Largest absolute entry of @a v across all ranks.
|
||||
static real_t GlobalNormlinf(MPI_Comm comm, const Vector &v)
|
||||
{
|
||||
real_t local = v.Normlinf(), global = 0.0;
|
||||
MPI_Allreduce(&local, &global, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX,
|
||||
comm);
|
||||
return global;
|
||||
}
|
||||
|
||||
/// With beta = 0 and f = 0 the operator is exactly
|
||||
/// DiffusionIntegrator(kappa) + MassIntegrator(alpha), which a single
|
||||
/// ParBilinearForm can build. Comparing against it checks that the Gradient and
|
||||
/// the Value output really do accumulate into the same residual vector: drop
|
||||
/// either FieldOperator from the Outputs tuple and these comparisons fail.
|
||||
///
|
||||
/// The operator used here is built without essential attributes, so no
|
||||
/// elimination happens on either side and the raw action and matrix are
|
||||
/// compared.
|
||||
template <int DIM>
|
||||
static void CheckLinearLimit(ParFiniteElementSpace &fes,
|
||||
const IntegrationRule &ir,
|
||||
real_t kappa, real_t alpha)
|
||||
{
|
||||
MPI_Comm comm = fes.GetComm();
|
||||
const int tvsize = fes.GetTrueVSize();
|
||||
|
||||
ReactionDiffusionOperator<DIM> lin_op(fes, ir, kappa, alpha, 0.0, 0.0);
|
||||
|
||||
ConstantCoefficient kappa_coeff(kappa), alpha_coeff(alpha);
|
||||
ParBilinearForm a(&fes);
|
||||
a.AddDomainIntegrator(new DiffusionIntegrator(kappa_coeff, &ir));
|
||||
a.AddDomainIntegrator(new MassIntegrator(alpha_coeff, &ir));
|
||||
a.Assemble();
|
||||
a.Finalize();
|
||||
|
||||
Vector U(tvsize);
|
||||
U.Randomize(0x9e3779b9);
|
||||
ParGridFunction u_gf(&fes), y_gf(&fes);
|
||||
u_gf.SetFromTrueDofs(U);
|
||||
|
||||
// 1. Action of the two-output q-function against the two integrators.
|
||||
Vector Y(tvsize), Yref(tvsize);
|
||||
lin_op.Mult(U, Y);
|
||||
a.Mult(u_gf, y_gf);
|
||||
fes.GetProlongationMatrix()->MultTranspose(y_gf, Yref);
|
||||
|
||||
Vector diff(Y);
|
||||
diff -= Yref;
|
||||
const real_t action_err = GlobalNormlinf(comm, diff);
|
||||
|
||||
// 2. Assembled tangent and its diagonal. Both assembly paths walk the
|
||||
// quadrature point cache written by DerivativeSetup, which is laid out
|
||||
// over every output FieldOperator, so both have to pick up the reaction
|
||||
// term as well as the diffusion term. The matrix comparison is done
|
||||
// through matrix-vector products, so that it does not depend on the
|
||||
// sparsity pattern or the ordering of the two assemblies.
|
||||
const auto &J = lin_op.GetGradient(U);
|
||||
const auto *Jop =
|
||||
dynamic_cast<const typename ReactionDiffusionOperator<DIM>::JacobianOperator *>
|
||||
(&J);
|
||||
MFEM_VERIFY(Jop, "expected a JacobianOperator");
|
||||
|
||||
HypreParMatrix *A = nullptr;
|
||||
Jop->AssembleJacobian(A);
|
||||
HypreParMatrix *Aref = a.ParallelAssemble();
|
||||
|
||||
real_t matrix_err = 0.0;
|
||||
Vector v(tvsize), w1(tvsize), w2(tvsize);
|
||||
for (int k = 0; k < 3; k++)
|
||||
{
|
||||
v.Randomize(0x01000193 + k);
|
||||
A->Mult(v, w1);
|
||||
Aref->Mult(v, w2);
|
||||
w1 -= w2;
|
||||
matrix_err = std::max(matrix_err, GlobalNormlinf(comm, w1));
|
||||
}
|
||||
|
||||
Vector d1(tvsize), d2(tvsize);
|
||||
Jop->AssembleDiagonal(d1);
|
||||
Aref->GetDiag(d2);
|
||||
d1 -= d2;
|
||||
const real_t diag_err = GlobalNormlinf(comm, d1);
|
||||
|
||||
delete A;
|
||||
delete Aref;
|
||||
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "linear limit (beta = 0, f = 0) vs "
|
||||
<< "DiffusionIntegrator + MassIntegrator:\n"
|
||||
<< " action max |dFEM - reference| = " << action_err << '\n'
|
||||
<< " tangent max |dFEM - reference| = " << matrix_err << '\n'
|
||||
<< " diagonal max |dFEM - reference| = " << diag_err
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
// Loose enough for the accumulated round-off of two different assembly
|
||||
// paths, tight enough that a missing output term cannot slip through.
|
||||
const real_t tol = 1e-10 * std::max(1.0_r, kappa + alpha);
|
||||
MFEM_VERIFY(action_err < tol, "two-output action does not match the "
|
||||
"DiffusionIntegrator + MassIntegrator reference");
|
||||
MFEM_VERIFY(matrix_err < tol, "assembled two-output tangent does not match "
|
||||
"the DiffusionIntegrator + MassIntegrator reference");
|
||||
MFEM_VERIFY(diag_err < tol, "diagonal of the two-output tangent does not "
|
||||
"match the DiffusionIntegrator + MassIntegrator reference");
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// 1. Initialize MPI and HYPRE
|
||||
Mpi::Init(argc, argv);
|
||||
Hypre::Init();
|
||||
|
||||
// 2. Parse command-line options
|
||||
int order = 1;
|
||||
int refinements = 3;
|
||||
const char *device_config = "cpu";
|
||||
real_t kappa = 1.0;
|
||||
real_t alpha = 1.0;
|
||||
real_t beta = 1.0;
|
||||
real_t source = 1.0;
|
||||
int prec_type = static_cast<int>(PreconditionerType::AMG);
|
||||
real_t krylov_tol = 1e-8;
|
||||
bool check = false;
|
||||
bool visualization = true;
|
||||
bool paraview = false;
|
||||
int visport = 19916;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
"Finite element order (polynomial degree).");
|
||||
args.AddOption(&refinements, "-r", "--refinements",
|
||||
"Number of uniform refinements.");
|
||||
args.AddOption(&device_config, "-d", "--device",
|
||||
"Device configuration string, see Device::Configure().");
|
||||
args.AddOption(&kappa, "-k", "--kappa", "Diffusion coefficient, positive.");
|
||||
args.AddOption(&alpha, "-a", "--alpha",
|
||||
"Linear reaction coefficient, non-negative.");
|
||||
args.AddOption(&beta, "-b", "--beta",
|
||||
"Cubic reaction coefficient, non-negative. Zero gives the "
|
||||
"linear screened Poisson problem.");
|
||||
args.AddOption(&source, "-s", "--source", "Constant source term.");
|
||||
args.AddOption(&prec_type, "-pc", "--preconditioner",
|
||||
"Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
|
||||
args.AddOption(&krylov_tol, "-tol", "--krylov-tol",
|
||||
"Relative tolerance for the linear solver.");
|
||||
args.AddOption(&check, "-check", "--check-linear", "-no-check",
|
||||
"--no-check-linear",
|
||||
"Compare the beta = 0 limit against a ParBilinearForm "
|
||||
"carrying DiffusionIntegrator and MassIntegrator.");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.AddOption(¶view, "-pv", "--paraview", "-no-pv", "--no-paraview",
|
||||
"Enable or disable ParaView DataCollection output.");
|
||||
args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
|
||||
args.ParseCheck();
|
||||
|
||||
// 3. Enable hardware devices such as GPUs, and programming models such as
|
||||
// CUDA
|
||||
Device device(device_config);
|
||||
if (Mpi::Root()) { device.Print(); }
|
||||
|
||||
if (Mpi::Root())
|
||||
{
|
||||
mfem::out << "operator: -div(" << kappa << " grad u) + " << alpha
|
||||
<< " u + " << beta << " u^3 = " << source
|
||||
<< ", outputs: Gradient<Solution>, Value<Solution>"
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
// 4. Create the unit square/cube mesh and refine it
|
||||
Mesh mesh = (dim == 2)
|
||||
? Mesh::MakeCartesian2D(4, 4, Element::QUADRILATERAL)
|
||||
: Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
|
||||
for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
|
||||
mesh.SetCurvature(order);
|
||||
|
||||
// 5. Define a parallel mesh
|
||||
ParMesh pmesh(MPI_COMM_WORLD, mesh);
|
||||
mesh.Clear();
|
||||
pmesh.EnsureNodes();
|
||||
|
||||
// 6. Define a scalar H1 space on the parallel mesh. GlobalTrueVSize is
|
||||
// collective, so every rank has to reach it.
|
||||
H1_FECollection fec(order, dim);
|
||||
ParFiniteElementSpace fes(&pmesh, &fec, 1);
|
||||
const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
|
||||
if (Mpi::Root()) { mfem::out << "#dofs: " << global_dofs << std::endl; }
|
||||
|
||||
// 7. Set up the integration rule. It has to resolve the cubic reaction
|
||||
// term, not just the bilinear part.
|
||||
const IntegrationRule &ir =
|
||||
IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * order + 2);
|
||||
|
||||
// 8. Optionally verify the two-output path against a classical assembly.
|
||||
// This uses its own operator, so it is independent of -b and -s.
|
||||
if (check) { CheckLinearLimit<dim>(fes, ir, kappa, alpha); }
|
||||
|
||||
// 9. Create the nonlinear operator
|
||||
ReactionDiffusionOperator<dim> rd_op(fes, ir, kappa, alpha, beta, source);
|
||||
|
||||
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
rd_op.SetEssentialAttributes(ess_bdr);
|
||||
|
||||
ParGridFunction u_gf(&fes);
|
||||
u_gf = 0.0;
|
||||
|
||||
Vector U(fes.GetTrueVSize());
|
||||
u_gf.GetTrueDofs(U);
|
||||
|
||||
// 10. Set up the linear solver used within Newton's method. For
|
||||
// non-negative alpha and beta the tangent is symmetric positive
|
||||
// definite, so CG works.
|
||||
CGSolver krylov(MPI_COMM_WORLD);
|
||||
krylov.SetAbsTol(0.0);
|
||||
krylov.SetRelTol(krylov_tol);
|
||||
krylov.SetMaxIter(2000);
|
||||
krylov.SetPrintLevel(0);
|
||||
|
||||
std::unique_ptr<Solver> pc =
|
||||
MakePreconditioner<dim>(static_cast<PreconditionerType>(prec_type));
|
||||
if (pc) { krylov.SetPreconditioner(*pc); }
|
||||
|
||||
// 11. Set up the nonlinear solver (Newton)
|
||||
NewtonSolver newton(MPI_COMM_WORLD);
|
||||
newton.SetOperator(rd_op);
|
||||
newton.SetSolver(krylov);
|
||||
newton.SetAbsTol(0.0);
|
||||
#ifdef MFEM_USE_SINGLE
|
||||
newton.SetRelTol(1e-6);
|
||||
#else
|
||||
newton.SetRelTol(1e-10);
|
||||
#endif
|
||||
newton.SetMaxIter(30);
|
||||
newton.SetPrintLevel(1);
|
||||
|
||||
// 12. Solve the nonlinear system using Newton's method
|
||||
Vector zero;
|
||||
newton.Mult(zero, U);
|
||||
MFEM_VERIFY(newton.GetConverged(), "Newton did not converge");
|
||||
|
||||
u_gf.Distribute(U);
|
||||
|
||||
// 13. Send the solution by socket to a GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
char vishost[] = "localhost";
|
||||
socketstream sol_sock(vishost, visport);
|
||||
sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
|
||||
<< "\n";
|
||||
sol_sock.precision(8);
|
||||
sol_sock << "solution\n" << pmesh << u_gf << std::flush;
|
||||
}
|
||||
|
||||
// 14. Save the solution in parallel using ParaView format
|
||||
if (paraview)
|
||||
{
|
||||
ParaViewDataCollection pd("dfem-reaction-diffusion-output", &pmesh);
|
||||
pd.RegisterField("solution", &u_gf);
|
||||
pd.SetDataFormat(VTKFormat::BINARY);
|
||||
if (order > 1)
|
||||
{
|
||||
pd.SetHighOrderOutput(true);
|
||||
pd.SetLevelsOfDetail(order);
|
||||
}
|
||||
pd.SetCycle(0);
|
||||
pd.SetTime(0.0);
|
||||
pd.Save();
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -21,8 +21,7 @@ MFEM_LIB_FILE = mfem_is_not_built
|
||||
-include $(CONFIG_MK)
|
||||
|
||||
SEQ_EXAMPLES =
|
||||
PAR_EXAMPLES = dfem-minimal-surface dfem-hyperelasticity dfem-nonlinear-poisson \
|
||||
dfem-hyperbolic-heat dfem-reaction-diffusion
|
||||
PAR_EXAMPLES = dfem-minimal-surface dfem-hyperelasticity_energy
|
||||
ifeq ($(MFEM_USE_MPI),NO)
|
||||
EXAMPLES = $(SEQ_EXAMPLES)
|
||||
else
|
||||
@@ -65,5 +64,4 @@ clean-build:
|
||||
rm -rf *.dSYM *.TVD.*breakpoints
|
||||
|
||||
clean-exec:
|
||||
@rm -rf dfem-minimal-surface-output dfem-nonlinear-poisson-output \
|
||||
dfem-hyperbolic-heat-output dfem-reaction-diffusion-output Output
|
||||
@rm -rf dfem-minimal-surface-output
|
||||
|
||||
@@ -69,7 +69,9 @@ if (MFEM_USE_BENCHMARK)
|
||||
add_benchmark(assembly_levels)
|
||||
add_benchmark(bpl)
|
||||
add_benchmark(ceed)
|
||||
add_benchmark(dfem)
|
||||
if (MFEM_USE_MPI)
|
||||
add_benchmark(dfem)
|
||||
endif(MFEM_USE_MPI)
|
||||
add_benchmark(dg_amr)
|
||||
add_benchmark(elasticity)
|
||||
add_benchmark(nlvc)
|
||||
|
||||
+416
-576
File diff suppressed because it is too large
Load Diff
@@ -174,56 +174,6 @@ struct mass_diffusion_local_qf
|
||||
}
|
||||
};
|
||||
|
||||
// Three outputs across two test fields: V carries mass + diffusion (two
|
||||
// outputs on one field), P carries a diffusion scaled by kappa. The
|
||||
// two blocks are different, so we should spot if the row blocks
|
||||
// get mismatched/corrupted inadvertedly
|
||||
constexpr real_t kappa = 2.0;
|
||||
|
||||
struct two_field_local_qf
|
||||
{
|
||||
inline MFEM_HOST_DEVICE
|
||||
void operator()(
|
||||
const dscalar_t &u,
|
||||
const tensor<dscalar_t, DIM> &dudxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
dscalar_t &out_v,
|
||||
tensor<dscalar_t, DIM> &out_dv,
|
||||
tensor<dscalar_t, DIM> &out_dp) const
|
||||
{
|
||||
const auto invJ = inv(J);
|
||||
const auto detJ = det(J);
|
||||
const auto grad = (dudxi * invJ) * transpose(invJ);
|
||||
out_v = u * detJ * w;
|
||||
out_dv = grad * (detJ * w);
|
||||
out_dp = grad * (kappa * detJ * w);
|
||||
}
|
||||
};
|
||||
|
||||
// Mass + diffusion on an FE test field, plus quadrature point data on a
|
||||
// VectorQuadratureSpace. The second row block has no basis to contract against,
|
||||
// so it is not assemblable while the first one still is.
|
||||
struct mass_diffusion_qdata_local_qf
|
||||
{
|
||||
inline MFEM_HOST_DEVICE
|
||||
void operator()(
|
||||
const dscalar_t &u,
|
||||
const tensor<dscalar_t, DIM> &dudxi,
|
||||
const tensor<real_t, DIM, DIM> &J,
|
||||
const real_t &w,
|
||||
dscalar_t &out1,
|
||||
tensor<dscalar_t, DIM> &out2,
|
||||
tensor<real_t, DIM, DIM> &out3) const
|
||||
{
|
||||
const auto invJ = inv(J);
|
||||
const auto detJ = det(J);
|
||||
out1 = u * detJ * w;
|
||||
out2 = (dudxi * invJ) * transpose(invJ) * (detJ * w);
|
||||
out3 = J;
|
||||
}
|
||||
};
|
||||
|
||||
TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
|
||||
{
|
||||
const bool all_tests = launch_all_non_regression_tests;
|
||||
@@ -546,7 +496,7 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
|
||||
ParBilinearForm blf_fa(&fes);
|
||||
blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
|
||||
blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
|
||||
blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACY);
|
||||
blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
|
||||
blf_fa.Assemble();
|
||||
blf_fa.Finalize();
|
||||
|
||||
@@ -556,9 +506,6 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
|
||||
{COORDINATES, nodes->ParFESpace()},
|
||||
};
|
||||
|
||||
// The test field is named V even though it is the same space as the
|
||||
// trial field U, which is the usual "V is the test function" idiom.
|
||||
// The diagonal of dR_V/dU is then named explicitly.
|
||||
const std::vector<FieldDescriptor> out_fds
|
||||
{
|
||||
{V, &fes},
|
||||
@@ -595,7 +542,7 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
|
||||
SECTION("Multiple Outputs Assemble Diagonal")
|
||||
{
|
||||
Vector diag(fes.GetTrueVSize()), diag_ref_l(fes.GetVSize());
|
||||
dRdU->AssembleDiagonal(V, diag);
|
||||
dRdU->AssembleDiagonal(diag);
|
||||
blf_fa.SpMat().GetDiag(diag_ref_l);
|
||||
|
||||
Vector diag_ref(fes.GetTrueVSize());
|
||||
@@ -609,312 +556,6 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
}
|
||||
|
||||
// Outputs spanning two test fields. dR/dU is a block column, one row
|
||||
// block per output field, all sharing the trial space of U, and each row
|
||||
// block assembles into its own matrix.
|
||||
// This would look smth like:
|
||||
//
|
||||
// dR/dU = [ dR_P/dU; dR_U/dU ]
|
||||
//
|
||||
{
|
||||
static constexpr int P = 5;
|
||||
|
||||
ParBilinearForm blf_u(&fes);
|
||||
blf_u.AddDomainIntegrator(new MassIntegrator(ir));
|
||||
blf_u.AddDomainIntegrator(new DiffusionIntegrator(ir));
|
||||
blf_u.SetAssemblyLevel(AssemblyLevel::LEGACY);
|
||||
blf_u.Assemble();
|
||||
blf_u.Finalize();
|
||||
|
||||
ConstantCoefficient kappa_coeff(kappa);
|
||||
ParBilinearForm blf_p(&fes);
|
||||
blf_p.AddDomainIntegrator(new DiffusionIntegrator(kappa_coeff, ir));
|
||||
blf_p.SetAssemblyLevel(AssemblyLevel::LEGACY);
|
||||
blf_p.Assemble();
|
||||
blf_p.Finalize();
|
||||
|
||||
const std::vector<FieldDescriptor> in_fds
|
||||
{
|
||||
{U, &fes},
|
||||
{COORDINATES, nodes->ParFESpace()},
|
||||
};
|
||||
|
||||
// U is deliberately *not* the first output field: A[f] and the
|
||||
// blocks of Mult are indexed by position in out_fds, and
|
||||
// AssembleDiagonal has to pick the block whose field is the
|
||||
// differentiated one rather than simply the first.
|
||||
const std::vector<FieldDescriptor> out_fds
|
||||
{
|
||||
{P, &fes},
|
||||
{U, &fes},
|
||||
};
|
||||
|
||||
DifferentiableOperator dop(in_fds, out_fds, pmesh);
|
||||
|
||||
auto qf = two_field_local_qf{};
|
||||
dop.AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
|
||||
tuple{Value<U>{}, Gradient<U>{}, Gradient<P>{}},
|
||||
*ir, all_domain_attr, Derivatives<U> {});
|
||||
|
||||
Vector nodestv;
|
||||
nodes->GetTrueDofs(nodestv);
|
||||
fes.GetRestrictionMatrix()->Mult(x, xtvec);
|
||||
MultiVector X{xtvec, nodestv};
|
||||
auto dRdU = dop.GetDerivative(U, X);
|
||||
|
||||
REQUIRE(dRdU->Height() == 2 * fes.GetTrueVSize());
|
||||
REQUIRE(dRdU->Width() == fes.GetTrueVSize());
|
||||
|
||||
SECTION("Multiple Fields SparseMatrix")
|
||||
{
|
||||
std::vector<SparseMatrix *> A;
|
||||
dRdU->Assemble(A);
|
||||
|
||||
REQUIRE(A.size() == 2);
|
||||
REQUIRE(A[0] != nullptr);
|
||||
REQUIRE(A[1] != nullptr);
|
||||
|
||||
// TestSameMatrices only goes thru the first matrix' sparsity pattern,
|
||||
// so if we compare both ways we can catch entries missing from either side.
|
||||
TestSameMatrices(*A[0], blf_p.SpMat());
|
||||
TestSameMatrices(blf_p.SpMat(), *A[0]);
|
||||
TestSameMatrices(*A[1], blf_u.SpMat());
|
||||
TestSameMatrices(blf_u.SpMat(), *A[1]);
|
||||
|
||||
delete A[0];
|
||||
delete A[1];
|
||||
|
||||
// The element matrix banks are reused, so assembling a second time
|
||||
// (what a Newton loop does every iteration) has to give the same
|
||||
// matrices instead of accumulating on top of the first ones.
|
||||
dRdU->Assemble(A);
|
||||
TestSameMatrices(*A[0], blf_p.SpMat());
|
||||
TestSameMatrices(*A[1], blf_u.SpMat());
|
||||
|
||||
delete A[0];
|
||||
delete A[1];
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
|
||||
SECTION("Multiple Fields HypreParMatrix")
|
||||
{
|
||||
std::vector<HypreParMatrix *> A;
|
||||
dRdU->Assemble(A);
|
||||
|
||||
REQUIRE(A.size() == 2);
|
||||
REQUIRE(A[0] != nullptr);
|
||||
REQUIRE(A[1] != nullptr);
|
||||
|
||||
std::unique_ptr<HypreParMatrix> ref_u(blf_u.ParallelAssemble());
|
||||
std::unique_ptr<HypreParMatrix> ref_p(blf_p.ParallelAssemble());
|
||||
|
||||
// The assembled blocks have to agree with the matrix free action
|
||||
// block by block, which is also what pins A[f] to output field f.
|
||||
Vector dir(fes.GetTrueVSize());
|
||||
dir.Randomize(7);
|
||||
|
||||
Vector yp(fes.GetTrueVSize()), yu(fes.GetTrueVSize());
|
||||
MultiVector Y{yp, yu};
|
||||
dRdU->Mult(dir, Y);
|
||||
|
||||
Vector a(fes.GetTrueVSize());
|
||||
A[0]->Mult(dir, a);
|
||||
a -= Y[0];
|
||||
real_t norm_l = a.Normlinf(), norm_g = norm_l;
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
Vector r(fes.GetTrueVSize());
|
||||
ref_p->Mult(dir, r);
|
||||
a = Y[0];
|
||||
a -= r;
|
||||
norm_l = a.Normlinf();
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
|
||||
A[1]->Mult(dir, a);
|
||||
a -= Y[1];
|
||||
norm_l = a.Normlinf();
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
ref_u->Mult(dir, r);
|
||||
a = Y[1];
|
||||
a -= r;
|
||||
norm_l = a.Normlinf();
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
|
||||
delete A[0];
|
||||
delete A[1];
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
|
||||
SECTION("Multiple Fields Assemble Diagonal")
|
||||
{
|
||||
// Both row blocks are square here, since P and U share a space, so
|
||||
// squareness alone cannot pick one. They carry different operators
|
||||
// though, so reading the wrong row block cannot pass unnoticed.
|
||||
const auto check_diag = [&](const Vector &diag,
|
||||
ParBilinearForm &ref)
|
||||
{
|
||||
REQUIRE(diag.Size() == fes.GetTrueVSize());
|
||||
|
||||
Vector diag_ref_l(fes.GetVSize());
|
||||
ref.SpMat().GetDiag(diag_ref_l);
|
||||
Vector diag_ref(fes.GetTrueVSize());
|
||||
fes.GetProlongationMatrix()->MultTranspose(diag_ref_l, diag_ref);
|
||||
|
||||
Vector d(diag);
|
||||
d -= diag_ref;
|
||||
real_t norm_l = d.Normlinf(), norm_g = norm_l;
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
};
|
||||
|
||||
// Default: the row block of the differentiated field, dR_U/dU.
|
||||
Vector diag_u;
|
||||
dRdU->AssembleDiagonal(diag_u);
|
||||
check_diag(diag_u, blf_u);
|
||||
|
||||
// Named: the other row block, dR_P/dU.
|
||||
Vector diag_p;
|
||||
dRdU->AssembleDiagonal(P, diag_p);
|
||||
check_diag(diag_p, blf_p);
|
||||
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
}
|
||||
|
||||
// The same two field structure, but now the two output fields live on
|
||||
// different spaces. That is what makes AssembleDiagonal's choice of row
|
||||
// block observable at all: with both fields on one space, reading the
|
||||
// wrong block still lands on an identically sized vector holding the
|
||||
// same numbers, so the check above cannot tell the two apart.
|
||||
{
|
||||
static constexpr int P = 5;
|
||||
|
||||
H1_FECollection fec_p(p + 1, DIM);
|
||||
ParFiniteElementSpace fes_p(&pmesh, &fec_p);
|
||||
|
||||
ParBilinearForm blf_u(&fes);
|
||||
blf_u.AddDomainIntegrator(new MassIntegrator(ir));
|
||||
blf_u.AddDomainIntegrator(new DiffusionIntegrator(ir));
|
||||
blf_u.SetAssemblyLevel(AssemblyLevel::LEGACY);
|
||||
blf_u.Assemble();
|
||||
blf_u.Finalize();
|
||||
|
||||
const std::vector<FieldDescriptor> in_fds
|
||||
{
|
||||
{U, &fes},
|
||||
{COORDINATES, nodes->ParFESpace()},
|
||||
};
|
||||
|
||||
// U is the second output field on purpose, and fes_p is a different
|
||||
// space, so picking out_fds[0] would give a differently sized vector.
|
||||
const std::vector<FieldDescriptor> out_fds
|
||||
{
|
||||
{P, &fes_p},
|
||||
{U, &fes},
|
||||
};
|
||||
|
||||
DifferentiableOperator dop(in_fds, out_fds, pmesh);
|
||||
|
||||
auto qf = two_field_local_qf{};
|
||||
dop.AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
|
||||
tuple{Value<U>{}, Gradient<U>{}, Gradient<P>{}},
|
||||
*ir, all_domain_attr, Derivatives<U> {});
|
||||
|
||||
Vector nodestv;
|
||||
nodes->GetTrueDofs(nodestv);
|
||||
fes.GetRestrictionMatrix()->Mult(x, xtvec);
|
||||
MultiVector X{xtvec, nodestv};
|
||||
auto dRdU = dop.GetDerivative(U, X);
|
||||
|
||||
SECTION("Assemble Diagonal Picks The Square Block")
|
||||
{
|
||||
Vector diag;
|
||||
dRdU->AssembleDiagonal(diag);
|
||||
|
||||
// Sized by U's space, not by the first output field's.
|
||||
REQUIRE(diag.Size() == fes.GetTrueVSize());
|
||||
|
||||
Vector diag_ref_l(fes.GetVSize());
|
||||
blf_u.SpMat().GetDiag(diag_ref_l);
|
||||
Vector diag_ref(fes.GetTrueVSize());
|
||||
fes.GetProlongationMatrix()->MultTranspose(diag_ref_l, diag_ref);
|
||||
|
||||
diag -= diag_ref;
|
||||
real_t norm_l = diag.Normlinf(), norm_g = norm_l;
|
||||
MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
|
||||
pmesh.GetComm());
|
||||
REQUIRE(norm_g == MFEM_Approx(0.0));
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
}
|
||||
|
||||
// One assemblable row block and one that is not: S lives on quadrature
|
||||
// points, so it has no basis to contract against and stays null, while
|
||||
// the mass + diffusion block on V assembles as usual.
|
||||
{
|
||||
ParBilinearForm blf_fa(&fes);
|
||||
blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
|
||||
blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
|
||||
blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
|
||||
blf_fa.Assemble();
|
||||
blf_fa.Finalize();
|
||||
|
||||
const std::vector<FieldDescriptor> in_fds
|
||||
{
|
||||
{U, &fes},
|
||||
{COORDINATES, nodes->ParFESpace()},
|
||||
};
|
||||
|
||||
const std::vector<FieldDescriptor> out_fds
|
||||
{
|
||||
{V, &fes},
|
||||
{S, &vqs},
|
||||
};
|
||||
|
||||
DifferentiableOperator dop(in_fds, out_fds, pmesh);
|
||||
|
||||
auto qf = mass_diffusion_qdata_local_qf{};
|
||||
dop.AddDomainIntegrator<LocalQFBackend>(
|
||||
qf,
|
||||
tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
|
||||
tuple{Value<V>{}, Gradient<V>{}, Identity<S>{}},
|
||||
*ir, all_domain_attr, Derivatives<U> {});
|
||||
|
||||
Vector nodestv;
|
||||
nodes->GetTrueDofs(nodestv);
|
||||
fes.GetRestrictionMatrix()->Mult(x, xtvec);
|
||||
MultiVector X{xtvec, nodestv};
|
||||
auto dRdU = dop.GetDerivative(U, X);
|
||||
|
||||
SECTION("Partly Assemblable Outputs SparseMatrix")
|
||||
{
|
||||
std::vector<SparseMatrix *> A;
|
||||
dRdU->Assemble(A);
|
||||
|
||||
REQUIRE(A.size() == 2);
|
||||
REQUIRE(A[0] != nullptr);
|
||||
REQUIRE(A[1] == nullptr);
|
||||
|
||||
TestSameMatrices(*A[0], blf_fa.SpMat());
|
||||
TestSameMatrices(blf_fa.SpMat(), *A[0]);
|
||||
|
||||
delete A[0];
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user