Implement 3D EPR simulations in Elmer

Implement a new solution type `ElmerEPR3DSolution` which:

- Keeps the original (dielectric) layers in Elmer instead of grouping them by material.
   - This enables the use of partition regions in Elmer 3D simulations
- Support for use of tls sheet interfaces in Elmer (`tls_sheet_approximation`)
   - Result energies for tls sheets are saved in normal and tangential components with prefixes `Ez` and `Exy`
   - Correction based on thickness and permittivity can be done in post-processing with `produce_epr_table`
- Solution type does not produce capacitance matrix but only energies
- Add a simulation script `tls_waveguide_sim_elmer` showcasing the new solution type

To enable new solution type, some other changes were also needed

- Add a new keyword `detach_tls_sheets_from_body` in  `Simulation` to control whether tls_layers are shifted from the metal by `tls_layer_thickness`
- Implement custom Elmer Solver module for saving the energies on the 2d tls sheets
    - As a temporary solution the module is automatically compiled at runtime when needed.
    - Compiling at runtime only supported in Linux or WSL
- Gmsh mesh is generation refactored to enable using partition regions
- When using the new solution type Gmsh writes signals and grounds as 3D bodies instead of 2D boundaries and does not write ports

Bugfixes and minor unrelated changes:
- bugfix: warnings from simulation export (stderr stream) is propagated correctly when usig `kqc sim`
- Energies written in project results also for `ElmerCapacitanceSolution`  when using `integrate_energies=True`

klayout_package/python/kqcircuits/simulations/export/elmer/elmer_export.py

@@ -22,6 +22,7 @@
 import logging
 import json
 import argparse
+import platform
 
 from pathlib import Path
 from typing import Sequence, Union, Tuple, Dict, Optional
@@ -34,10 +35,10 @@
 )
 from kqcircuits.simulations.export.util import export_layers
 from kqcircuits.util.export_helper import write_commit_reference_file
-from kqcircuits.defaults import ELMER_SCRIPT_PATHS, KQC_REMOTE_ACCOUNT
+from kqcircuits.defaults import ELMER_SCRIPT_PATHS, KQC_REMOTE_ACCOUNT, SIM_SCRIPT_PATH
 from kqcircuits.simulations.simulation import Simulation
 from kqcircuits.simulations.cross_section_simulation import CrossSectionSimulation
-from kqcircuits.simulations.export.elmer.elmer_solution import ElmerSolution, get_elmer_solution
+from kqcircuits.simulations.export.elmer.elmer_solution import ElmerEPR3DSolution, ElmerSolution, get_elmer_solution
 from kqcircuits.simulations.post_process import PostProcess
 
 
@@ -118,9 +119,11 @@ def export_elmer_script(
     json_filenames,
     path: Path,
     workflow=None,
-    file_prefix="simulation",
-    execution_script="scripts/run.py",
+    script_folder: str = "scripts",
+    file_prefix: str = "simulation",
+    script_file: str = "run.py",
     post_process=None,
+    compile_elmer_modules=False,
 ):
     """
     Create script files for running one or more simulations.
@@ -130,9 +133,11 @@ def export_elmer_script(
         json_filenames: List of paths to json files to be included into the script.
         path: Location where to write the script file.
         workflow: Parameters for simulation workflow
-        file_prefix: Name of the script file to be created.
-        execution_script: The script file to be executed.
+        script_folder: Path to the Elmer-scripts folder.
+        file_prefix: File prefix of the script file to be created.
+        script_file: Name of the script file to run.
         post_process: List of PostProcess objects, a single PostProcess object, or None to be executed after simulations
+        compile_elmer_modules: Compile custom Elmer energy integration module at runtime. Not supported on Windows.
 
     Returns:
 
@@ -145,6 +150,7 @@ def export_elmer_script(
 
     python_executable = workflow.get("python_executable", "python")
     main_script_filename = str(path.joinpath(file_prefix + ".sh"))
+    execution_script = Path(script_folder).joinpath(script_file)
 
     path.joinpath("log_files").mkdir(parents=True, exist_ok=True)
 
@@ -309,6 +315,12 @@ def _divup(a, b):
             main_file.write("export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n")
             main_file.write("\n")
 
+            if compile_elmer_modules:
+                main_file.write('echo "Compiling Elmer modules"\n')
+                main_file.write(
+                    f"elmerf90 -fcheck=all {script_folder}/SaveBoundaryEnergy.F90 -o SaveBoundaryEnergy > /dev/null\n"
+                )
+
             for i, json_filename in enumerate(json_filenames):
                 with open(json_filename) as f:
                     json_data = json.load(f)
@@ -453,9 +465,17 @@ def _divup(a, b):
         with open(main_script_filename, "w") as main_file:
             main_file.write("#!/bin/bash\n")
 
+            if compile_elmer_modules:
+                main_file.write('echo "Compiling Elmer modules"\n')
+                main_file.write(
+                    f"elmerf90 -fcheck=all {script_folder}/SaveBoundaryEnergy.F90 -o SaveBoundaryEnergy > /dev/null\n"
+                )
+
             if parallelize_workload:
                 main_file.write("export OMP_NUM_THREADS={}\n".format(workflow["elmer_n_threads"]))
-                main_file.write("{} scripts/simple_workload_manager.py {}".format(python_executable, n_workers))
+                main_file.write(
+                    "{} {}/simple_workload_manager.py {}".format(python_executable, script_folder, n_workers)
+                )
 
             for i, json_filename in enumerate(json_filenames):
                 with open(json_filename) as f:
@@ -574,8 +594,13 @@ def export_elmer(
 
     workflow = _update_elmer_workflow(simulations, common_sol, workflow)
 
+    # If doing 3D epr simulations the custom Elmer energy integration module is compiled at runtime
+    epr_sim = _is_epr_sim(simulations, common_sol)
+    script_paths = ELMER_SCRIPT_PATHS + [SIM_SCRIPT_PATH / "elmer_modules"] if epr_sim else ELMER_SCRIPT_PATHS
+
     write_commit_reference_file(path)
-    copy_content_into_directory(ELMER_SCRIPT_PATHS, path, script_folder)
+    copy_content_into_directory(script_paths, path, script_folder)
+
     json_filenames = []
 
     for sim_sol in simulations:
@@ -601,12 +626,28 @@ def export_elmer(
         json_filenames,
         path,
         workflow,
+        script_folder=script_folder,
         file_prefix=file_prefix,
-        execution_script=Path(script_folder).joinpath(script_file),
+        script_file=script_file,
         post_process=post_process,
+        compile_elmer_modules=epr_sim,
     )
 
 
+def _is_epr_sim(simulations, common_sol):
+    """Helper to check if doing 3D epr simulation"""
+    epr_sim = False
+    if common_sol is None:
+        if any(isinstance(simsol[1], ElmerEPR3DSolution) for simsol in simulations):
+            epr_sim = True
+    elif isinstance(common_sol, ElmerEPR3DSolution):
+        epr_sim = True
+
+    if epr_sim and platform.system() == "Windows":
+        logging.warning("Elmer 3D EPR Simulations are not supported on Windows")
+    return epr_sim
+
+
 def _update_elmer_workflow(simulations, common_solution, workflow):
     """
     Modify workflow based on number of simulations and available computing resources

klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py

@@ -173,14 +173,36 @@ class ElmerCrossSectionSolution(ElmerSolution):
     run_inductance_sim: bool = True
 
 
+@dataclass
+class ElmerEPR3DSolution(ElmerSolution):
+    """
+    Class for Elmer 3D EPR simulations. Similar to electrostatics simulations done with ElmerCapacitanceSolution,
+    but supports separating energies by PartitionRegions and produces no capacitance matrix
+
+    Args:
+        p_element_order: polynomial order of p-elements
+        linear_system_method: Options: 1. mg (multigrid), 2. bicgstab or any other iterative solver mentioned in
+                              ElmerSolver manual section 4.3.1
+        convergence_tolerance: Convergence tolerance of the iterative solver.
+        max_iterations: Maximum number of iterations for the iterative solver.
+    """
+
+    tool: ClassVar[str] = "epr_3d"
+
+    p_element_order: int = 3
+    linear_system_method: str = "bicgstab"
+    convergence_tolerance: float = 1.0e-9
+    max_iterations: int = 1000
+
+
 def get_elmer_solution(tool="capacitance", **solution_params):
     """Returns an instance of ElmerSolution subclass.
 
     Args:
         tool: Determines the subclass of ElmerSolution.
         solution_params: Arguments passed for  ElmerSolution subclass.
     """
-    for c in [ElmerVectorHelmholtzSolution, ElmerCapacitanceSolution, ElmerCrossSectionSolution]:
+    for c in [ElmerVectorHelmholtzSolution, ElmerCapacitanceSolution, ElmerCrossSectionSolution, ElmerEPR3DSolution]:
         if tool == c.tool:
             return c(**solution_params)
     raise ValueError(f"No ElmerSolution found for tool={tool}.")

klayout_package/python/kqcircuits/simulations/export/export_and_run.py

@@ -74,10 +74,10 @@ def run_export_script(export_script: Path, export_path: Path, quiet: bool = Fals
     )
     # Run export script and capture stdout to be processed
     with subprocess.Popen(export_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) as process:
-        process_stdout, _ = process.communicate()
+        process_stdout, process_stderr = process.communicate()
         print(process_stdout)
+        print(process_stderr, file=sys.stderr)
         if process.returncode:
-            # This provides full traceback so error doesn't need to be captured in process.communicate()
             raise subprocess.CalledProcessError(process.returncode, export_cmd)
 
     # Parse export paths from stdout printed in `create_or_empty_tmp_directory`

klayout_package/python/kqcircuits/simulations/simulation.py

@@ -240,6 +240,13 @@ class Simulation:
         docstring="Use only keywords introduced in material_dict.",
     )
     tls_sheet_approximation = Param(pdt.TypeBoolean, "Approximate TLS interface layers as sheets", False)
+    detach_tls_sheets_from_body = Param(
+        pdt.TypeBoolean,
+        "TLS interface layers are created `tls_layer_thickness` above or below the interface of 3D bodies",
+        True,
+        docstring="Only has an effect when tls_sheet_approximation=True."
+        "Setting to False when using `ElmerEPR3DSolution` significantly improves simulation performance",
+    )
 
     minimum_point_spacing = Param(pdt.TypeDouble, "Tolerance for merging adjacent points in polygon", 0.01, unit="µm")
     polygon_tolerance = Param(pdt.TypeDouble, "Tolerance for merging adjacent polygons in a layer", 0.004, unit="µm")
@@ -699,7 +706,10 @@ def create_simulation_layers(self):
                 layer_top_z = layer_z + [sign, -sign, -sign][layer_num] * thickness
                 material = self.ith_value(self.tls_layer_material, layer_num)
                 if self.tls_sheet_approximation:
-                    z_params = {"z0": layer_top_z, "z1": layer_top_z}
+                    if self.detach_tls_sheets_from_body:
+                        z_params = {"z0": layer_top_z, "z1": layer_top_z}
+                    else:
+                        z_params = {"z0": layer_z, "z1": layer_z}
                 elif thickness != 0.0:
                     z_params = {"z0": layer_z, "z1": layer_top_z}
 

klayout_package/python/scripts/simulations/elmer/cross_section_helpers.py

@@ -16,7 +16,6 @@
 # Please see our contribution agreements for individuals (meetiqm.com/iqm-individual-contributor-license-agreement)
 # and organizations (meetiqm.com/iqm-organization-contributor-license-agreement).
 
-import re
 from pathlib import Path
 
 import numpy as np
@@ -317,30 +316,3 @@ def get_cross_section_capacitance_and_inductance(json_data, folder_path):
         return {"Cs": c_matrix.tolist(), "Ls": None}
 
     return {"Cs": c_matrix.tolist(), "Ls": l_matrix.tolist()}
-
-
-def get_energy_integrals(path):
-    """
-    Return electric energy integrals
-
-    Args:
-        path(Path): folder path of the model result files
-
-    Returns:
-        (dict): energies stored
-    """
-    try:
-        energy_data, energy_layer_data = Path(path) / "energy.dat", Path(path) / "energy.dat.names"
-        energies = pd.read_csv(energy_data, sep=r"\s+", header=None).values.flatten()
-
-        with open(energy_layer_data) as fp:
-            energy_layers = [
-                match.group(1)
-                for line in fp
-                for match in re.finditer("diffusive energy: potential mask ([a-z_0-9]+)", line)
-            ]
-
-        return {f"E_{k}": energy for k, energy in zip(energy_layers, energies)}
-
-    except FileNotFoundError:
-        return {"total_energy": None}