Skip to content

Commit

Permalink
Use new function integration without linear forms
Browse files Browse the repository at this point in the history
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Jul 1, 2024
1 parent 68ea834 commit 1c3dc85
Show file tree
Hide file tree
Showing 8 changed files with 134 additions and 209 deletions.
160 changes: 60 additions & 100 deletions palace/models/lumpedportoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,88 +159,24 @@ double LumpedPortData::GetExcitationVoltage() const
}
}

void LumpedPortData::InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespace) const
{
const auto &mesh = *nd_fespace.GetParMesh();
mfem::Array<int> attr_marker;
if (!s || !v)
{
mfem::Array<int> attr_list;
for (const auto &elem : elems)
{
attr_list.Append(elem->GetAttrList());
}
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mesh::AttrToMarker(bdr_attr_max, attr_list, attr_marker);
}

// The port S-parameter, or the projection of the field onto the port mode, is computed
// as: (E x H_inc) ⋅ n = E ⋅ (E_inc / Z_s), integrated over the port surface.
if (!s)
{
SumVectorCoefficient fb(mesh.SpaceDimension());
for (const auto &elem : elems)
{
const double Rs = R * GetToSquare(*elem);
const double Hinc = (std::abs(Rs) > 0.0)
? 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
elem->GetGeometryLength() * elems.size())
: 0.0;
fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
}
s = std::make_unique<mfem::LinearForm>(&nd_fespace);
s->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
s->UseFastAssembly(false);
s->UseDevice(false);
s->Assemble();
s->UseDevice(true);
}

// The voltage across a port is computed using the electric field solution.
// We have:
// V = ∫ E ⋅ l̂ dl = 1/w ∫ E ⋅ l̂ dS (for rectangular ports)
// or,
// V = 1/(2π) ∫ E ⋅ r̂ / r dS (for coaxial ports).
// We compute the surface integral via an inner product between the linear form with the
// averaging function as a vector coefficient and the solution expansion coefficients.
if (!v)
{
SumVectorCoefficient fb(mesh.SpaceDimension());
for (const auto &elem : elems)
{
fb.AddCoefficient(
elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * elems.size())));
}
v = std::make_unique<mfem::LinearForm>(&nd_fespace);
v->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
v->UseFastAssembly(false);
v->UseDevice(false);
v->Assemble();
v->UseDevice(true);
}
}

std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B) const
{
// Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface using
// the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation (into the
// domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an outward normal,
// so we multiply by -1. The linear form is reconstructed from scratch each time due to
// changing H.
// Compute port power, (E x H⋆) ⋅ n = E ⋅ (-n x H⋆), integrated over the port surface
// using the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation
// (into the domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an
// outward normal, so we multiply by -1.
MFEM_VERIFY((E.HasImag() && B.HasImag()) || (!E.HasImag() && !B.HasImag()),
"Mismatch between real- and complex-valued E and B fields in port power "
"calculation!");
const bool has_imag = E.HasImag();
auto &nd_fespace = *E.ParFESpace();
const auto &mesh = *nd_fespace.GetParMesh();
const auto &mesh = *E.ParFESpace()->GetParMesh();
SumVectorCoefficient fbr(mesh.SpaceDimension()), fbi(mesh.SpaceDimension());
mfem::Array<int> attr_list;
for (const auto &elem : elems)
{
fbr.AddCoefficient(
std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
elem->GetAttrList(), B.Real(), mat_op));
if (has_imag)
if (E.HasImag())
{
fbi.AddCoefficient(
std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
Expand All @@ -250,57 +186,81 @@ std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B)
}
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
std::complex<double> dot;
{
mfem::LinearForm pr(&nd_fespace);
pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbr), attr_marker);
pr.UseFastAssembly(false);
pr.UseDevice(false);
pr.Assemble();
pr.UseDevice(true);
dot = -(pr * E.Real()) + (has_imag ? -1i * (pr * E.Imag()) : 0.0);
}
if (has_imag)
if (E.HasImag())
{
mfem::LinearForm pi(&nd_fespace);
pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi), attr_marker);
pi.UseFastAssembly(false);
pi.UseDevice(false);
pi.Assemble();
pi.UseDevice(true);
dot += -(pi * E.Imag()) + 1i * (pi * E.Real());
Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
BdrInnerProductCoefficient prr(E.Real(), fbr), pir(E.Imag(), fbr), pri(E.Real(), fbi),
pii(E.Imag(), fbi);
mfem::SumCoefficient pr(prr, pii), pi(pir, pri, 1.0, -1.0);
std::complex<double> dot = {-fem::IntegrateFunctionLocal(mesh, attr_marker, true, pr),
-fem::IntegrateFunctionLocal(mesh, attr_marker, true, pi)};
Mpi::GlobalSum(1, &dot, mesh.GetComm());
return dot;
}
else
{
double rdot = dot.real();
Mpi::GlobalSum(1, &rdot, E.ParFESpace()->GetComm());
return rdot;
BdrInnerProductCoefficient pr(E.Real(), fbr);
double dot = -fem::IntegrateFunctionLocal(mesh, attr_marker, true, pr);
Mpi::GlobalSum(1, &dot, mesh.GetComm());
return dot;
}
}

std::complex<double> LumpedPortData::GetSParameter(GridFunction &E) const
{
// Compute port S-parameter, or the projection of the field onto the port mode.
InitializeLinearForms(*E.ParFESpace());
std::complex<double> dot((*s) * E.Real(), 0.0);
// The port S-parameter, or the projection of the field onto the port mode, is computed
// as: (E x H_inc⋆) ⋅ n = E ⋅ (E_inc / Z_s), integrated over the port surface.
const auto &mesh = *E.ParFESpace()->GetParMesh();
SumVectorCoefficient fb(mesh.SpaceDimension());
mfem::Array<int> attr_list;
for (const auto &elem : elems)
{
const double Rs = R * GetToSquare(*elem);
const double Hinc = (std::abs(Rs) > 0.0)
? 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
elem->GetGeometryLength() * elems.size())
: 0.0;
fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
attr_list.Append(elem->GetAttrList());
}
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
BdrInnerProductCoefficient sr(E.Real(), fb);
std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, sr), 0.0);
if (E.HasImag())
{
dot.imag((*s) * E.Imag());
BdrInnerProductCoefficient si(E.Imag(), fb);
dot.imag(fem::IntegrateFunctionLocal(mesh, attr_marker, true, si));
}
Mpi::GlobalSum(1, &dot, E.GetComm());
return dot;
}

std::complex<double> LumpedPortData::GetVoltage(GridFunction &E) const
{
// Compute the average voltage across the port.
InitializeLinearForms(*E.ParFESpace());
std::complex<double> dot((*v) * E.Real(), 0.0);
// The voltage across a port is computed using the electric field solution.
// We have:
// V = ∫ E ⋅ l̂ dl = 1/w ∫ E ⋅ l̂ dS (for rectangular ports)
// or,
// V = 1/(2π) ∫ E ⋅ r̂ / r dS (for coaxial ports).
// We compute the surface integral via an inner product between the linear form with the
// averaging function as a vector coefficient and the solution expansion coefficients.
const auto &mesh = *E.ParFESpace()->GetParMesh();
SumVectorCoefficient fb(mesh.SpaceDimension());
mfem::Array<int> attr_list;
for (const auto &elem : elems)
{
fb.AddCoefficient(
elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * elems.size())));
attr_list.Append(elem->GetAttrList());
}
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
BdrInnerProductCoefficient vr(E.Real(), fb);
std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, vr), 0.0);
if (E.HasImag())
{
dot.imag((*v) * E.Imag());
BdrInnerProductCoefficient vi(E.Imag(), fb);
dot.imag(fem::IntegrateFunctionLocal(mesh, attr_marker, true, vi));
}
Mpi::GlobalSum(1, &dot, E.GetComm());
return dot;
Expand Down
6 changes: 0 additions & 6 deletions palace/models/lumpedportoperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,12 +44,6 @@ class LumpedPortData
double R, L, C;
bool excitation, active;

private:
// Linear forms for postprocessing integrated quantities on the port.
mutable std::unique_ptr<mfem::LinearForm> s, v;

void InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespace) const;

public:
LumpedPortData(const config::LumpedPortData &data, const MaterialOperator &mat_op,
const mfem::ParMesh &mesh);
Expand Down
14 changes: 7 additions & 7 deletions palace/models/postoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ auto CreateParaviewPath(const IoData &iodata, const std::string &name)
PostOperator::PostOperator(const IoData &iodata, SpaceOperator &space_op,
const std::string &name)
: mat_op(space_op.GetMaterialOp()),
surf_post_op(iodata, space_op.GetMaterialOp(), space_op.GetH1Space()),
surf_post_op(iodata, space_op.GetMaterialOp(), space_op.GetMesh()),
dom_post_op(iodata, space_op.GetMaterialOp(), space_op.GetNDSpace(),
space_op.GetRTSpace()),
E(std::make_unique<GridFunction>(space_op.GetNDSpace(),
Expand Down Expand Up @@ -91,15 +91,15 @@ PostOperator::PostOperator(const IoData &iodata, SpaceOperator &space_op,
PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplace_op,
const std::string &name)
: mat_op(laplace_op.GetMaterialOp()),
surf_post_op(iodata, laplace_op.GetMaterialOp(), laplace_op.GetH1Space()),
surf_post_op(iodata, laplace_op.GetMaterialOp(), laplace_op.GetMesh()),
dom_post_op(iodata, laplace_op.GetMaterialOp(), laplace_op.GetH1Space()),
E(std::make_unique<GridFunction>(laplace_op.GetNDSpace())),
V(std::make_unique<GridFunction>(laplace_op.GetH1Space())), lumped_port_init(false),
wave_port_init(false),
paraview(CreateParaviewPath(iodata, name), &laplace_op.GetNDSpace().GetParMesh()),
paraview(CreateParaviewPath(iodata, name), &laplace_op.GetH1Space().GetParMesh()),
paraview_bdr(CreateParaviewPath(iodata, name) + "_boundary",
&laplace_op.GetNDSpace().GetParMesh()),
interp_op(iodata, laplace_op.GetNDSpace().GetParMesh())
&laplace_op.GetH1Space().GetParMesh()),
interp_op(iodata, laplace_op.GetH1Space().GetParMesh())
{
// Note: When using this constructor, you should not use any of the magnetic field related
// postprocessing functions (magnetic field energy, inductor energy, surface currents,
Expand All @@ -121,7 +121,7 @@ PostOperator::PostOperator(const IoData &iodata, LaplaceOperator &laplace_op,
PostOperator::PostOperator(const IoData &iodata, CurlCurlOperator &curlcurl_op,
const std::string &name)
: mat_op(curlcurl_op.GetMaterialOp()),
surf_post_op(iodata, curlcurl_op.GetMaterialOp(), curlcurl_op.GetH1Space()),
surf_post_op(iodata, curlcurl_op.GetMaterialOp(), curlcurl_op.GetMesh()),
dom_post_op(iodata, curlcurl_op.GetMaterialOp(), curlcurl_op.GetNDSpace()),
B(std::make_unique<GridFunction>(curlcurl_op.GetRTSpace())),
A(std::make_unique<GridFunction>(curlcurl_op.GetNDSpace())), lumped_port_init(false),
Expand Down Expand Up @@ -223,7 +223,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
}

// Extract energy density field for electric field energy 1/2 Dᴴ E or magnetic field
// energy 1/2 Hᴴ B. Also Poynting vector S = E x H⋆.
// energy 1/2 Hᴴ B. Also Poynting vector S = Re{E x H⋆}.
if (U_e)
{
paraview.RegisterCoeffField("U_e", U_e.get());
Expand Down
41 changes: 15 additions & 26 deletions palace/models/surfacepostoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,11 +190,10 @@ SurfacePostOperator::InterfaceDielectricData::GetCoefficient(

SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
const MaterialOperator &mat_op,
mfem::ParFiniteElementSpace &h1_fespace)
: mat_op(mat_op), h1_fespace(h1_fespace)
const mfem::ParMesh &mesh)
: mat_op(mat_op)
{
// Check that boundary attributes have been specified correctly.
const auto &mesh = *h1_fespace.GetParMesh();
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> bdr_attr_marker;
if (!iodata.boundaries.postpro.flux.empty() ||
Expand All @@ -219,7 +218,7 @@ SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
data.type == config::SurfaceFluxData::Type::MAGNETIC,
"Electric field or power surface flux postprocessing are not available "
"for electrostatic problems!");
flux_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
flux_surfs.try_emplace(idx, data, mesh, bdr_attr_marker);
}

// Interface dielectric postprocessing.
Expand All @@ -229,7 +228,7 @@ SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
"magnetostatic problems!");
for (const auto &[idx, data] : iodata.boundaries.postpro.dielectric)
{
eps_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
eps_surfs.try_emplace(idx, data, mesh, bdr_attr_marker);
}
}

Expand All @@ -243,14 +242,17 @@ std::complex<double> SurfacePostOperator::GetSurfaceFlux(int idx, const GridFunc
MFEM_VERIFY(it != flux_surfs.end(),
"Unknown surface flux postprocessing index requested!");
const bool has_imag = (E) ? E->HasImag() : B->HasImag();
const auto &mesh = (E) ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
auto f =
it->second.GetCoefficient(E ? &E->Real() : nullptr, B ? &B->Real() : nullptr, mat_op);
std::complex<double> dot(GetLocalSurfaceIntegral(*f, it->second.attr_list), 0.0);
std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f), 0.0);
if (has_imag)
{
f = it->second.GetCoefficient(E ? &E->Imag() : nullptr, B ? &B->Imag() : nullptr,
mat_op);
double doti = GetLocalSurfaceIntegral(*f, it->second.attr_list);
double doti = fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f);
if (it->second.type == SurfaceFluxType::POWER)
{
dot += doti;
Expand All @@ -260,7 +262,7 @@ std::complex<double> SurfacePostOperator::GetSurfaceFlux(int idx, const GridFunc
dot.imag(doti);
}
}
Mpi::GlobalSum(1, &dot, (E) ? E->GetComm() : B->GetComm());
Mpi::GlobalSum(1, &dot, mesh.GetComm());
return dot;
}

Expand All @@ -278,26 +280,13 @@ double SurfacePostOperator::GetInterfaceElectricFieldEnergy(int idx,
auto it = eps_surfs.find(idx);
MFEM_VERIFY(it != eps_surfs.end(),
"Unknown interface dielectric postprocessing index requested!");
const auto &mesh = *E.ParFESpace()->GetParMesh();
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
auto f = it->second.GetCoefficient(E, mat_op);
double dot = GetLocalSurfaceIntegral(*f, it->second.attr_list);
Mpi::GlobalSum(1, &dot, E.GetComm());
double dot = fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f);
Mpi::GlobalSum(1, &dot, mesh.GetComm());
return dot;
}

double SurfacePostOperator::GetLocalSurfaceIntegral(mfem::Coefficient &f,
const mfem::Array<int> &attr_list) const
{
// Integrate the coefficient over the boundary attributes making up this surface index.
const auto &mesh = *h1_fespace.GetParMesh();
int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
mfem::LinearForm s(&h1_fespace);
s.AddBoundaryIntegrator(new BoundaryLFIntegrator(f), attr_marker);
s.UseFastAssembly(false);
s.UseDevice(false);
s.Assemble();
s.UseDevice(true);
return linalg::LocalSum(s);
}

} // namespace palace
9 changes: 1 addition & 8 deletions palace/models/surfacepostoperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,20 +69,13 @@ class SurfacePostOperator
// Reference to material property operator (not owned).
const MaterialOperator &mat_op;

// Reference to scalar finite element space used for computing surface integrals (not
// owned).
mfem::ParFiniteElementSpace &h1_fespace;

double GetLocalSurfaceIntegral(mfem::Coefficient &f,
const mfem::Array<int> &attr_list) const;

public:
// Data structures for postprocessing the surface with the given type.
std::map<int, SurfaceFluxData> flux_surfs;
std::map<int, InterfaceDielectricData> eps_surfs;

SurfacePostOperator(const IoData &iodata, const MaterialOperator &mat_op,
mfem::ParFiniteElementSpace &h1_fespace);
const mfem::ParMesh &mesh);

// Get surface integrals computing electric or magnetic field flux through a boundary.
std::complex<double> GetSurfaceFlux(int idx, const GridFunction *E,
Expand Down
2 changes: 1 addition & 1 deletion palace/models/timeoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ class TimeDependentCurlCurlOperator : public mfem::SecondOrderTimeDependentOpera

// Bindings to SpaceOperator functions to get the system matrix and preconditioner, and
// construct the linear solver.
std::function<void(double a0, double a1)> ConfigureLinearSolver;
std::function<void(double, double)> ConfigureLinearSolver;

public:
TimeDependentCurlCurlOperator(const IoData &iodata, SpaceOperator &space_op,
Expand Down
Loading

0 comments on commit 1c3dc85

Please sign in to comment.