Skip to content

Commit 6fd02b2

Browse files
AlejandroMunozManterolaAlejandroMunozManterola
committed
Working on tensor-style PML / Volumetric simple PML implemented
1 parent fd84cb6 commit 6fd02b2

21 files changed

Lines changed: 4670 additions & 35 deletions

‎PML.md‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
# Volumetric PML in dgtd
2+
3+
This document describes the first volumetric PML implementation in `dgtd`: what is implemented, how the main pieces interact, and what a user must provide to build a valid PML case.
4+
5+
## Scope
6+
7+
The current implementation is intentionally narrow:
8+
9+
- Only the `GlobalOperator` path is supported.
10+
- The PML region must be an explicitly meshed outer shell.
11+
- The shell must be vacuum-matched.
12+
- The shell must be an axis-aligned wrapper around the physical region.
13+
- The first formulation is sigma-only.
14+
- `kappa = 1` and `alpha = 0` are kept internal.
15+
- Serial and MPI execution are supported.
16+
17+
Not implemented in this first version:
18+
19+
- `HesthavenEvolution`
20+
- `MaxwellEvolution`
21+
- arbitrary PML shapes
22+
- automatic shell generation
23+
- user-exposed `kappa`, `alpha`, or `sigma_max`
24+
25+
## What Was Implemented
26+
27+
### Model and parsing
28+
29+
- `driver.cpp` recognizes `model.materials[]` entries with `"type": "PML"`.
30+
- `Material` stays a simple scalar material holder. PML data is stored separately instead of creating a `PML_Material` subclass.
31+
- `Model` stores PML metadata in `GeomTagToPMLRegion` and provides geometry helpers through `PMLBoxGeometry`.
32+
- `Model::inferPMLBoxGeometry()` checks that the tagged region is a valid outer shell and computes inner and outer extents plus thickness per axis.
33+
34+
### Submesh and runtime state
35+
36+
- `VolumetricPMLSubMesher` builds the PML-only submesh using `SubMesh::CreateFromDomain` in serial and `ParSubMesh::CreateFromDomain` in MPI.
37+
- `PMLWrapper` owns the PML submesh, the PML finite element space, the parent/submesh vdof mapping, the graded sigma profile, and the local correction operator.
38+
- `PMLRegionState` stores the PML auxiliary state used across time steps.
39+
40+
### Operator path
41+
42+
- `DGOperatorFactory` can build a matched conductive operator on the PML submesh.
43+
- The full-domain `buildGlobalOperator()` still provides the main Maxwell update.
44+
- The PML operator is not a replacement for the global operator; it is a local correction applied only on the PML region.
45+
46+
### Time stepping
47+
48+
- `GlobalEvolution` builds and owns `PMLWrapper` when the model contains PML volumes.
49+
- `Solver::step()` mirrors the SGBC pattern and now calls a PML checkpoint/finalize pair around the main ODE step.
50+
- After the usual explicit update, `PMLWrapper` gathers the PML dofs from the parent fields, solves a local implicit correction on the PML submesh, and scatters the corrected values back into the parent fields.
51+
52+
### Diagnostics and protection
53+
54+
- Startup validation is handled in `driver.cpp`.
55+
- The solver prints inferred thickness per active axis.
56+
- When the source spectrum can be estimated, the solver also prints maximum frequency, shortest wavelength, and thickness-to-wavelength ratios.
57+
- The solver prints normal PML cell size, cells across the thickness, effective dofs across the thickness, and the internal `sigma_max`.
58+
- Invalid geometry hard-fails.
59+
- Clearly underresolved shells hard-fail.
60+
- Marginal shells warn.
61+
- Very thick shells only report an informational message.
62+
63+
## How The Pieces Interact
64+
65+
The flow is:
66+
67+
1. `driver.cpp` parses a material entry tagged as `PML`.
68+
2. `Model` stores the PML region metadata and infers the shell geometry.
69+
3. `GlobalEvolution` creates a `PMLWrapper` when PML volumes are present.
70+
4. `PMLWrapper` builds the PML submesh, submesh finite element space, sigma profile, and local matched conductive operator.
71+
5. `Solver::step()` runs the usual explicit global step.
72+
6. Right after that step, `PMLWrapper` applies an implicit correction only on the PML dofs.
73+
7. The corrected PML values are written back into the parent field vectors.
74+
75+
In short: the main solver still advances the full Maxwell system as before, and the PML is added as a local post-step correction on a dedicated submesh.
76+
77+
## What The User Must Provide
78+
79+
To build a valid PML case, the user must provide all of the following:
80+
81+
- A mesh where the physical region and the PML shell are separate domain tags.
82+
- A shell that wraps the physical region as an outer box.
83+
- A shell that is axis-aligned.
84+
- A shell tagged in `model.materials[]` with `"type": "PML"`.
85+
- A physical region that remains a normal material entry, usually vacuum.
86+
87+
The user must not provide:
88+
89+
- manual PML thickness in JSON
90+
- manual `sigma_max`
91+
- manual `kappa` or `alpha`
92+
- non-vacuum constitutive values on the PML material entry
93+
94+
The current JSON contract is based on `model.materials[]`, not on boundary conditions. A minimal pattern is:
95+
96+
```json
97+
{
98+
"model": {
99+
"materials": [
100+
{
101+
"tags": [1],
102+
"type": "vacuum"
103+
},
104+
{
105+
"tags": [2],
106+
"type": "PML"
107+
}
108+
]
109+
}
110+
}
111+
```
112+
113+
Here `1` is the physical region and `2` is the outer shell.
114+
115+
## Expectations For Mesh Design
116+
117+
When preparing a case, the user should expect these rules to matter:
118+
119+
- The shell must be the outermost region.
120+
- The physical region must remain inside the shell with a clean box-like interface.
121+
- The shell must exist in the mesh already; the solver will not create it.
122+
- The shell should have enough cells across the normal thickness for the chosen polynomial order.
123+
- If the source definition includes enough information to estimate frequency content, the solver will also judge thickness relative to wavelength.
124+
125+
Practical guidance:
126+
127+
- Use more than the bare minimum cell count across the shell.
128+
- Keep the shell thickness consistent on opposite sides when possible.
129+
- Start with vacuum physical material plus vacuum-matched PML shell.
130+
- Prefer source descriptions with `spread` and, when relevant, `frequency`, so the wavelength-based diagnostics are available.
131+
132+
## Current Validation Surface
133+
134+
The implementation is covered by focused tests rather than a full production PML case mesh:
135+
136+
- parser and model tests for `type: PML`
137+
- geometry inference tests
138+
- serial and MPI submesh extraction tests
139+
- direct `PMLWrapper` construction and damping tests
140+
- a 1D end-to-end solver comparison showing stronger damping with the PML shell than with the same shell treated as ordinary vacuum
141+
142+
This means the core path is implemented and exercised, but a full external case still depends on the user supplying a purpose-built volumetric PML mesh.
143+
144+
## Summary
145+
146+
The current volumetric PML support is a real runtime feature, not only parser scaffolding. The solver can now:
147+
148+
- recognize PML regions from `model.materials[]`
149+
- verify that the tagged region is a valid outer shell
150+
- build a PML-only submesh
151+
- assemble a local matched conductive correction
152+
- apply an implicit PML update after the explicit global step
153+
- protect the user with startup diagnostics before time stepping begins
154+
155+
To use it successfully, the user must provide a valid outer-shell mesh and tag that shell as a PML material region.

‎README.md‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,8 @@ An OpenSEMBA/dgtd JSON file must have the following structure. Legend: **[REQUIR
136136
- Integer. Number of uniform refinement levels to apply to the loaded mesh. Optional.
137137
- **materials**:
138138
- Array. At least one entry is required. Each entry defines the electromagnetic properties of one or more mesh domains:
139+
- type:
140+
- String. Optional material class selector. Omit it, or use `"vacuum"`, for the normal scalar-material path. Use `"PML"` to tag an explicitly meshed volumetric PML shell.
139141
- **tags**:
140142
- Array of integers. Mesh volume/surface/segment attribute IDs that share these material properties. (Volumes in 3D, Surfaces in 2D, Segments in 1D.)
141143
- relative_permittivity:
@@ -144,6 +146,10 @@ An OpenSEMBA/dgtd JSON file must have the following structure. Legend: **[REQUIR
144146
- Double. Relative permeability $\mu_r$. Defaults to `1.0`.
145147
- bulk_conductivity:
146148
- Double. Bulk electrical conductivity in S/m. Defaults to `0.0`. Internally scaled by the free-space impedance.
149+
- *(For `type: "PML"` only)*
150+
- First-version volumetric PML support is limited to the GlobalOperator path, vacuum-matched regions, and an axis-aligned outer shell explicitly present in the mesh.
151+
- `relative_permittivity`, `relative_permeability`, and `bulk_conductivity` must be omitted. Thickness is inferred from the mesh geometry; `sigma_max`, grading, and validation thresholds remain internal in v1.
152+
- At startup the solver prints inferred shell thickness, estimated maximum source frequency when available, shortest wavelength, normal cell size, cells and effective DOFs across thickness, and the internally chosen `sigma_max`. Invalid or clearly underresolved shells hard-fail before time stepping.
147153
- **boundaries**:
148154
- Array. At least one entry is required. Each entry defines a boundary or interface condition:
149155
- **tags**:
@@ -213,6 +219,10 @@ An OpenSEMBA/dgtd JSON file must have the following structure. Legend: **[REQUIR
213219
- All entries in `sources` are instantiated and superimposed during the run. Multiple simultaneous planewave entries are supported.
214220
- **type**: String. Can be `"initial"`, `"planewave"`, or `"dipole"`.
215221

222+
### Volumetric PML workflow
223+
224+
For the first volumetric PML implementation, mesh the physical vacuum region as an inner box and the absorber as a separate outer shell in Gmsh. Tag only the outer shell with a `type: "PML"` material entry in `model.materials`; keep the physical region on ordinary vacuum material tags. The solver infers the shell thickness from the mesh, builds the graded sigma profile internally, and prints startup diagnostics so you can confirm the shell is adequately resolved before running longer simulations.
225+
216226
*For `type: "initial"` — volumetric initial condition:*
217227
- **field_type**: String. Can be `"electric"` or `"magnetic"`.
218228
- **polarization**: Array of 3 doubles. Polarization direction vector.

‎src/components/DGOperatorFactory.h‎

Lines changed: 48 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -299,9 +299,15 @@ namespace maxwell
299299

300300
// Methors for complete Global Operators //
301301

302+
template <typename BF>
303+
std::unique_ptr<BF> buildWeightedMassOperator(mfem::Coefficient& coeff);
304+
302305
template <typename BF>
303306
std::unique_ptr<BF> buildSigmaMassOperator();
304307

308+
template <typename BF>
309+
std::unique_ptr<mfem::SparseMatrix> buildMatchedConductiveOperator(mfem::Coefficient& coeff);
310+
305311
template <typename BF>
306312
void addGlobalZeroNormalIBFIOperators(mfem::SparseMatrix* global);
307313
template <typename BF>
@@ -909,18 +915,55 @@ namespace maxwell
909915

910916
template <typename FES>
911917
template <typename BF>
912-
std::unique_ptr<BF> DGOperatorFactory<FES>::buildSigmaMassOperator()
918+
std::unique_ptr<BF> DGOperatorFactory<FES>::buildWeightedMassOperator(mfem::Coefficient& coeff)
913919
{
914-
Vector sigma = pd_.model.buildSigmaPiecewiseVector();
915-
PWConstCoefficient SigCoeff(sigma);
916-
917920
auto bf = std::make_unique<BF>(&fes_);
918-
bf->AddDomainIntegrator(new MassIntegrator(SigCoeff));
921+
bf->AddDomainIntegrator(new MassIntegrator(coeff));
919922
bf->Assemble();
920923
bf->Finalize();
921924
return bf;
922925
}
923926

927+
template <typename FES>
928+
template <typename BF>
929+
std::unique_ptr<BF> DGOperatorFactory<FES>::buildSigmaMassOperator()
930+
{
931+
Vector sigma = pd_.model.buildSigmaPiecewiseVector();
932+
PWConstCoefficient SigCoeff(sigma);
933+
934+
return buildWeightedMassOperator<BF>(SigCoeff);
935+
}
936+
937+
template <typename FES>
938+
template <typename BF>
939+
std::unique_ptr<mfem::SparseMatrix> DGOperatorFactory<FES>::buildMatchedConductiveOperator(mfem::Coefficient& coeff)
940+
{
941+
auto res = std::make_unique<mfem::SparseMatrix>(6 * fes_.GetNDofs(), 6 * fes_.GetNDofs());
942+
auto MInv = buildMaxwellInverseMassMatrixOperator<BF>();
943+
auto MCond = buildWeightedMassOperator<BF>(coeff);
944+
auto ACondE = buildByMult<FES, BF>(MInv[E]->SpMat(), MCond->SpMat(), fes_);
945+
auto ACondH = buildByMult<FES, BF>(MInv[H]->SpMat(), MCond->SpMat(), fes_);
946+
947+
GlobalIndices gid(fes_.GetNDofs(), 0, true);
948+
for (auto d : { X, Y, Z }) {
949+
loadBlockInGlobalAtIndices(
950+
ACondE->SpMat(),
951+
*res,
952+
std::make_pair(*gid.offsets[E][d].get(), *gid.offsets[E][d].get()),
953+
-1.0
954+
);
955+
loadBlockInGlobalAtIndices(
956+
ACondH->SpMat(),
957+
*res,
958+
std::make_pair(*gid.offsets[H][d].get(), *gid.offsets[H][d].get()),
959+
-1.0
960+
);
961+
}
962+
963+
res->Finalize();
964+
return res;
965+
}
966+
924967
template <typename FES>
925968
template <typename BF>
926969
void DGOperatorFactory<FES>::addGlobalZeroNormalIBFIOperators(SparseMatrix* global)

0 commit comments

Comments
 (0)