Merge pull request #72 from OpenSEMBA/dev

Up-to-date-ing
Tests pass, I bypass

Author: AlejandroMunozManterola

.gitignore

@@ -2,6 +2,7 @@
 .vscode/
 .github/agents/
 .github/instructions/
+.opencode/
 
 build/
 out/

CMakeLists.txt

@@ -17,6 +17,7 @@ option(SEMBA_DGTD_ENABLE_MFEM_AS_SUBDIRECTORY "Use MFEM as a subdirectory" ON )
 option(SEMBA_DGTD_ENABLE_EXTENSIVE_CASE_TESTS "Enable extensive case tests" ON )
 option(SEMBA_DGTD_ENABLE_EXTENSIVE_SOLVER_TESTS "Enable extensive solver tests" ON )
 option(SEMBA_DGTD_ENABLE_EXTENSIVE_RCS_TESTS "Enable extensive RCS tests" ON )
+option(SEMBA_DGTD_ENABLE_L2_ANALYSIS_TESTS "Enable L2 Analysis Tests" OFF)
 
 option(SEMBA_DGTD_ENABLE_TIMER_INFORMATION "Enable timer information" ON)
 

README.md

@@ -102,7 +102,7 @@ MFEM is included as a submodule (`external/mfem-geg`) and built automatically. T
 ### Defining a JSON file
 
 Cases can be defined by parsing a JSON file with the problem information. See a complete [example](https://github.com/OpenSEMBA/dgtd/blob/main/testData/maxwellInputs/1D_PEC/1D_PEC.json) of a valid JSON file.
-An OpenSEMBA/dgtd JSON file must have the following structure; **bold** entries are **required**:
+An OpenSEMBA/dgtd JSON file must have the following structure. Legend: **[REQUIRED]** = must be present; **[OPTIONAL]** = has a default value (shown).
 
 - solver_options:
 	- Object. User can customise solver settings. If undefined, all defaults apply.
@@ -114,18 +114,18 @@ An OpenSEMBA/dgtd JSON file must have the following structure; **bold** entries
 			- Double. Total simulation duration in natural units (1 meter/c). If undefined, defaults to `2.0`.
 		- time_step:
 			- Double. Fixed time step in natural units. Must be defined for 2D and 3D problems. Overrides `cfl` for 1D if both are defined. If undefined or `0.0` in 1D, an automatic step is computed via CFL.
-		- cfl:
-			- Double. Courant–Friedrichs–Lewy condition used to compute the automatic time step in 1D. Not available for 2D or 3D. Ignored if `time_step` is also defined. Must be in the range (0.0, 1.0].
-		- order:
-			- Integer. Polynomial order of the finite element basis. If undefined, defaults to `3`.
+		- **cfl**:
+			- Double. Courant–Friedrichs–Lewy condition used to compute the automatic time step in 1D. Not available for 2D or 3D. Ignored if `time_step` is also defined. Must be in the range (0.0, 1.0]. Defaults to `1.0`.
+		- **order**:
+			- Integer. Polynomial order of the finite element basis. Defaults to `2`.
 		- spectral:
 			- Boolean. Use a spectral evolution operator that assembles the full E/H system matrix and derives the time step from its eigenvalues. High computational cost; does not support all features. If undefined, defaults to `false`.
 		- export_operator:
 			- Boolean. Write the assembled evolution operator matrix to disk for inspection. If undefined, defaults to `false`.
-		- basis_type:
-			- String. Finite element basis type passed to MFEM.
-		- ode_type:
-			- String. ODE time-integration method passed to MFEM.
+		- **basis_type**:
+			- Integer. Finite element basis type passed to MFEM. Defaults to `1` (GaussLobatto). Values: `0` (GaussLegendre), `1` (GaussLobatto), `2` (Bernstein), `3` (OpenUniform), `4` (CloseUniform), `5` (OpenHalfUniform).
+		- **ode_type**:
+			- Integer. ODE time-integration method. Defaults to `0` (RK4). Values: `0` (RK4), `1` (BackwardEuler), `2` (Trapezoidal), `3` (ImplicitMidpoint), `4` (SDIRK33), `5` (SDIRK23), `6` (SDIRK34).
 
 - **model**:
 	- Object. Contains geometry, material, and boundary information.
@@ -168,30 +168,30 @@ An OpenSEMBA/dgtd JSON file must have the following structure; **bold** entries
 	- Object. Controls data extraction. If undefined, no data is recorded.
 		- exporter:
 			- Object. Enables ParaView (VisIt) field export.
-				- name: String. Output dataset name. Defaults to the mesh filename stem.
+				- name: String. Output dataset name. Defaults to mesh filename (e.g. `mesh.msh` → `mesh`).
 				- steps: Integer. Export every N time steps. Mutually exclusive with `saves`.
 				- saves: Integer. Total number of exports over the whole simulation. The step interval is computed automatically. Mutually exclusive with `steps`.
 		- point:
 			- Array. Each entry records all E and H field components at a single point every interval.
 				- **position**: Array of doubles. Spatial coordinates. Must match the mesh dimension (e.g. `[x, y]` for 2D). ***Warning:*** If the point lies outside the mesh the simulation will crash.
-				- steps / saves: See exporter above.
+				- steps / saves: (steps and saves are mutually exclusive; see `exporter` above for details)
 		- field:
 			- Array. Each entry records a single scalar field component at a point.
 				- **field_type**: String. Can be `"electric"` or `"magnetic"`.
 				- **polarization**: String. Component to record. Can be `"X"`, `"Y"`, or `"Z"`.
 				- **position**: Array of doubles. Spatial coordinates. Same constraints as for `point`.
-				- steps / saves: See exporter above.
+				- steps / saves: (steps and saves are mutually exclusive; see `exporter` above for details)
 		- domain_snapshot:
 			- Object. Periodic full-domain field snapshot (alternative to the incremental exporter).
-				- name: String. Output name. Defaults to the mesh filename stem.
-				- steps / saves: See exporter above.
+				- name: String. Output name. Defaults to mesh filename (e.g. `mesh.msh` → `mesh`).
+				- steps / saves: (steps and saves are mutually exclusive; see `exporter` above for details)
 		- rcssurface:
 			- Array. Each entry exports surface E/H field snapshots projected onto a boundary submesh to a compact binary file (`surface_data.bin`). The binary contains a geometry header (quadrature point positions, outward normals, and weights) followed by time-stamped E/H field blocks. Intended for offline RCS post-processing.
 				- **tags**: Array of integers. Mesh surface tags that define the integration surface.
-				- name: String. Output name.
-				- steps / saves: See exporter above.
+				- name: String. [OPTIONAL] Output name. Defaults to `"RCSSurfaceProbe"`.
+				- steps / saves: (steps and saves are mutually exclusive; see `exporter` above for details)
 		- mor_state:
-			- Array. Each entry periodically saves the full DG state vector (all E and H DOFs) to disk over a specified time window. The saved vectors are compatible with the exported global operator matrix (`{name}_global.csr`), so that `y = A * x` can be evaluated offline (e.g. in Python). Each snapshot is written as a plain-text file `x_0`, `x_1`, … containing the vector size on the first line followed by one DOF value per line (16-digit precision). The vector layout is `[Ex₀…ExN | Ey | Ez | Hx | Hy | Hz]`.
+			- Array. Each entry periodically saves the full DG state vector (all E and H DOFs) to disk over a specified time window. The saved vectors are compatible with the exported global operator matrix (`{name}_global.csr`), so that `y = A * x` can be evaluated offline (e.g. in Python). Each snapshot is written as a plain-text file `x_0`, `x_1`, … with: first line = absolute simulation time, second line = vector size, followed by one DOF value per line (16-digit precision). The vector layout is `[Ex₀…ExN | Ey | Ez | Hx | Hy | Hz]`.
 				- **record_time_start**: Double. Simulation time at which recording begins.
 				- **record_time_final**: Double. Simulation time at which recording ends.
 				- **saves**: Integer. Number of snapshots to record, distributed uniformly between `record_time_start` and `record_time_final`.
@@ -208,8 +208,8 @@ An OpenSEMBA/dgtd JSON file must have the following structure; **bold** entries
 			- **type**: String. Can be `"gaussian"`, `"resonant"`, `"besselj6_2D"`, or `"besselj6_3D"`.
 			- *(For `type: "gaussian"`)* **spread**: Double. Standard deviation $\sigma$ of the Gaussian pulse.
 			- *(For `type: "resonant"`)* **modes**: Array of integers. Number of standing waves along each spatial axis.
-		- *(Required for `magnitude.type: "gaussian"`)* **center**: Array of n doubles. Spatial position of the Gaussian centroid.
-		- *(Required for `magnitude.type: "gaussian"`)* **dimension**: Integer. Number of active spatial dimensions in the Gaussian exponent.
+		- **center** (Required for `magnitude.type: "gaussian"`): Array of n doubles. Spatial position of the Gaussian centroid.
+		- **dimension** (Required for `magnitude.type: "gaussian"`): Integer. Number of active spatial dimensions in the Gaussian exponent.
 
 		*For `type: "planewave"` — total-field/scattered-field (TFSF) plane wave:*
 		- **tags**: Array of integers. Mesh boundary tags that form the TFSF interface surface.
@@ -227,6 +227,59 @@ An OpenSEMBA/dgtd JSON file must have the following structure; **bold** entries
 			- **spread**: Double. Gaussian spread parameter.
 			- **mean**: Double. Gaussian center position along the dipole axis.
 
+
+## MOR To ParaView Post-Processing
+
+`opensemba_mor2paraview` replays previously exported `mor_state` snapshots (`x_0`, `x_1`, ...) into a ParaView time series dataset.
+
+Expected folder layout (run command from this folder):
+
+```text
+.
+|-- <case>.json
+|-- <mesh>.msh
+`-- x/
+		|-- x_0
+		|-- x_1
+		`-- ...
+```
+
+Snapshot format per `x_k` file:
+
+1. first line: absolute simulation time
+2. second line: vector size
+3. remaining lines: state values in `[Ex | Ey | Ez | Hx | Hy | Hz]` packed order
+
+Default command:
+
+```sh
+./build/gnu-release-mpi/bin/opensemba_mor2paraview
+```
+
+Default behavior:
+
+- Detects exactly one `.json` in current directory (or errors if none/multiple).
+- Detects exactly one `.msh` in current directory (or errors if none/multiple).
+- Reads snapshots from `./x`.
+- Writes ParaView output to `./output`.
+- Uses dataset name `mor_paraview.vtk` (ParaView collection name).
+
+Optional arguments:
+
+```sh
+./build/gnu-release-mpi/bin/opensemba_mor2paraview \
+	--case ./my_case.json \
+	--mesh ./my_mesh.msh \
+	--xdir ./x \
+	--out ./output \
+	--name mor_paraview.vtk
+```
+
+Notes:
+
+- `opensemba_mor2paraview` supports only single-rank execution.
+- Output is the same ParaView-style time-series format used by the exporter probe (collection plus time-step data), suitable for direct loading in ParaView.
+
 ## Funding
 
 - Spanish Ministry of Science and Innovation (MICIN/AEI) (Grant Number: PID2022-137495OB-C31).

src/components/DGOperatorFactory.h

@@ -1566,7 +1566,9 @@ namespace maxwell
 		// Free sub-operator blocks before applying threshold.
 		blocks.clear();
 
-		auto threshold = 1e-20;
+		// Threshold set to sqrt(eps_machine) ~ 1e-8, the standard criterion for distinguishing
+		// genuine matrix entries from floating-point assembly noise relative to unit-scale quantities.
+		auto threshold = 1e-8;
 		res->Threshold(threshold);
 
 		if(this->pd_.opts.export_evolution_operator){

src/driver/driver.cpp

@@ -1649,6 +1649,12 @@ maxwell::Solver buildSolver(const json& case_data, const std::string& case_path,
 
 	maxwell::SolverOptions solverOpts{ buildSolverOptions(case_data) };
 	maxwell::Probes probes{ buildProbes(case_data) };
+
+	// If MOR state probes are defined, force export_operator to true
+	if (!probes.morStateProbes.empty()) {
+		solverOpts.setExportEO(true);
+	}
+
 	// Model (mesh) must be built before sources so that auto-mean computation
 	// can inspect the TFSF surface geometry to determine the correct delay.
 	maxwell::Model model{ buildModel(case_data, case_path, isTest) };

src/evolution/GlobalEvolution.cpp

@@ -29,8 +29,31 @@ SGBCHelperFields initSGBCHelperFields(const int size)
     return res;
 }
 
+static mfem::Array<int> buildSurfaceMarkerFromTags(const std::vector<int>& tags, const mfem::ParFiniteElementSpace& fes)
+{
+    mfem::Array<int> marker(fes.GetMesh()->bdr_attributes.Max());
+    marker = 0;
+    for (const int t : tags) {
+        marker[t - 1] = 1;
+    }
+    return marker;
+}
+
+static std::vector<int> collectRCSSurfaceTags(const Probes& probes)
+{
+    std::unordered_set<int> unique_tags;
+    for (const auto& p : probes.rcsSurfaceProbes) {
+        for (const int t : p.tags) {
+            unique_tags.insert(t);
+        }
+    }
+    std::vector<int> tags(unique_tags.begin(), unique_tags.end());
+    std::sort(tags.begin(), tags.end());
+    return tags;
+}
+
 GlobalEvolution::GlobalEvolution(
-    mfem::ParFiniteElementSpace& fes, Model& model, SourcesManager& srcmngr, EvolutionOptions& options) :
+    mfem::ParFiniteElementSpace& fes, Model& model, SourcesManager& srcmngr, EvolutionOptions& options, const Probes& probes) :
     mfem::TimeDependentOperator(numberOfFieldComponents* numberOfMaxDimensions* fes.GetNDofs()),
     fes_{ fes },
     model_{ model },
@@ -312,7 +335,7 @@ GlobalEvolution::GlobalEvolution(
         numberOfFieldComponents * numberOfMaxDimensions * (fes_.GetNDofs() + fes_.num_face_nbr_dofs)
     );
 
-    Probes probes; 
+    Probes emptyProbes; 
 
     // Keep TFSF submesh infrastructure for planewave source evaluation
     if (model_.getTotalFieldScatteredFieldToMarker().find(BdrCond::TotalFieldIn) != model_.getTotalFieldScatteredFieldToMarker().end()) {
@@ -323,7 +346,7 @@ GlobalEvolution::GlobalEvolution(
     }
 
     // Build all operators on the global mesh
-    ProblemDescription pd(model_, probes, srcmngr_.sources, opts_);
+    ProblemDescription pd(model_, emptyProbes, srcmngr_.sources, opts_);
     DGOperatorFactory<mfem::ParFiniteElementSpace> dgops(pd, fes_);
     globalOperator_ = dgops.buildGlobalOperator();
 
@@ -344,9 +367,88 @@ GlobalEvolution::GlobalEvolution(
                 TFSFOperator_->PrintCSR2(ofs);
                 ofs.close();
                 std::cout << "TFSF operator exported to " << file_path << std::endl;
+
+                // Export TFSF mapping matrix T (maps TFSF submesh DOFs to parent mesh DOFs)
+                const int tfsf_ndofs = tfsf_sub_to_parent_ids_.Size();
+                const int parent_ndofs = fes_.GetNDofs();
+                auto mapping_matrix = std::make_unique<mfem::SparseMatrix>(6 * parent_ndofs, 6 * tfsf_ndofs);
+
+                for (int comp = 0; comp < 6; ++comp) {
+                    for (int tfsf_dof = 0; tfsf_dof < tfsf_ndofs; ++tfsf_dof) {
+                        int parent_dof = tfsf_sub_to_parent_ids_[tfsf_dof];
+                        int row = comp * parent_ndofs + parent_dof;
+                        int col = comp * tfsf_ndofs + tfsf_dof;
+                        mapping_matrix->Set(row, col, 1.0);
+                    }
+                }
+                mapping_matrix->Finalize();
+
+                std::filesystem::path mapping_path = export_dir / (model_.meshName_ + "_tfsf_mapping.csr");
+                std::ofstream ofs_map(mapping_path);
+                if (!ofs_map.is_open()) {
+                    throw std::runtime_error("Could not open file for writing: " + mapping_path.string());
+                }
+                mapping_matrix->PrintCSR2(ofs_map);
+                ofs_map.close();
+                std::cout << "TFSF mapping matrix exported to " << mapping_path << std::endl;
+            }
+        }
+    }
+
+    const bool has_rcs_probe = !probes.rcsSurfaceProbes.empty();
+    const bool has_mor_probe = !probes.morStateProbes.empty();
+    if (opts_.export_evolution_operator && has_rcs_probe && has_mor_probe) {
+        if (Mpi::WorldSize() == 1) {
+            const auto rcs_tags = collectRCSSurfaceTags(probes);
+            if (!rcs_tags.empty()) {
+                auto marker = buildSurfaceMarkerFromTags(rcs_tags, fes_);
+                NearToFarFieldSubMesher ntff_submesher(model_.getConstMesh(), fes_, marker);
+
+                auto* ntff_submesh = ntff_submesher.getSubMesh();
+                auto dgfec = dynamic_cast<const mfem::DG_FECollection*>(fes_.FEColl());
+                if (!dgfec) {
+                    throw std::runtime_error("The FiniteElementCollection in the FiniteElementSpace is not DG.");
+                }
+
+                mfem::FiniteElementSpace ntff_fes(ntff_submesh, dgfec);
+                mfem::Array<int> farfield_sub_to_parent_ids;
+                mfem::SubMeshUtils::BuildVdofToVdofMap(ntff_fes,
+                                                       fes_,
+                                                       ntff_submesh->GetFrom(),
+                                                       ntff_submesh->GetParentElementIDMap(),
+                                                       farfield_sub_to_parent_ids);
+
+                const int farfield_ndofs = farfield_sub_to_parent_ids.Size();
+                const int parent_ndofs = fes_.GetNDofs();
+                auto farfield_mapping_matrix = std::make_unique<mfem::SparseMatrix>(6 * farfield_ndofs, 6 * parent_ndofs);
+
+                for (int comp = 0; comp < 6; ++comp) {
+                    for (int sub_dof = 0; sub_dof < farfield_ndofs; ++sub_dof) {
+                        const int parent_dof = farfield_sub_to_parent_ids[sub_dof];
+                        const int row = comp * farfield_ndofs + sub_dof;
+                        const int col = comp * parent_ndofs + parent_dof;
+                        farfield_mapping_matrix->Set(row, col, 1.0);
+                    }
+                }
+                farfield_mapping_matrix->Finalize();
+
+                std::filesystem::path export_dir = std::filesystem::path("Exports") / "Operators" / model_.meshName_;
+                if (!std::filesystem::exists(export_dir)) {
+                    std::filesystem::create_directories(export_dir);
+                }
+
+                std::filesystem::path farfield_path = export_dir / (model_.meshName_ + "_farfield.csr");
+                std::ofstream ofs_farfield(farfield_path);
+                if (!ofs_farfield.is_open()) {
+                    throw std::runtime_error("Could not open file for writing: " + farfield_path.string());
+                }
+                farfield_mapping_matrix->PrintCSR2(ofs_farfield);
+                ofs_farfield.close();
+                std::cout << "Farfield mapping matrix exported to " << farfield_path << std::endl;
             }
         }
     }
+
     if (model_.getSGBCToMarker().find(BdrCond::SGBC) != model_.getSGBCToMarker().end()) {
         SGBCOperator_ = dgops.buildSGBCGlobalOperator();
     }

src/evolution/GlobalEvolution.h

@@ -3,6 +3,7 @@
 #include "solver/SourcesManager.h"
 #include "solver/SolverExtension.h"
 #include "components/DGOperatorFactory.h"
+#include "components/Probes.h"
 #include <unordered_map>
 
 namespace maxwell {
@@ -16,7 +17,7 @@ class GlobalEvolution : public mfem::TimeDependentOperator {
     static const int numberOfFieldComponents = 2;
     static const int numberOfMaxDimensions = 3;
 
-    GlobalEvolution(mfem::ParFiniteElementSpace&, Model&, SourcesManager&, EvolutionOptions&);
+    GlobalEvolution(mfem::ParFiniteElementSpace&, Model&, SourcesManager&, EvolutionOptions&, const Probes&);
     ~GlobalEvolution();
 
     virtual void Mult(const mfem::Vector& x, mfem::Vector& y) const;
@@ -28,6 +29,8 @@ class GlobalEvolution : public mfem::TimeDependentOperator {
 
     const mfem::SparseMatrix& getConstGlobalOperator() { return *globalOperator_.get(); }
 
+    const mfem::Array<int>& getTFSFMapping() const { return tfsf_sub_to_parent_ids_; }
+
     bool hasSGBC() const { return !sgbc_states_.empty(); }
 
 private:

src/launcher/CMakeLists.txt

@@ -30,4 +30,18 @@ target_link_libraries(opensemba_rcs
 if (SEMBA_DGTD_ENABLE_OPENMP)
 	message(STATUS "Linking opensemba_rcs with OpenMP libraries")
 	target_link_libraries(opensemba_rcs OpenMP::OpenMP_C OpenMP::OpenMP_CXX)
+endif()
+
+message(STATUS "Creating build system for opensemba_mor2paraview")
+
+add_executable(opensemba_mor2paraview "mor2paraview.cpp")
+
+target_link_libraries(opensemba_mor2paraview
+		maxwell-driver
+		mfem
+)
+
+if (SEMBA_DGTD_ENABLE_OPENMP)
+	message(STATUS "Linking opensemba_mor2paraview with OpenMP libraries")
+	target_link_libraries(opensemba_mor2paraview OpenMP::OpenMP_C OpenMP::OpenMP_CXX)
 endif()