diff --git a/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_export.py b/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_export.py index f6d759eea..21a147af3 100644 --- a/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_export.py +++ b/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 diff --git a/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py b/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py index 58446d87d..2aa6de148 100644 --- a/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py +++ b/klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py @@ -173,6 +173,28 @@ 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. @@ -180,7 +202,7 @@ def get_elmer_solution(tool="capacitance", **solution_params): 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}.") diff --git a/klayout_package/python/kqcircuits/simulations/export/export_and_run.py b/klayout_package/python/kqcircuits/simulations/export/export_and_run.py index e72c015b2..1172452b5 100644 --- a/klayout_package/python/kqcircuits/simulations/export/export_and_run.py +++ b/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` diff --git a/klayout_package/python/kqcircuits/simulations/simulation.py b/klayout_package/python/kqcircuits/simulations/simulation.py index 7b040eb0c..e3bedea13 100644 --- a/klayout_package/python/kqcircuits/simulations/simulation.py +++ b/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} diff --git a/klayout_package/python/scripts/simulations/elmer/cross_section_helpers.py b/klayout_package/python/scripts/simulations/elmer/cross_section_helpers.py index 78a09586d..321f76e47 100644 --- a/klayout_package/python/scripts/simulations/elmer/cross_section_helpers.py +++ b/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} diff --git a/klayout_package/python/scripts/simulations/elmer/elmer_helpers.py b/klayout_package/python/scripts/simulations/elmer/elmer_helpers.py index 8de45371a..d60d3560d 100644 --- a/klayout_package/python/scripts/simulations/elmer/elmer_helpers.py +++ b/klayout_package/python/scripts/simulations/elmer/elmer_helpers.py @@ -17,6 +17,7 @@ # and organizations (meetiqm.com/iqm-organization-contributor-license-agreement). # pylint: disable=too-many-lines import csv +import re import json import logging import time @@ -26,10 +27,11 @@ from scipy.constants import epsilon_0 from scipy.signal import find_peaks import numpy as np +import pandas as pd def read_mesh_names(path): - """Returns names from mesh.names file""" + """Returns all names from mesh.names file""" list_of_names = [] with open(path.joinpath("mesh.names")) as file: for line in file: @@ -40,6 +42,41 @@ def read_mesh_names(path): return list_of_names +def read_mesh_bodies(path): + """Returns names of bodies from mesh.names file""" + list_of_names = [] + with open(path.joinpath("mesh.names")) as file: + for line in file: + if "! ----- names for boundaries -----" in line: + break + if line.startswith("$ "): + eq_sign = line.find(" =") + if eq_sign > 2: + list_of_names.append(line[2:eq_sign]) + return list_of_names + + +def read_mesh_boundaries(path): + """Returns names of boundaries from mesh.names file""" + with open(path.joinpath("mesh.names")) as file: + lines = [line.strip() for line in file] + + i = 0 + while i < len(lines) and "! ----- names for boundaries -----" not in lines[i]: + i += 1 + + list_of_names = [] + for line in lines[i + 1 :]: + if line.startswith("$ "): + eq_sign = line.find(" =") + if eq_sign > 2: + list_of_names.append(line[2:eq_sign]) + else: + logging.warning(f"Unexpected mesh boundary name: {line}") + break + return list_of_names + + def coordinate_scaling(json_data: Dict[str, Any]) -> float: """ Returns coordinate scaling, which is determined by parameters 'units' in json_data. @@ -61,6 +98,7 @@ def sif_common_header( def_file=None, dim="3", discontinuous_boundary=False, + constraint_modes_analysis=True, ) -> str: """ Returns common header and simulation blocks of a sif file in string format. @@ -81,7 +119,7 @@ def sif_common_header( res += sif_block( "Run Control", [ - "Constraint Modes Analysis = True", + f"Constraint Modes Analysis = {constraint_modes_analysis}", ] + reset_adaptive_remesh_str, ) @@ -434,6 +472,7 @@ def get_electrostatics_solver( minimum_passes=1, convergence_tolerance=-1, max_iterations=-1, + c_matrix_output=True, ): """ Returns electrostatics solver in sif file format. @@ -448,6 +487,7 @@ def get_electrostatics_solver( max_outlier_fraction(float): Maximum fraction of outliers from the total number of elements minimum_passes(int): Maximum number of adaptive meshing iterations. minimum_passes(int): Minimum number of adaptive meshing iterations. + c_matrix_output(bool): Can be used to turn off capacitance matrix output Returns: (str): electrostatics solver in sif file format @@ -459,7 +499,7 @@ def get_electrostatics_solver( "Equation = Electro Statics", f'Procedure = "{solver}" "StatElecSolver"', "Variable = Potential", - "Calculate Capacitance Matrix = True", + f"Calculate Capacitance Matrix = {c_matrix_output}", "Calculate Electric Field = True", "Calculate Elemental Fields = True", "Average Within Materials = False", @@ -669,7 +709,7 @@ def get_save_data_solver(ordinate, result_file="results.dat"): return sif_block(f"Solver {ordinate}", solver_lines) -def get_save_energy_solver(ordinate, energy_file, bodies): +def get_save_energy_solver(ordinate, energy_file, bodies, sheet_bodies=None): """ Returns save energy solver in sif file format. @@ -677,6 +717,7 @@ def get_save_energy_solver(ordinate, energy_file, bodies): ordinate(int): solver ordinate energy_file(str): data file name for energy results bodies(list(str)): body names for energy calculation + sheet_bodies(list(str)): boundary names for energy calculation in 3D simulation Returns: (str): save energy solver in sif file format @@ -687,21 +728,26 @@ def get_save_energy_solver(ordinate, energy_file, bodies): 'Procedure = "SaveData" "SaveScalars"', f"Filename = {energy_file}", "Parallel Reduce = Logical True", - # Add all target bodies to the solver - *( - line - for layer_props in ( - ( - f"Variable {i} = Potential", - f"Operator {i} = body diffusive energy", - f"Mask Name {i} = {interface}", - f"Coefficient {i} = Relative Permittivity", - ) - for i, interface in enumerate(bodies, 1) - ) - for line in layer_props - ), ] + # Add all target bodies to the solver + for i, interface in enumerate(bodies, 1): + solver_lines += [ + f"Variable {i} = Potential", + f"Operator {i} = body diffusive energy", + f"Mask Name {i} = {interface}", + f"Coefficient {i} = Relative Permittivity", + ] + + # Add sheet body energies using a custom energy solver separating the energy into normal and tangential components + if sheet_bodies is not None: + i = len(bodies) + 1 + for interface in sheet_bodies: + solver_lines += [ + f"Variable {i} = {interface}_norm_component", + f"Variable {i+1} = {interface}_tan_component", + ] + i += 2 + return sif_block(f"Solver {ordinate}", solver_lines) @@ -801,15 +847,17 @@ def produce_sif_files(json_data: dict, path: Path) -> List[Path]: """ path.mkdir(exist_ok=True, parents=True) sif_names = json_data["sif_names"] - - if json_data["tool"] == "capacitance" and len(sif_names) != 1: + tool = json_data["tool"] + if tool == "capacitance" and len(sif_names) != 1: logging.warning(f"Capacitance tool only supports 1 sif name, given {len(sif_names)}") sif_filepaths = [] for ind, sif in enumerate(sif_names): - if json_data["tool"] == "capacitance": + if tool == "capacitance": content = sif_capacitance(json_data, path, vtu_name=path, angular_frequency=0, dim=3, with_zero=False) - elif json_data["tool"] == "wave_equation": + elif tool == "epr_3d": + content = sif_epr_3d(json_data, path, vtu_name=path) + elif tool == "wave_equation": freqs = json_data["frequency"] if len(freqs) != len(sif_names): logging.warning( @@ -817,7 +865,7 @@ def produce_sif_files(json_data: dict, path: Path) -> List[Path]: ) content = sif_wave_equation(json_data, path, frequency=freqs[ind]) else: - logging.warning(f"Unkown tool: {json_data['tool']}. No sif file created") + logging.warning(f"Unkown tool: {tool}. No sif file created") return [] sif_filepath = path.joinpath(f"{sif}.sif") @@ -929,6 +977,154 @@ def get_grounds(json_data, dim, mesh_names): return [] +def sif_placeholder_boundaries(groups, n_boundaries): + boundary_conditions = "" + for i, s in enumerate(groups, 1): + boundary_conditions += sif_boundary_condition( + ordinate=i + n_boundaries, + target_boundaries=[s], + conditions=[ + "! This BC does not do anything, but", + "! MMG does not conserve GeometryIDs if there is no BC defined.", + ], + ) + return boundary_conditions + + +def sif_epr_3d(json_data: dict, folder_path: Path, vtu_name: str): + """ + Returns 3D EPR simulation sif + + This solution assumes that signals and grounds are 3d bodies + Also mesh boundaries are assumed to contain only tls layers and no ports + + Args: + json_data: all the model data produced by `export_elmer_json` + folder_path: folder path of the model files + vtu_name: name of the paraview file + """ + + header = sif_common_header(json_data, folder_path, angular_frequency=0, dim=3, constraint_modes_analysis=False) + constants = sif_block("Constants", [f"Permittivity Of Vacuum = {epsilon_0}"]) + + solvers = get_electrostatics_solver( + ordinate=1, + capacitance_file="none", + method=json_data["linear_system_method"], + p_element_order=json_data["p_element_order"], + maximum_passes=json_data["maximum_passes"], + minimum_passes=json_data["minimum_passes"], + percent_error=json_data["percent_error"], + max_error_scale=json_data["max_error_scale"], + max_outlier_fraction=json_data["max_outlier_fraction"], + convergence_tolerance=json_data["convergence_tolerance"], + max_iterations=json_data["max_iterations"], + c_matrix_output=False, + ) + solvers += get_result_output_solver( + ordinate=2, + output_file_name=vtu_name, + exec_solver="Always" if json_data.get("vtu_output", True) else "Never", + ) + equations = get_equation( + ordinate=1, + solver_ids=[1], + ) + + mesh_bodies = read_mesh_bodies(folder_path) + mesh_boundaries = read_mesh_boundaries(folder_path) + + grounds = [n for n in mesh_bodies if n.startswith("ground")] + signals = [n for n in mesh_bodies if n.startswith("signal")] + mesh_bodies = [n for n in mesh_bodies if n not in grounds + signals] + + permittivity_list = [] + for b in mesh_bodies: + # mesh names have gmsh prefix if layer starts with number + mat = json_data["layers"].get(b.removeprefix("gmsh_"), {}).get("material", "vacuum") + perm = 1.0 + if mat in json_data["material_dict"]: + perm = json_data["material_dict"][mat]["permittivity"] + elif mat != "vacuum": + logging.warning(f"Material {mat} not in material_dict. Using permittivity 1.0") + permittivity_list.append(perm) + + # Solver(s) with masks for saving energy + solver_lines = [ + "Exec Solver = " + ("Always" if len(mesh_boundaries) > 0 else "Never"), + 'Equation = "SaveBoundaryEnergy"', + 'Procedure = "SaveBoundaryEnergy" "SaveBoundaryEnergyComponents"', + ] + solvers += sif_block("Solver 3", solver_lines) + + solvers += get_save_energy_solver( + ordinate=4, + energy_file="energy.dat", + bodies=mesh_bodies, + sheet_bodies=mesh_boundaries, + ) + + bodies = "" + materials = "" + n_bodies = 0 + for i, (body, perm) in enumerate(zip(mesh_bodies, permittivity_list), 1): + bodies += sif_body( + ordinate=i, target_bodies=[body], equation=1, material=i, keywords=[f"{body} = Logical True"] + ) + materials += sif_block(f"Material {i}", [f"Relative Permittivity = {perm}"]) + n_bodies += 1 + + # add material for pec + pec_material_index = n_bodies + 1 + materials += sif_block(f"Material {pec_material_index}", ["Relative Permittivity = 1.0"]) + + if len(signals) == 0: + raise RuntimeError("No signals in the system!") + if len(signals) > 1: + logging.warning("Multiple signals in the system. All of them will be set to 1V at once") + + # grounds + bodies += sif_body( + ordinate=n_bodies + 1, + target_bodies=grounds, + equation=1, + material=pec_material_index, + keywords=["Body Force = 2"], + ) + n_bodies += 1 + + # signals + bodies += sif_body( + ordinate=n_bodies + 1, + target_bodies=signals, + equation=1, + material=pec_material_index, + keywords=["Body Force = 1"], + ) + n_bodies += 1 + + body_forces = "" + body_forces += sif_block("Body Force 1", ["Potential = Real 1.0"]) + body_forces += sif_block("Body Force 2", ["Potential = Real 0.0"]) + + boundary_conditions = "" + # tls bcs + for i, s in enumerate(mesh_boundaries, 1): + boundary_conditions += sif_boundary_condition( + ordinate=i, + target_boundaries=[s], + conditions=[f'Boundary Energy Name = String "{s}"'], + ) + + # If there are no boundary conditions, add one to suppress warnings + if not boundary_conditions: + boundary_conditions += sif_block( + "Boundary Condition 1", ["Target Boundaries(0)", "! Placeholder boundary to suppress warnings"] + ) + + return header + constants + solvers + equations + materials + bodies + body_forces + boundary_conditions + + def sif_capacitance( json_data: dict, folder_path: Path, vtu_name: str, angular_frequency: float, dim: int, with_zero: bool = False ): @@ -1044,15 +1240,7 @@ def sif_capacitance( for n in mesh_names if n not in body_list + ground_boundaries + signals_boundaries + outer_bc_names and not n.startswith("port_") ] - for i, s in enumerate(other_groups, 1): - boundary_conditions += sif_boundary_condition( - ordinate=i + n_boundaries, - target_boundaries=[s], - conditions=[ - "! This BC does not do anything, but", - "! MMG does not conserve GeometryIDs if there is no BC defined.", - ], - ) + boundary_conditions += sif_placeholder_boundaries(other_groups, n_boundaries) n_boundaries += len(other_groups) return header + constants + solvers + equations + materials + bodies + boundary_conditions @@ -1154,16 +1342,7 @@ def sif_inductance(json_data, folder_path, angular_frequency, circuit_definition # Add place-holder boundaries (if additional physical groups are given) other_groups = [n for n in mesh_names if n not in body_list and not n.startswith("port_")] - boundary_conditions = "" - for i, s in enumerate(other_groups, 1): - boundary_conditions += sif_boundary_condition( - ordinate=i, - target_boundaries=[s], - conditions=[ - "! This BC does not do anything, but", - "! MMG does not conserve GeometryIDs if there is no BC defined.", - ], - ) + boundary_conditions = "" + sif_placeholder_boundaries(other_groups, 0) return header + equations + solvers + materials + bodies + components + body_force + boundary_conditions @@ -1553,6 +1732,46 @@ def read_result_smatrix(s_matrix_filename: str, path: Path = None, polar_form: b return smatrix_full +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() + + energy_layers = [] + with open(energy_layer_data) as fp: + reached_data = False + for line in fp: + if not reached_data and "Variables in columns of matrix:" in line: + reached_data = True + continue + + if reached_data: + matches = [ + match.group(1) for match in re.finditer("diffusive energy: potential mask ([a-z_0-9]+)", line) + ] + if matches: + energy_layers += matches + else: + matches = line.strip().partition(": ")[2] + if matches: + energy_layers += [matches] + + return {f"E_{k}": energy for k, energy in zip(energy_layers, energies)} + + except FileNotFoundError: + logging.warning(f"Energy file not found in {path}") + return {"total_energy": None} + + def write_project_results_json(json_data: dict, path: Path, msh_filepath, polar_form: bool = True): """ Writes the solution data in '_project_results.json' format for one Elmer simulation. @@ -1586,17 +1805,50 @@ def write_project_results_json(json_data: dict, path: Path, msh_filepath, polar_ for net_j in range(len(c_matrix)) for net_i in range(len(c_matrix)) } + results = { + "CMatrix": c_matrix, + "Cdata": c_data, + "Frequency": [0], + } + if json_data["integrate_energies"]: + results.update(get_energy_integrals(sif_folder)) with open(json_filename, "w") as outfile: json.dump( - { - "CMatrix": c_matrix, - "Cdata": c_data, - "Frequency": [0], - }, + results, outfile, indent=4, ) + elif tool == "epr_3d": + + def _rename_energy_key(e_name): + # Remove "gmsh_" prefix from layer names + e_name = e_name.replace("_gmsh_", "_") + # Change 2-component energies to same naming as used in Ansys + if e_name.endswith("_norm_component"): + e_name = "Ez" + e_name.removesuffix("_norm_component").removeprefix("E") + elif e_name.endswith("_tan_component"): + e_name = "Exy" + e_name.removesuffix("_tan_component").removeprefix("E") + + # Elmer forces all keys to be lowercase so let's capitalise the interface abbreviations + # to correspond to Ansys naming + for k in ("MA", "MS", "SA"): + elmer_int_key = "_layer" + k.lower() + if elmer_int_key in e_name: + e_name = e_name.replace(elmer_int_key, "_layer" + k) + break + + return e_name + + energies = {_rename_energy_key(k): v for k, v in get_energy_integrals(sif_folder).items()} + + with open(json_filename, "w") as outfile: + json.dump( + energies, + outfile, + indent=4, + ) + elif tool == "wave_equation": frequencies = json_data["frequency"] diff --git a/klayout_package/python/scripts/simulations/elmer/gmsh_helpers.py b/klayout_package/python/scripts/simulations/elmer/gmsh_helpers.py index 03a974cc6..e4197c4e4 100644 --- a/klayout_package/python/scripts/simulations/elmer/gmsh_helpers.py +++ b/klayout_package/python/scripts/simulations/elmer/gmsh_helpers.py @@ -156,20 +156,36 @@ def produce_mesh(json_data, msh_file): gmsh_n_threads = int(n_threads_dict.get("gmsh_n_threads", 1)) set_meshing_options(mesh_field_ids, mesh_global_max_size, gmsh_n_threads) - # Group dim tags by material + tls_sheet_layers = { + name: data + for name, data in layers.items() + if any(("_layer" + layer_key in name for layer_key in ["MA", "MS", "SA"])) and data["thickness"] == 0.0 + } + accepted_thin_materials = ["pec"] + accepted_thin_materials += list({data["material"] for _, data in tls_sheet_layers.items()}) + materials = list(set(accepted_thin_materials + list(json_data["material_dict"].keys()) + ["vacuum"])) - material_dim_tags = {m: [] for m in materials} + + filtered_tags_with_material = {} # store material as the 3rd component of the tuple for name, data in layers.items(): material = data.get("material", None) if material in accepted_thin_materials: - material_dim_tags[material] += [(d, t) for d, t in new_tags.get(name, []) if d in [2, 3]] + filtered_tags_with_material[name] = [(d, t, material) for d, t in new_tags.get(name, []) if d in [2, 3]] elif material in materials: - material_dim_tags[material] += [(d, t) for d, t in new_tags.get(name, []) if d == 3] + filtered_tags_with_material[name] = [(d, t, material) for d, t in new_tags.get(name, []) if d == 3] edge_material = data.get("edge_material", None) if edge_material in accepted_thin_materials: - material_dim_tags[edge_material] += [(d, t) for d, t in new_tags.get("&" + name, []) if d == 2] + filtered_tags_with_material["gmsh_" + "&" + name] += [ + (d, t, material) for d, t in new_tags.get("&" + name, []) if d == 2 + ] + + # Group dim tags by material + material_dim_tags = {m: [] for m in materials} + for _, tags in filtered_tags_with_material.items(): + material = tags[0][2] + material_dim_tags[material] += [(d, t) for d, t, _ in tags] # Sort boundaries of material_dim_tags['pec'] into pec_islands and leave the bodies into material_dim_tags['pec'] pec_with_boundaries = get_recursive_children(material_dim_tags["pec"]).union(material_dim_tags["pec"]) @@ -214,37 +230,67 @@ def _find_island_at(_location, _islands): material_dim_tags[f"ground_{counter}"] = pec_island counter += 1 - # set domain boundary as ground for wave equation simulations ports_dts = [dt for port in ports for dt in new_tags.get(f'port_{port["number"]}', [])] - if json_data.get("tool", "capacitance") == "wave_equation": + # set domain boundary as ground for wave equation simulations + if json_data["tool"] == "wave_equation": solid_dts = [(d, t) for dts in material_dim_tags.values() for d, t in dts if d == 3] face_dts = [(d, t) for dt in solid_dts for d, t in get_recursive_children([dt]) if d == 2] material_dim_tags[f"ground_{counter}"] = [d for d in face_dts if face_dts.count(d) == 1 and d not in ports_dts] - # Add physical groups - for mat, dts in material_dim_tags.items(): + if json_data["tool"] == "epr_3d": + # replace 2d pec boundaries with the corresponding 3d bodies + # also remove 3d pec body containing all of these + + phys_group_tags = { + # Elmer doesn't support layer names starting with number so prefix them with gmsh_ + (k if k[0].isalpha() else "gmsh_" + k): [(d, t) for d, t, _ in dts] + for k, dts in filtered_tags_with_material.items() + if dts[0][2] != "pec" + } + + if material_dim_tags["pec"]: + for name, dt_list in material_dim_tags.items(): + if name.startswith("signal_") or name.startswith("ground_"): + signal_tags_3d = set() + for dim, tag in dt_list: + parent_pec_tags_3d, _ = gmsh.model.getAdjacencies(dim, tag) + parent_pec_tags_3d = [ + (3, tag) for tag in parent_pec_tags_3d if (3, tag) in material_dim_tags["pec"] + ] + signal_tags_3d.update(parent_pec_tags_3d) + phys_group_tags[name] = list(signal_tags_3d) + + else: + phys_group_tags = material_dim_tags + + # add physical groups + for name, dts in phys_group_tags.items(): no_port_dts = [dt for dt in dts if dt not in ports_dts] if no_port_dts: - gmsh.model.addPhysicalGroup(max(dt[0] for dt in no_port_dts), [dt[1] for dt in no_port_dts], name=f"{mat}") - - for port in ports: - port_name = f'port_{port["number"]}' - if port_name in new_tags: - if port["type"] == "EdgePort": - port_dts = set(new_tags[port_name]) - key_dts = {"signal": [], "ground": []} - for mat, dts in material_dim_tags.items(): - if mat.startswith("signal"): - key_dts["signal"] += [(d, t) for d, t in port_dts.intersection(dts) if d == 2] - elif mat.startswith("ground"): - key_dts["ground"] += [(d, t) for d, t in port_dts.intersection(dts) if d == 2] - elif mat != "pec": - key_dts[mat] = [(d, t) for d, t in port_dts.intersection(get_recursive_children(dts)) if d == 2] - for key, dts in key_dts.items(): - if dts: - gmsh.model.addPhysicalGroup(2, [dt[1] for dt in dts], name=f"{port_name}_{key}") - else: - gmsh.model.addPhysicalGroup(2, [dt[1] for dt in new_tags[port_name]], name=port_name) + gmsh.model.addPhysicalGroup(max(dt[0] for dt in no_port_dts), [dt[1] for dt in no_port_dts], name=f"{name}") + + if json_data["tool"] != "epr_3d": + # port physical groups + for port in ports: + port_name = f'port_{port["number"]}' + if port_name in new_tags: + if port["type"] == "EdgePort": + port_dts = set(new_tags[port_name]) + key_dts = {"signal": [], "ground": []} + for mat, dts in material_dim_tags.items(): + if mat.startswith("signal"): + key_dts["signal"] += [(d, t) for d, t in port_dts.intersection(dts) if d == 2] + elif mat.startswith("ground"): + key_dts["ground"] += [(d, t) for d, t in port_dts.intersection(dts) if d == 2] + elif mat != "pec": + key_dts[mat] = [ + (d, t) for d, t in port_dts.intersection(get_recursive_children(dts)) if d == 2 + ] + for key, dts in key_dts.items(): + if dts: + gmsh.model.addPhysicalGroup(2, [dt[1] for dt in dts], name=f"{port_name}_{key}") + else: + gmsh.model.addPhysicalGroup(2, [dt[1] for dt in new_tags[port_name]], name=port_name) # Generate and save mesh gmsh.model.mesh.generate(3) diff --git a/klayout_package/python/scripts/simulations/elmer/run.py b/klayout_package/python/scripts/simulations/elmer/run.py index 601e1cace..2287eea9f 100644 --- a/klayout_package/python/scripts/simulations/elmer/run.py +++ b/klayout_package/python/scripts/simulations/elmer/run.py @@ -21,13 +21,12 @@ from interpolating_frequency_sweep import interpolating_frequency_sweep from gmsh_helpers import produce_mesh -from elmer_helpers import produce_sif_files, write_project_results_json +from elmer_helpers import produce_sif_files, write_project_results_json, get_energy_integrals from run_helpers import run_elmer_grid, run_elmer_solver, run_paraview, write_simulation_machine_versions_file from cross_section_helpers import ( produce_cross_section_mesh, produce_cross_section_sif_files, get_cross_section_capacitance_and_inductance, - get_energy_integrals, ) parser = argparse.ArgumentParser(description="Run script for Gmsh-Elmer workflow") diff --git a/klayout_package/python/scripts/simulations/elmer_modules/SaveBoundaryEnergy.F90 b/klayout_package/python/scripts/simulations/elmer_modules/SaveBoundaryEnergy.F90 new file mode 100644 index 000000000..7264126b9 --- /dev/null +++ b/klayout_package/python/scripts/simulations/elmer_modules/SaveBoundaryEnergy.F90 @@ -0,0 +1,252 @@ + !------------------------------------------------------------------------------ + ! Compile with "elmerf90 SaveBoundaryEnergy.F90 -o obj/SaveBoundaryEnergy" + ! To check for warnings/debugging also add -fcheck=all -Wall -Wextra -g + !------------------------------------------------------------------------------ +SUBROUTINE SaveBoundaryEnergyComponents(Model,Solver,dt,Transient) + !------------------------------------------------------------------------------ + USE DefUtils + + IMPLICIT NONE + !------------------------------------------------------------------------------ + TYPE(Solver_t) :: Solver + TYPE(Model_t) :: Model + REAL(KIND=dp) :: dt + LOGICAL :: Transient + !------------------------------------------------------------------------------ + TYPE(Variable_t), POINTER :: evar + TYPE(Element_t), POINTER :: Element + INTEGER ::i, Active, iMode, NoModes + CHARACTER(LEN=MAX_NAME_LEN) :: bc_name !, debug_str + LOGICAL :: Found + TYPE(Mesh_t), POINTER :: Mesh + REAL(KIND=dp), ALLOCATABLE :: IntNorm(:), IntTan(:) + TYPE(ValueList_t), POINTER :: BC + CHARACTER(*), PARAMETER :: Caller = 'SaveBoundaryEnergyComponents' + + SAVE IntNorm, IntTan, NoModes + !------------------------------------------------------------------------------------------- + + CALL Info(Caller,'----------------------------------------------------------',Level=6 ) + CALL Info(Caller,'Computing Boundary energy integrals for normal and tangential components of Electric field',Level=4 ) + + Mesh => GetMesh() + + NoModes = Model % NumberOfBCs + ALLOCATE(IntNorm(NoModes), IntTan(NoModes)) + + CALL Info(Caller,'Assuming electric field living on nodes') + + evar => VariableGet( Mesh % Variables, 'Electric field e', ThisOnly = .TRUE.) + IF(.NOT. ASSOCIATED( evar ) ) THEN + evar => VariableGet( Mesh % Variables, 'Electric field', ThisOnly = .TRUE.) + END IF + IF(.NOT. ASSOCIATED(evar) ) THEN + CALL Fatal(Caller,'Could not find nodal electric field!') + END IF + + IntNorm = 0.0_dp + IntTan = 0.0_dp + + ! integration + Active = GetNOFBoundaryElements() + DO i=1,Active + Element => GetBoundaryElement(i) + + iMode = GetBCId( Element ) + IF ( iMode > 0 ) THEN + BC => Model % BCs(iMode) % Values + + IF (.NOT. ASSOCIATED(BC) ) CYCLE + + Found = .FALSE. + bc_name = ListGetString(BC, 'Boundary Energy Name', Found ) + + IF( Found ) THEN + CALL LocalIntegBC(Element) + END IF + END IF + END DO + + ! parallel sum + IF( ParEnv % PEs > 1 ) THEN + DO i=1,Model % NumberOfBCs + IntNorm(i) = ParallelReduction(IntNorm(i)) + IntTan(i) = ParallelReduction(IntTan(i)) + END DO + END IF + + IntNorm = 0.5 * IntNorm + IntTan = 0.5 * IntTan + + ! save data to be used from other solvers + DO i=1,Model % NumberOfBCs + BC => Model % BCs(i) % Values + + IF (.NOT. ASSOCIATED(BC) ) CYCLE + + Found = .FALSE. + bc_name = ListGetString(BC,'Boundary Energy Name', Found ) + IF( Found ) THEN + ! Make data available for SaveData / SaveScalars solver + CALL ListAddConstReal(GetSimulation(), TRIM(bc_name) // "_norm_component", IntNorm(i)) + CALL ListAddConstReal(GetSimulation(), TRIM(bc_name) // "_tan_component", IntTan(i)) + END IF + END DO + +CONTAINS + !----------------------------------------------------------------------------- + SUBROUTINE LocalIntegBC(Element) + !------------------------------------------------------------------------------ + TYPE(Element_t), POINTER :: Element + !------------------------------------------------------------------------------ + REAL(KIND=dp) :: e_ip(3), e_ip_norm, e_ip_tan(3) + REAL(KIND=dp), ALLOCATABLE :: Basis(:), e_local(:,:) + REAL(KIND=dp), ALLOCATABLE :: e_local_parent(:,:), e_local_parentL(:,:), e_local_parentR(:,:) + REAL(KIND=dp) :: weight, DetJ, Normal(3) + LOGICAL :: Stat, Found, E_R_zero, E_L_zero + TYPE(GaussIntegrationPoints_t) :: IP + INTEGER :: t, i, m, ndofs, n, parent_node_ind, boundary_node_ind + TYPE(Nodes_t), SAVE :: Nodes + LOGICAL :: AllocationsDone = .FALSE. + TYPE(Nodes_t) :: ParentNodesR, ParentNodesL + TYPE(Element_t), POINTER :: ParentElementR, ParentElementL + TYPE(ValueList_t), POINTER :: ParentMaterialR, ParentMaterialL + REAL(KIND=dp) :: epsr_R, epsr_L + SAVE AllocationsDone, Basis, e_local, e_local_parent, e_local_parentR, e_local_parentL + REAL(KIND=dp), PARAMETER :: zero_threshold = 1e-25 ! TODO find a good threshold + + ndofs = evar % dofs + IF(.NOT. AllocationsDone ) THEN + m = Mesh % MaxElementDOFs + ALLOCATE( Basis(m), e_local(ndofs,m), e_local_parent(ndofs,m), & + e_local_parentR(ndofs,m), e_local_parentL(ndofs,m)) + AllocationsDone = .TRUE. + END IF + + e_local = 0.0_dp + e_local_parentR = 0.0_dp + e_local_parentL = 0.0_dp + + ParentElementR => Element % BoundaryInfo % Right + ParentElementL => Element % BoundaryInfo % Left + + IF ( .NOT. (ASSOCIATED(ParentElementR) .AND. ASSOCIATED(ParentElementL))) THEN + CALL FATAL(Caller,'No parent elements found') + END IF + + CALL GetElementNodes( Nodes, Element ) + n = Element % TYPE % NumberOfNodes + + CALL GetElementNodes( ParentNodesR, ParentElementR ) + CALL GetElementNodes( ParentNodesL, ParentElementL ) + + ParentMaterialR => GetMaterial(ParentElementR) + ParentMaterialL => GetMaterial(ParentElementL) + + Found = .FALSE. + epsr_R = GetConstReal(ParentMaterialR,'Relative Permittivity', Found) + IF (.NOT. FOUND) THEN + CALL FATAL(Caller, 'Did not find R parent permittivity') + END IF + + Found = .FALSE. + epsr_L = GetConstReal(ParentMaterialL,'Relative Permittivity', Found) + IF (.NOT. FOUND) THEN + CALL FATAL(Caller, 'Did not find L parent permittivity') + END IF + + e_local_parent = 0.0_dp + Found = .FALSE. + CALL GetVectorLocalSolution( e_local_parent, uelement = ParentElementR, uvariable = evar, Found=Found) + IF (.NOT. FOUND) THEN + CALL FATAL(Caller, 'Did not find ParentR e-field') + END IF + + ! Extract e-field from parent element common nodes on the boundary + outer_loop_R: DO boundary_node_ind=1,n + Found = .FALSE. + inner_loop_R: DO parent_node_ind=1, ParentElementR % TYPE % NumberOfNodes + IF ( ParentElementR % NodeIndexes(parent_node_ind) == Element % NodeIndexes(boundary_node_ind) ) THEN + e_local_parentR(:, boundary_node_ind) = e_local_parent(:, parent_node_ind) + Found = .TRUE. + EXIT inner_loop_R + END IF + END DO inner_loop_R + IF (.NOT. Found) THEN + CALL FATAL(Caller, 'Did not find a matching node on the parent element (R)') + END IF + END DO outer_loop_R + + e_local_parent = 0.0_dp + Found = .FALSE. + CALL GetVectorLocalSolution( e_local_parent, uelement = ParentElementL, uvariable = evar, Found=Found) + IF (.NOT. FOUND) THEN + CALL FATAL(Caller, 'Did not find ParentL e-field') + END IF + + outer_loop_L: DO boundary_node_ind=1,n + Found = .FALSE. + inner_loop_L: DO parent_node_ind=1, ParentElementL % TYPE % NumberOfNodes + IF ( ParentElementL % NodeIndexes(parent_node_ind) == Element % NodeIndexes(boundary_node_ind) ) THEN + e_local_parentL(:, boundary_node_ind) = e_local_parent(:, parent_node_ind) + Found = .TRUE. + EXIT inner_loop_L + END IF + END DO inner_loop_L + IF (.NOT. Found) THEN + CALL FATAL(Caller, 'Did not find a matching node on the parent element (L)') + END IF + END DO outer_loop_L + + E_R_zero = SUM(ABS(e_local_parentR(:, 1:n))) < zero_threshold + E_L_zero = SUM(ABS(e_local_parentL(:, 1:n))) < zero_threshold + + IF ( E_R_zero ) THEN + IF ( E_L_zero ) THEN + e_local(:, 1:n) = 0.0_dp + ELSE + ! choose E_L + e_local(:, 1:n) = e_local_parentL(:, 1:n) + END IF + ELSE + IF ( E_L_zero ) THEN + ! choose E_R + e_local(:, 1:n) = e_local_parentR(:, 1:n) + ELSE + ! choose the one with lower permittivity + IF (epsr_R < epsr_L) THEN + e_local(:, 1:n) = e_local_parentR(:, 1:n) + END IF + e_local(:, 1:n) = e_local_parentL(:, 1:n) + END IF + END IF + + ! Sign of normal vector does not matter as we are only interested in the absolute flux through each element + Normal = NormalVector(Element, Nodes, Check=.TRUE.) + + ! Numerical integration: + !----------------------- + IP = GaussPoints(Element) + DO t=1,IP % n + stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), & + IP % W(t), detJ, Basis) + weight = IP % s(t) * detJ + + DO i=1,3 + e_ip(i) = SUM( Basis(1:n) * e_local(i,1:n) ) + END DO + + e_ip_norm = SUM(e_ip*Normal) + e_ip_tan = e_ip - e_ip_norm * Normal + + IntNorm(iMode) = IntNorm(iMode) + weight * e_ip_norm*e_ip_norm + IntTan(iMode) = IntTan(iMode) + weight * SUM(e_ip_tan*e_ip_tan) + END DO + + !------------------------------------------------------------------------------ + END SUBROUTINE LocalIntegBC + !------------------------------------------------------------------------------ + !------------------------------------------------------------------------ +END SUBROUTINE SaveBoundaryEnergyComponents + !------------------------------------------------------------------------ + diff --git a/klayout_package/python/scripts/simulations/tls_waveguide_sim_elmer.py b/klayout_package/python/scripts/simulations/tls_waveguide_sim_elmer.py new file mode 100644 index 000000000..5b26092b0 --- /dev/null +++ b/klayout_package/python/scripts/simulations/tls_waveguide_sim_elmer.py @@ -0,0 +1,134 @@ +# This code is part of KQCircuits +# Copyright (C) 2024 IQM Finland Oy +# +# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public +# License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied +# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with this program. If not, see +# https://www.gnu.org/licenses/gpl-3.0.html. +# +# The software distribution should follow IQM trademark policy for open-source software +# (meetiqm.com/iqm-open-source-trademark-policy). IQM welcomes contributions to the code. +# Please see our contribution agreements for individuals (meetiqm.com/iqm-individual-contributor-license-agreement) +# and organizations (meetiqm.com/iqm-organization-contributor-license-agreement). +import ast +import logging +import sys + +# from itertools import product +from pathlib import Path + +from kqcircuits.elements.waveguide_coplanar import WaveguideCoplanar +from kqcircuits.pya_resolver import pya +from kqcircuits.simulations.export.elmer.elmer_export import export_elmer + +# from kqcircuits.simulations.export.elmer.elmer_solution import ElmerEPR3DSolution + +from kqcircuits.simulations.export.simulation_export import export_simulation_oas, sweep_simulation +from kqcircuits.simulations.port import EdgePort +from kqcircuits.simulations.post_process import PostProcess +from kqcircuits.simulations.simulation import Simulation +from kqcircuits.util.export_helper import ( + create_or_empty_tmp_directory, + get_active_or_new_layout, + open_with_klayout_or_default_application, +) + + +class TlsWaveguideSim(Simulation): + """A very short segment of waveguide.""" + + def build(self): + self.insert_cell( + WaveguideCoplanar, + path=pya.DPath( + [pya.DPoint(self.box.left, self.box.center().y), pya.DPoint(self.box.right, self.box.center().y)], 0 + ), + ) + self.ports.append(EdgePort(1, pya.DPoint(self.box.left, self.box.center().y), face=0)) + self.ports.append(EdgePort(2, pya.DPoint(self.box.right, self.box.center().y), face=0)) + + +# Prepare output directory +dir_path = create_or_empty_tmp_directory(Path(__file__).stem + "_output") + + +# If True models the tls layers as sheets. Leads to computationally easier system especially with +# small interface thicknesses. sheet_approximations need to be set in post-processing to get correct EPRs. +# Uses a custom Elmer energy integration module +sheet_interfaces = True + + +dielectric_surfaces = { + "MA": {"thickness": 1e-8, "eps_r": 8, "background_eps_r": 1.0}, + "SA": {"thickness": 1e-8, "eps_r": 4, "background_eps_r": 11.45}, + "MS": {"thickness": 1e-8, "eps_r": 11.4, "background_eps_r": 11.45}, +} + +# Simulation parameters +sim_class = TlsWaveguideSim # pylint: disable=invalid-name +sim_parameters = { + "name": "tls_waveguide_sim", + "face_stack": ["1t1"], # single chip + "box": pya.DBox(pya.DPoint(0, 0), pya.DPoint(10, 100)), + "substrate_height": 50, # limited simulation domain + "upper_box_height": 50, # limited simulation domain + "metal_height": [0.2], + "partition_regions": [ + {"name": "mer", "metal_edge_dimensions": 1.0, "vertical_dimensions": 1.0, "face": "1t1", "visualise": True} + ], + "tls_layer_thickness": [dielectric_surfaces[layer]["thickness"] * 1e6 for layer in ("MA", "MS", "SA")], + "tls_layer_material": ["oxideMA", "oxideMS", "oxideSA"], + "material_dict": { + **ast.literal_eval(Simulation.material_dict), + "oxideMA": {"permittivity": dielectric_surfaces["MA"]["eps_r"]}, + "oxideMS": {"permittivity": dielectric_surfaces["MS"]["eps_r"]}, + "oxideSA": {"permittivity": dielectric_surfaces["SA"]["eps_r"]}, + }, + "tls_sheet_approximation": sheet_interfaces, + "detach_tls_sheets_from_body": not sheet_interfaces, +} + +sol_parameters = { + "tool": "epr_3d", + "mesh_size": {"1t1_layerMAmer": 0.5, "1t1_layerMSmer": 0.5, "1t1_layerSAmer": 0.5}, + "linear_system_method": "mg", +} + +if sheet_interfaces: + post_process = PostProcess( + "produce_epr_table.py", + sheet_approximations=dielectric_surfaces, + ) +else: + # Refine MA wall if using 3D interfaces + sol_parameters["mesh_size"]["1t1_layerMAwallmer"] = 0.3 + + post_process = PostProcess("produce_epr_table.py") + + +workflow = { + "run_gmsh_gui": False, + "run_elmergrid": True, + "run_elmer": True, + "run_paraview": False, + "python_executable": "python", + "gmsh_n_threads": -1, # Number of omp threads in gmsh + "elmer_n_processes": -1, # Number of dependent processes (tasks) in elmer + "elmer_n_threads": 1, # Number of omp threads per process in elmer +} + +# Get layout +logging.basicConfig(level=logging.WARN, stream=sys.stdout) +layout = get_active_or_new_layout() + +simulations = sweep_simulation(layout, sim_class, sim_parameters, {"a": [2, 10]}) + +export_elmer(simulations, path=dir_path, workflow=workflow, post_process=post_process, **sol_parameters) + +# Write and open oas file +open_with_klayout_or_default_application(export_simulation_oas(simulations, dir_path)) diff --git a/singularity/create_links.sh b/singularity/create_links.sh index d8db2e3ee..4b84f3d5e 100755 --- a/singularity/create_links.sh +++ b/singularity/create_links.sh @@ -4,7 +4,7 @@ function create_link () { ln -sfv "$PWD"/libexec/kqclib.sh bin/"$1" } -EXECUTABLES=("ElmerSolver" "ElmerSolver_mpi" "ElmerGrid" "klayout" "kqclib" "paraview" "python") +EXECUTABLES=("ElmerSolver" "ElmerSolver_mpi" "ElmerGrid" "klayout" "kqclib" "paraview" "python" "elmerf90") echo "Creating symbolic links to the singularity image software" mkdir -p bin