Skip to content

Commit

Permalink
Merge pull request #223 from awslabs/sjg/lumped-port-postpro-fix
Browse files Browse the repository at this point in the history 
Fix current postprocessing for lumped ports
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg authored Apr 4, 2024
2 parents d04debc + 6de1468 commit 25aa677
Show file tree
Hide file tree
Showing 5 changed files with 131 additions and 93 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@ The format of this changelog is based on
- Fixed a bug related to mesh cleaning for unspecified domains and mesh partitioning.
- Change computation of domain energy postprocessing for electrostatic and magnetostatic
simulation types in order to improve performance.
- Fixed a bug when computing the energy associated with lumped elements with more than
one nonzero R, L, or C. This also affects the inductive EPR for lumped inductors with
and associated parallel capacitance.

## [0.12.0] - 2023-12-21

Expand Down
15 changes: 10 additions & 5 deletions palace/models/lumpedportoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,21 +107,26 @@ LumpedPortData::LumpedPortData(const config::LumpedPortData &data,
}
}

std::complex<double> LumpedPortData::GetCharacteristicImpedance(double omega) const
std::complex<double>
LumpedPortData::GetCharacteristicImpedance(double omega,
LumpedPortData::Branch branch) const
{
MFEM_VERIFY((L == 0.0 && C == 0.0) || omega > 0.0,
MFEM_VERIFY((L == 0.0 && C == 0.0) || branch == Branch::R || omega > 0.0,
"Lumped port with nonzero reactance requires frequency in order to define "
"characteristic impedance!");
std::complex<double> Y = 0.0;
if (std::abs(R) > 0.0)
if (std::abs(R) > 0.0 && (branch == Branch::TOTAL || branch == Branch::R))
{
Y += 1.0 / R;
}
if (std::abs(L) > 0.0)
if (std::abs(L) > 0.0 && (branch == Branch::TOTAL || branch == Branch::L))
{
Y += 1.0 / (1i * omega * L);
}
Y += 1i * omega * C;
if (std::abs(C) > 0.0 && (branch == Branch::TOTAL || branch == Branch::C))
{
Y += 1i * omega * C;
}
MFEM_VERIFY(std::abs(Y) > 0.0,
"Characteristic impedance requested for lumped port with zero admittance!")
return 1.0 / Y;
Expand Down
10 changes: 9 additions & 1 deletion palace/models/lumpedportoperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,15 @@ class LumpedPortData
return elem.GetGeometryWidth() / elem.GetGeometryLength() * elems.size();
}

std::complex<double> GetCharacteristicImpedance(double omega = 0.0) const;
enum class Branch
{
TOTAL,
R,
L,
C
};
std::complex<double> GetCharacteristicImpedance(double omega = 0.0,
Branch branch = Branch::TOTAL) const;

double GetExcitationPower() const;
double GetExcitationVoltage() const;
Expand Down
155 changes: 92 additions & 63 deletions palace/models/postoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
#include "fem/errorindicator.hpp"
#include "models/curlcurloperator.hpp"
#include "models/laplaceoperator.hpp"
#include "models/lumpedportoperator.hpp"
#include "models/materialoperator.hpp"
#include "models/spaceoperator.hpp"
#include "models/surfacecurrentoperator.hpp"
Expand Down Expand Up @@ -340,53 +339,6 @@ void PostOperator::SetAGridFunction(const Vector &a, bool exchange_face_nbr_data
}
}

void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega)
{
MFEM_VERIFY(E && B, "Incorrect usage of PostOperator::UpdatePorts!");
if (lumped_port_init)
{
return;
}
for (const auto &[idx, data] : lumped_port_op)
{
auto &vi = lumped_port_vi[idx];
vi.P = data.GetPower(*E, *B);
vi.V = data.GetVoltage(*E);
if (HasImag())
{
MFEM_VERIFY(
omega > 0.0,
"Frequency domain lumped port postprocessing requires nonzero frequency!");
vi.S = data.GetSParameter(*E);
vi.Z = data.GetCharacteristicImpedance(omega);
}
else
{
vi.S = vi.Z = 0.0;
}
}
lumped_port_init = true;
}

void PostOperator::UpdatePorts(const WavePortOperator &wave_port_op, double omega)
{
MFEM_VERIFY(HasImag() && E && B, "Incorrect usage of PostOperator::UpdatePorts!");
if (wave_port_init)
{
return;
}
for (const auto &[idx, data] : wave_port_op)
{
MFEM_VERIFY(omega > 0.0,
"Frequency domain wave port postprocessing requires nonzero frequency!");
auto &vi = wave_port_vi[idx];
vi.P = data.GetPower(*E, *B);
vi.S = data.GetSParameter(*E);
vi.V = vi.Z = 0.0; // Not yet implemented (Z = V² / P, I = V / Z)
}
wave_port_init = true;
}

double PostOperator::GetEFieldEnergy() const
{
if (V)
Expand Down Expand Up @@ -439,16 +391,80 @@ double PostOperator::GetHFieldEnergy(int idx) const
}
}

void PostOperator::UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega)
{
MFEM_VERIFY(E && B, "Incorrect usage of PostOperator::UpdatePorts!");
if (lumped_port_init)
{
return;
}
for (const auto &[idx, data] : lumped_port_op)
{
auto &vi = lumped_port_vi[idx];
vi.P = data.GetPower(*E, *B);
vi.V = data.GetVoltage(*E);
if (HasImag())
{
// Compute current from the port impedance, separate contributions for R, L, C
// branches.
MFEM_VERIFY(
omega > 0.0,
"Frequency domain lumped port postprocessing requires nonzero frequency!");
vi.I[0] =
(std::abs(data.R) > 0.0)
? vi.V / data.GetCharacteristicImpedance(omega, LumpedPortData::Branch::R)
: 0.0;
vi.I[1] =
(std::abs(data.L) > 0.0)
? vi.V / data.GetCharacteristicImpedance(omega, LumpedPortData::Branch::L)
: 0.0;
vi.I[2] =
(std::abs(data.C) > 0.0)
? vi.V / data.GetCharacteristicImpedance(omega, LumpedPortData::Branch::C)
: 0.0;
vi.S = data.GetSParameter(*E);
}
else
{
// Compute current from P = V I⋆ (no scattering parameter output).
vi.I[0] = (std::abs(vi.V) > 0.0) ? std::conj(vi.P / vi.V) : 0.0;
vi.I[1] = vi.I[2] = vi.S = 0.0;
}
}
lumped_port_init = true;
}

void PostOperator::UpdatePorts(const WavePortOperator &wave_port_op, double omega)
{
MFEM_VERIFY(HasImag() && E && B, "Incorrect usage of PostOperator::UpdatePorts!");
if (wave_port_init)
{
return;
}
for (const auto &[idx, data] : wave_port_op)
{
MFEM_VERIFY(omega > 0.0,
"Frequency domain wave port postprocessing requires nonzero frequency!");
auto &vi = wave_port_vi[idx];
vi.P = data.GetPower(*E, *B);
vi.S = data.GetSParameter(*E);
vi.V = vi.I[0] = vi.I[1] = vi.I[2] = 0.0; // Not yet implemented
// (Z = V² / P, I = V / Z)
}
wave_port_init = true;
}

double PostOperator::GetLumpedInductorEnergy(const LumpedPortOperator &lumped_port_op) const
{
// Add contribution due to all capacitive lumped boundaries in the model:
// Add contribution due to all inductive lumped boundaries in the model:
// E_ind = ∑_j 1/2 L_j I_mj².
double U = 0.0;
for (const auto &[idx, data] : lumped_port_op)
{
if (std::abs(data.L) > 0.0)
{
std::complex<double> Ij = GetPortCurrent(lumped_port_op, idx);
std::complex<double> Ij =
GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::L);
U += 0.5 * std::abs(data.L) * std::real(Ij * std::conj(Ij));
}
}
Expand Down Expand Up @@ -555,24 +571,37 @@ std::complex<double> PostOperator::GetPortVoltage(const LumpedPortOperator &lump
return it->second.V;
}

std::complex<double> PostOperator::GetPortCurrent(const LumpedPortOperator &lumped_port_op,
std::complex<double> PostOperator::GetPortVoltage(const WavePortOperator &wave_port_op,
int idx) const
{
MFEM_ABORT("GetPortVoltage is not yet implemented for wave port boundaries!");
return 0.0;
}

std::complex<double> PostOperator::GetPortCurrent(const LumpedPortOperator &lumped_port_op,
int idx,
LumpedPortData::Branch branch) const
{
MFEM_VERIFY(lumped_port_init,
"Lumped port quantities not defined until ports are initialized!");
const auto it = lumped_port_vi.find(idx);
MFEM_VERIFY(it != lumped_port_vi.end(),
"Could not find lumped port when calculating lumped port current!");
if (std::abs(it->second.Z) > 0.0)
{
// Compute from V = I Z when impedance is available.
return it->second.V / it->second.Z;
}
else if (std::abs(it->second.V) > 0.0)
{
// Compute from P = V I⋆.
return std::conj(it->second.P / it->second.V);
}
return ((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::R)
? it->second.I[0]
: 0.0) +
((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::L)
? it->second.I[1]
: 0.0) +
((branch == LumpedPortData::Branch::TOTAL || branch == LumpedPortData::Branch::C)
? it->second.I[2]
: 0.0);
}

std::complex<double> PostOperator::GetPortCurrent(const WavePortOperator &wave_port_op,
int idx) const
{
MFEM_ABORT("GetPortCurrent is not yet implemented for wave port boundaries!");
return 0.0;
}

Expand All @@ -588,7 +617,7 @@ double PostOperator::GetInductorParticipation(const LumpedPortOperator &lumped_p
// An element with no assigned inductance will be treated as having zero admittance and
// thus zero current.
const LumpedPortData &data = lumped_port_op.GetPort(idx);
std::complex<double> Imj = GetPortCurrent(lumped_port_op, idx);
std::complex<double> Imj = GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::L);
return std::copysign(0.5 * std::abs(data.L) * std::real(Imj * std::conj(Imj)) / Em,
Imj.real()); // mean(I²) = (I_r² + I_i²) / 2
}
Expand All @@ -603,7 +632,7 @@ double PostOperator::GetExternalKappa(const LumpedPortOperator &lumped_port_op,
// from which the mode coupling quality factor is computed as:
// Q_mj = ω_m / κ_mj.
const LumpedPortData &data = lumped_port_op.GetPort(idx);
std::complex<double> Imj = GetPortCurrent(lumped_port_op, idx);
std::complex<double> Imj = GetPortCurrent(lumped_port_op, idx, LumpedPortData::Branch::R);
return std::copysign(0.5 * std::abs(data.R) * std::real(Imj * std::conj(Imj)) / Em,
Imj.real()); // mean(I²) = (I_r² + I_i²) / 2
}
Expand Down
41 changes: 17 additions & 24 deletions palace/models/postoperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "linalg/operator.hpp"
#include "linalg/vector.hpp"
#include "models/domainpostoperator.hpp"
#include "models/lumpedportoperator.hpp"
#include "models/surfacepostoperator.hpp"

namespace palace
Expand All @@ -25,7 +26,6 @@ class CurlCurlOperator;
class ErrorIndicator;
class IoData;
class LaplaceOperator;
class LumpedPortOperator;
class MaterialOperator;
class SpaceOperator;
class SurfaceCurrentOperator;
Expand Down Expand Up @@ -60,7 +60,7 @@ class PostOperator
// the grid functions are set.
struct PortPostData
{
std::complex<double> S, P, V, Z;
std::complex<double> P, V, I[3], S;
};
std::map<int, PortPostData> lumped_port_vi, wave_port_vi;
bool lumped_port_init, wave_port_init;
Expand Down Expand Up @@ -96,16 +96,6 @@ class PostOperator
void SetVGridFunction(const Vector &v, bool exchange_face_nbr_data = true);
void SetAGridFunction(const Vector &a, bool exchange_face_nbr_data = true);

// Update cached port voltages and currents for lumped and wave port operators.
void UpdatePorts(const LumpedPortOperator &lumped_port_op,
const WavePortOperator &wave_port_op, double omega = 0.0)
{
UpdatePorts(lumped_port_op, omega);
UpdatePorts(wave_port_op, omega);
}
void UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega = 0.0);
void UpdatePorts(const WavePortOperator &wave_port_op, double omega = 0.0);

// Postprocess the total electric and magnetic field energies in the electric and magnetic
// fields.
double GetEFieldEnergy() const;
Expand All @@ -116,6 +106,16 @@ class PostOperator
double GetEFieldEnergy(int idx) const;
double GetHFieldEnergy(int idx) const;

// Update cached port voltages and currents for lumped and wave port operators.
void UpdatePorts(const LumpedPortOperator &lumped_port_op,
const WavePortOperator &wave_port_op, double omega = 0.0)
{
UpdatePorts(lumped_port_op, omega);
UpdatePorts(wave_port_op, omega);
}
void UpdatePorts(const LumpedPortOperator &lumped_port_op, double omega = 0.0);
void UpdatePorts(const WavePortOperator &wave_port_op, double omega = 0.0);

// Postprocess the energy in lumped capacitor or inductor port boundaries with index in
// the provided set.
double GetLumpedInductorEnergy(const LumpedPortOperator &lumped_port_op) const;
Expand All @@ -136,18 +136,11 @@ class PostOperator
std::complex<double> GetPortPower(const WavePortOperator &wave_port_op, int idx) const;
std::complex<double> GetPortVoltage(const LumpedPortOperator &lumped_port_op,
int idx) const;
std::complex<double> GetPortVoltage(const WavePortOperator &wave_port_op, int idx) const
{
MFEM_ABORT("GetPortVoltage is not yet implemented for wave port boundaries!");
return 0.0;
}
std::complex<double> GetPortCurrent(const LumpedPortOperator &lumped_port_op,
int idx) const;
std::complex<double> GetPortCurrent(const WavePortOperator &wave_port_op, int idx) const
{
MFEM_ABORT("GetPortCurrent is not yet implemented for wave port boundaries!");
return 0.0;
}
std::complex<double> GetPortVoltage(const WavePortOperator &wave_port_op, int idx) const;
std::complex<double>
GetPortCurrent(const LumpedPortOperator &lumped_port_op, int idx,
LumpedPortData::Branch branch = LumpedPortData::Branch::TOTAL) const;
std::complex<double> GetPortCurrent(const WavePortOperator &wave_port_op, int idx) const;

// Postprocess the EPR for the electric field solution and lumped port index.
double GetInductorParticipation(const LumpedPortOperator &lumped_port_op, int idx,
Expand Down

0 comments on commit 25aa677

Please sign in to comment.