Functionalities for flux trapping analysis by dnpham23 · Pull Request #461 · awslabs/palace

dnpham23 · 2025-08-22T00:31:36Z

This PR adds methods to enable flux trapping analysis in magnetostatic simulations. Key additions are:

GetFluxLoopExcitationVector in curlcurloperator.cpp to compute the RHS for flux excitation
2D solver for the surface curl problem in surfacecurlsolver.cpp, along with postprocessing function for verification
Modified magnetostaticsolver.cpp to enable 3D solve using flux BC
Update PostprocessTerminals to enable inductance calculations in flux-only and mix current-flux setups.
Supporting methods in geodata.cpp to find loop boundary edges on submesh and their orientations
Update configfile.cpp to look for flux loop terminals for magnetostatics
Example config files and meshes
MFEM patch

hughcars

I haven't managed to review the MFEM portion directly yet, but this should be a lot to be going on with. At a high level,

check running on more than 2 ranks, I was getting errors for 3 or more from geom factors not being initialized etc.
run with order >= 3, I was getting errors where the computed flux was not matching the target flux, not sure what is causing this directly.

On the interface, broadly think more about how to a) avoid any 3D object allocation within the surface pde solving function and b) pass such objects in as buffers. There are already existing A sized gridfunctions, and even a true sized A buffer that can be used to compute the RHS (The next solution vector has already been allocated, you can use that!)

simlapointe · 2025-09-25T18:59:25Z

In light of the recently merged #459, let's not forget to save all .msh files in binary format (just add gmsh.option.setNumber("Mesh.Binary", 1) to the mesh-generating julia scripts and regenerate the meshes)

dnpham23 · 2025-10-01T21:12:04Z

@hughcars All your requested changes are addressed. Sorry for the delay. Please have a look and feel free to take over if you feel it is more convenient and efficient that way.

hughcars

Thanks, this is looking close. I left a few small inline suggestions for the remaining cleanup/fixes I found while building and running the current PR head through Spack.

The easiest path is probably to cherry-pick the two fix-up commits from:

hughcars/flux-trapping-dev

Those commits are:

11b8ba92f Fix flux loop units and rebase issues
285b2fc3b Add flux loop regression coverage

They include the inline-suggested changes, plus a few things that are not easy to express as GitHub suggestions because they touch unchanged lines or add new files:

config validation for unsupported mixed SurfaceCurrent + FluxLoop use,
duplicate Index validation including FluxLoop,
C++/Catch2 regression coverage for the circular-hole flux-loop example and reference outputs,
SchemaVer bump to 1-2-0.

If applying fixes manually instead of cherry-picking, please also make sure the schema version is bumped as a REVISION bump (1-2-0), since this PR extends the accepted config model.

hughcars · 2026-07-14T18:47:03Z

+  MFEM_VERIFY(n_step > 0, "No surface current boundaries or flux loops specified for "
+                          "magnetostatic simulation!");
+  MFEM_VERIFY(


Add the solver-side backstop for the unsupported mixed case as well as config validation. This keeps accidental programmatic construction from reaching the mixed postprocessing path.

Suggested change

MFEM_VERIFY(n_step > 0, "No surface current boundaries or flux loops specified for "

"magnetostatic simulation!");

MFEM_VERIFY(

MFEM_VERIFY(n_step > 0, "No surface current boundaries or flux loops specified for "

"magnetostatic simulation!");

MFEM_VERIFY(n_current_steps == 0 || n_flux_steps == 0,

"Combining SurfaceCurrent and FluxLoop excitations in the same "

"magnetostatic simulation is not yet supported!");

MFEM_VERIFY(

hughcars · 2026-07-14T18:47:03Z

+  if (n_current > 0 && n_flux > 0)
+  {
+    // Mixed current-flux case: use constraint system M×R = I
+    mfem::DenseMatrix R(n);
+
+    // Diagonal terms from energy
+    for (int i = 0; i < n; i++)
+    {
+      if (is_flux_loop[i])
+        R(i, i) = cross_energy(i, i) / (Phi_inc[i] * Phi_inc[i]);
+      else
+        M(i, i) = cross_energy(i, i) / (I_inc[i] * I_inc[i]);
+    }
+
+    int n_off = n * (n - 1) / 2;  // Number of off-diagonal elements
+    mfem::DenseMatrix A_sys(n * n, 2 * n_off);
+    mfem::Vector b_sys(n * n), x_sol(2 * n_off);
+    A_sys = 0.0;
+    b_sys = 0.0;
+
+    // Set up M×R = I constraint equations
+    for (int i = 0; i < n; i++)
+    {
+      for (int j = 0; j < n; j++)
+      {
+        int eq = i * n + j;
+        b_sys[eq] = (i == j) ? 1.0 : 0.0;
+
+        for (int k = 0; k < n; k++)
+        {
+          if (i != k && k != j)  // Off-diagonal terms
+          {
+            int M_idx = (i < k) ? i * n + k - (i + 1) * (i + 2) / 2
+                                : k * n + i - (k + 1) * (k + 2) / 2;
+            int R_idx = (k < j) ? k * n + j - (k + 1) * (k + 2) / 2
+                                : j * n + k - (j + 1) * (j + 2) / 2;
+            A_sys(eq, M_idx) += (k == j) ? 1.0 : 0.0;
+            A_sys(eq, n_off + R_idx) += (i == k) ? 1.0 : 0.0;
+          }
+        }
+        // Diagonal contributions
+        if (i == j)
+          b_sys[eq] -= M(i, i) * R(j, j);
+      }


For this PR, I would delete the mixed SurfaceCurrent/FluxLoop reconstruction. The current MR = I loop builds a singular system, and the correct mixed formulation needs a separate hybrid/Legendre treatment.

Suggested change

if (n_current > 0 && n_flux > 0)

{

// Mixed current-flux case: use constraint system M×R = I

mfem::DenseMatrix R(n);

// Diagonal terms from energy

for (int i = 0; i < n; i++)

{

if (is_flux_loop[i])

R(i, i) = cross_energy(i, i) / (Phi_inc[i] * Phi_inc[i]);

else

M(i, i) = cross_energy(i, i) / (I_inc[i] * I_inc[i]);

}

int n_off = n * (n - 1) / 2; // Number of off-diagonal elements

mfem::DenseMatrix A_sys(n * n, 2 * n_off);

mfem::Vector b_sys(n * n), x_sol(2 * n_off);

A_sys = 0.0;

b_sys = 0.0;

// Set up M×R = I constraint equations

for (int i = 0; i < n; i++)

{

for (int j = 0; j < n; j++)

{

int eq = i * n + j;

b_sys[eq] = (i == j) ? 1.0 : 0.0;

for (int k = 0; k < n; k++)

{

if (i != k && k != j) // Off-diagonal terms

{

int M_idx = (i < k) ? i * n + k - (i + 1) * (i + 2) / 2

: k * n + i - (k + 1) * (k + 2) / 2;

int R_idx = (k < j) ? k * n + j - (k + 1) * (k + 2) / 2

: j * n + k - (j + 1) * (j + 2) / 2;

A_sys(eq, M_idx) += (k == j) ? 1.0 : 0.0;

A_sys(eq, n_off + R_idx) += (i == k) ? 1.0 : 0.0;

}

}

// Diagonal contributions

if (i == j)

b_sys[eq] -= M(i, i) * R(j, j);

}

hughcars · 2026-07-14T18:47:03Z

+    // Solve system
+    mfem::DenseMatrixInverse A_inv(A_sys);
+    A_inv.Mult(b_sys, x_sol);
+
+    // Extract off-diagonal elements
+    int idx = 0;
+    for (int i = 0; i < n; i++)


Delete the remainder of the mixed solve block.

Suggested change

// Solve system

mfem::DenseMatrixInverse A_inv(A_sys);

A_inv.Mult(b_sys, x_sol);

// Extract off-diagonal elements

int idx = 0;

for (int i = 0; i < n; i++)

hughcars · 2026-07-14T18:47:03Z

+      for (int j = i + 1; j < n; j++)
+      {
+        M(i, j) = M(j, i) = x_sol[idx];
+        R(i, j) = R(j, i) = x_sol[n_off + idx];
+        idx++;
+      }
+    }
+
+    // Compute Minv = R
+    Minv = R;
+  }


Delete the remainder of the mixed solve block.

Suggested change

for (int j = i + 1; j < n; j++)

{

M(i, j) = M(j, i) = x_sol[idx];

R(i, j) = R(j, i) = x_sol[n_off + idx];

idx++;

}

}

// Compute Minv = R

Minv = R;

}

hughcars · 2026-07-14T18:47:03Z

+    // Compute Minv = R
+    Minv = R;
+  }
+  else if (n_flux == n)


After deleting the mixed branch, the pure-flux path becomes the first branch.

Suggested change

else if (n_flux == n)

if (n_flux == n)

hughcars · 2026-07-14T18:47:03Z

+  if constexpr (true)
+  {
+    // Use Palace's KSP solver (production approach)
+    X = 0.0;
+    ksp.Mult(RHS, X);
+
+    // Set solution directly in MFEM GridFunction
+    A.SetFromTrueDofs(X);
+  }
+  else
+  {
+    // Alternative debugging approach using direct MFEM solver
+    // Set up and solve system
+    mfem::ParBilinearForm a(&nd_fespace_submesh);
+
+    // Add curl term: ∫ (curl A) · (curl v) dΩ
+    mfem::ConstantCoefficient curl_reg(1.0);
+    a.AddDomainIntegrator(new mfem::CurlCurlIntegrator(curl_reg));
+
+    // Add small regularization for stability
+    mfem::ConstantCoefficient reg_param(1e-6);
+    a.AddDomainIntegrator(new mfem::VectorFEMassIntegrator(reg_param));
+
+    nd_fespace_submesh.GetEssentialTrueDofs(combined_inner_bdr_marker,
+                                            submesh_ess_tdof_list);
+    a.Assemble();
+    mfem::ParLinearForm b(&nd_fespace_submesh);
+    b.Assemble();
+
+    mfem::HypreParMatrix A_mat;
+    Vector B, X;
+    a.FormLinearSystem(submesh_ess_tdof_list, A, b, A_mat, X, B);
+
+    mfem::GMRESSolver gmres(MPI_COMM_WORLD);
+    gmres.SetOperator(A_mat);
+    gmres.SetRelTol(1e-8);
+    gmres.SetMaxIter(1000);
+    gmres.SetPrintLevel(0);
+    gmres.Mult(B, X);
+
+    a.RecoverFEMSolution(X, b, A);
+  }


I would delete the unreachable debug/direct-solver branch before merging and keep only the production KSP path.

Suggested change

if constexpr (true)

{

// Use Palace's KSP solver (production approach)

X = 0.0;

ksp.Mult(RHS, X);

// Set solution directly in MFEM GridFunction

A.SetFromTrueDofs(X);

}

else

{

// Alternative debugging approach using direct MFEM solver

// Set up and solve system

mfem::ParBilinearForm a(&nd_fespace_submesh);

// Add curl term: ∫ (curl A) · (curl v) dΩ

mfem::ConstantCoefficient curl_reg(1.0);

a.AddDomainIntegrator(new mfem::CurlCurlIntegrator(curl_reg));

// Add small regularization for stability

mfem::ConstantCoefficient reg_param(1e-6);

a.AddDomainIntegrator(new mfem::VectorFEMassIntegrator(reg_param));

nd_fespace_submesh.GetEssentialTrueDofs(combined_inner_bdr_marker,

submesh_ess_tdof_list);

a.Assemble();

mfem::ParLinearForm b(&nd_fespace_submesh);

b.Assemble();

mfem::HypreParMatrix A_mat;

Vector B, X;

a.FormLinearSystem(submesh_ess_tdof_list, A, b, A_mat, X, B);

mfem::GMRESSolver gmres(MPI_COMM_WORLD);

gmres.SetOperator(A_mat);

gmres.SetRelTol(1e-8);

gmres.SetMaxIter(1000);

gmres.SetPrintLevel(0);

gmres.Mult(B, X);

a.RecoverFEMSolution(X, b, A);

}

// Use Palace's KSP solver (production approach)

X = 0.0;

ksp.Mult(RHS, X);

// Set solution directly in MFEM GridFunction

A.SetFromTrueDofs(X);

hughcars · 2026-07-14T18:47:04Z

+The first configuration (`circular_hole.json`) models a circular metal plate of radius
+``R = 3\,\mu\text{m}`` with a concentric hole of radius ``r = 1\,\mu\text{m}``. One
+flux quantum ``\Phi_0 = 2.068 \times 10^{-15}\,\text{Wb}`` is prescribed through the hole.


Separate the internal solve normalization from the physical one-Φ₀ interpretation.

Suggested change

The first configuration (`circular_hole.json`) models a circular metal plate of radius

``R = 3\,\mu\text{m}`` with a concentric hole of radius ``r = 1\,\mu\text{m}``. One

flux quantum ``\Phi_0 = 2.068 \times 10^{-15}\,\text{Wb}`` is prescribed through the hole.

The first configuration (`circular_hole.json`) models a circular metal plate of radius

``R = 3\,\mu\text{m}`` with a concentric hole of radius ``r = 1\,\mu\text{m}``. A unit

nondimensional flux-loop excitation amplitude is prescribed through the hole.

hughcars · 2026-07-14T18:47:04Z

+    problem is solved.
+  - `"HoleAttributes"`: boundary attributes of the hole perimeters where integral
+    constraints are applied.
+  - `"FluxAmounts"`: prescribed flux through each hole, in units of ``\Phi_0``.


This should match the internal-amplitude semantics used by the solver.

Suggested change

- `"FluxAmounts"`: prescribed flux through each hole, in units of ``\Phi_0``.

- `"FluxAmounts"`: prescribed nondimensional flux-loop excitation amplitudes through

each hole.

hughcars · 2026-07-14T18:47:04Z

+For the single-hole configuration, the solver extracts a self-inductance of
+``M = 2.808\,\text{pH}`` for one flux quantum trapped in the hole, with a stored magnetic
+energy of ``E_\text{mag} = 1.680 \times 10^{-13}\,\text{J}``.


The extracted inductance is independent of the internal amplitude normalization. If we want to quote the physical energy for one trapped flux quantum, compute it from the extracted M afterward.

Suggested change

For the single-hole configuration, the solver extracts a self-inductance of

``M = 2.808\,\text{pH}`` for one flux quantum trapped in the hole, with a stored magnetic

energy of ``E_\text{mag} = 1.680 \times 10^{-13}\,\text{J}``.

For the single-hole configuration, the solver extracts a self-inductance of

``M = 1.902\,\text{pH}``. This corresponds to a stored magnetic energy of

``E_\text{mag} = Φ₀^2 / (2M) = 1.124 \times 10^{-18}\,\text{J}`` for one physical flux

quantum ``Φ₀ = 2.0678 \times 10^{-15}\,\text{Wb}`` trapped in the hole.

hughcars · 2026-07-14T18:47:04Z

+When independent flux excitations are configured, the solver extracts the full inductance
+matrix from the stored magnetic energy:
+
+```math
+M_{ij} = \frac{\mathbf{A}_j^T K \mathbf{A}_i}{\Phi_i \Phi_j},
+```
+
+where ``K`` is the curl-curl stiffness matrix. The diagonal entries ``M_{ii}`` give the


M_ij is still the inductance entry in Φ = M I; the issue is what the energy normalization computes directly for flux-controlled solves. It gives reluctance first, then M is recovered by inversion.

Suggested change

When independent flux excitations are configured, the solver extracts the full inductance

matrix from the stored magnetic energy:

```math

M_{ij} = \frac{\mathbf{A}_j^T K \mathbf{A}_i}{\Phi_i \Phi_j},

```

where ``K`` is the curl-curl stiffness matrix. The diagonal entries ``M_{ii}`` give the

When independent flux excitations are configured, the solver extracts the full inductance

matrix by first computing the reluctance matrix from the stored magnetic energy:

```math

R_{ij} = \frac{\mathbf{A}_j^T K \mathbf{A}_i}{\Phi_i \Phi_j}, \qquad M = R^{-1},

```

where ``K`` is the curl-curl stiffness matrix, ``R`` is the reluctance matrix, and ``M`` is

the inductance matrix written to `terminal-M.csv`. The diagonal entries ``M_{ii}`` give the

…ssing

…xample

…the meshes used.

…ation

This reverts commit d702d6c.

…ank runs

…rate reference data

…atch

hughcars force-pushed the main branch from fb2f1e6 to 1a4175e Compare August 25, 2025 13:03

hughcars requested changes Aug 25, 2025

View reviewed changes

hughcars reviewed Sep 24, 2025

View reviewed changes

Comment thread examples/double_circular_hole/mesh/two_square_sheets.jl Outdated

dnpham23 marked this pull request as ready for review October 1, 2025 21:12

dnpham23 force-pushed the main branch 3 times, most recently from 0d0272c to c7581bf Compare June 24, 2026 07:38

dnpham23 requested a review from hughcars June 29, 2026 21:45

dnpham23 force-pushed the main branch from 7e90093 to cedb55b Compare June 30, 2026 20:18

hughcars force-pushed the main branch from 437a8ff to 6569376 Compare July 1, 2026 23:09

dnpham23 force-pushed the main branch from 6569376 to 15de023 Compare July 9, 2026 00:40

hughcars reviewed Jul 14, 2026

View reviewed changes

dnpham23 added 16 commits July 14, 2026 14:32

Add MFEM patch for boundary edge DOF extraction and ldof synchronization

7ef7cdf

Add FluxLoopData config parsing for flux loop boundaries

999c939

Add FluxLoop JSON schema definition for boundary configuration

7b466fc

Add documentation for flux loop analysis and boundaries

1daa607

Add SurfaceFluxOperator for flux loop boundary management

973502d

Add SurfaceCurlSolver for computing flux loop excitation vectors

508aba2

Extend CurlCurlOperator with flux loop support and boundary properties

3abfd2b

Add ref_dir parameter to flux coefficient for magnetic flux direction

65ada8b

Add geodata utilities for submesh boundary edge orientation

a962429

Update magnetostatic solver for mixed current-flux terminal postproce…

aa48537

…ssing

Add example configurations and meshes for flux loop simulations

5177473

formatting fixes

ec1dcbe

Reorganize examples and mesh files

aa0daad

Update mesh file for multiple flux terminals on disjointed geometry e…

889b1f1

…xample

WIP - documentation for flux trapping analysis examples

37db26c

Update documentation and examples for flux trapping analysis

5157888

dnpham23 and others added 8 commits July 14, 2026 14:32

Update .md file explaning flux trapping examples and results

d0c425c

Fix CI: bump schema version, add MFEM patch to spack build

719982b

Fix default output filename in sheet_w_two_holes mesh script

7599776

Rename MFEM patch to mfem_pr4983.diff to match upstream PR convention

2a2b458

Fix CI: clang-format-19 style and broken doc links

cec1048

update the results on M matrix in the two-hole example to match with …

b38e3a3

…the meshes used.

Fix flux loop units and rebase issues

7fc4793

Add flux loop regression coverage

f324706

dnpham23 force-pushed the main branch from 6aa82c8 to f324706 Compare July 14, 2026 21:34

dnpham23 added 7 commits July 14, 2026 15:48

Fix GPU regression: disable device vectors in surface curl solver

0c6f910

Fix multi-rank segfault: call ExchangeFaceNbrData before flux verific…

ce92f9e

…ation

Update regression reference data for 32-rank parallel run

d702d6c

Revert "Update regression reference data for 32-rank parallel run"

3f37e1a

This reverts commit d702d6c.

Fix parallel boundary marker sizing and flux verification for multi-r…

a719e6f

…ank runs

Update circular hole mesh for robust parallel partitioning and regene…

91960a3

…rate reference data

Restore UseDevice for submesh vectors to fix GPU operator-vector mism…

89e5ec2

…atch

dnpham23 force-pushed the main branch from 85f9289 to 89e5ec2 Compare July 15, 2026 07:28

	- `"FluxAmounts"`: prescribed flux through each hole, in units of ``\Phi_0``.
	- `"FluxAmounts"`: prescribed nondimensional flux-loop excitation amplitudes through
	each hole.

Uh oh!

Conversation

dnpham23 commented Aug 22, 2025

Uh oh!

hughcars left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

simlapointe commented Sep 25, 2025

Uh oh!

dnpham23 commented Oct 1, 2025

Uh oh!

hughcars left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants