From 31ae1245a7eef7761d4d63f0b7cbe27a5c180bda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Mon, 27 Oct 2025 11:23:37 +0100 Subject: [PATCH 1/4] Add semi-abstract class for linear simulations Update importing system to use __all__. Update basic example. Charge calculation moved out of simulation functions. Moved enums to separate module. --- plutho/__init__.py | 19 +- plutho/coupled_sim.py | 297 ++++--- plutho/enums.py | 37 + plutho/io.py | 82 +- plutho/materials.py | 84 +- plutho/mesh/mesh.py | 9 + plutho/postprocessing.py | 16 +- plutho/simulation/__init__.py | 7 - plutho/simulation/nonlinear/__init__.py | 3 - plutho/simulation/nonlinear/base.py | 410 --------- plutho/simulation/nonlinear/piezo_hb.py | 505 ----------- plutho/simulation/nonlinear/piezo_time.py | 410 --------- plutho/simulation/piezo_time.py | 479 ---------- plutho/simulations/__init__.py | 5 + plutho/simulations/helpers.py | 224 +++++ .../base.py => simulations/integrals.py} | 701 +++++++-------- .../{simulation => simulations}/piezo_freq.py | 242 ++---- plutho/simulations/piezo_time.py | 247 ++++++ plutho/simulations/solver.py | 394 +++++++++ .../thermo_piezo_time.py | 284 ++---- .../thermo_time.py | 392 ++++----- plutho/single_sim.py | 818 ------------------ scripts/basic_example.py | 136 ++- tests/test_fields.py | 131 ++- 24 files changed, 1911 insertions(+), 4021 deletions(-) create mode 100644 plutho/enums.py delete mode 100644 plutho/simulation/__init__.py delete mode 100644 plutho/simulation/nonlinear/__init__.py delete mode 100644 plutho/simulation/nonlinear/base.py delete mode 100644 plutho/simulation/nonlinear/piezo_hb.py delete mode 100644 plutho/simulation/nonlinear/piezo_time.py delete mode 100644 plutho/simulation/piezo_time.py create mode 100644 plutho/simulations/__init__.py create mode 100644 plutho/simulations/helpers.py rename plutho/{simulation/base.py => simulations/integrals.py} (56%) rename plutho/{simulation => simulations}/piezo_freq.py (50%) create mode 100644 plutho/simulations/piezo_time.py create mode 100644 plutho/simulations/solver.py rename plutho/{simulation => simulations}/thermo_piezo_time.py (57%) rename plutho/{simulation => simulations}/thermo_time.py (51%) delete mode 100644 plutho/single_sim.py diff --git a/plutho/__init__.py b/plutho/__init__.py index 861fa9d..f6fba04 100644 --- a/plutho/__init__.py +++ b/plutho/__init__.py @@ -1,12 +1,7 @@ -from .simulation import PiezoSimTime, ThermoPiezoSimTime, MeshData, \ - SimulationData, MaterialData, SimulationType, ModelType, \ - ThermoSimTime, NonlinearType, NLPiezoHB, NLPiezoTime, \ - Nonlinearity -from .io import parse_charge_hist_file, parse_displacement_hist_file, \ - create_vector_field_as_csv, create_scalar_field_as_csv -from .postprocessing import calculate_impedance, \ - calculate_electrical_input_energy, calculate_stored_thermal_energy -from .materials import MaterialManager -from .single_sim import SingleSimulation, FieldType -from .mesh import Mesh -from .coupled_sim import CoupledThermoPiezoThermoSim, CoupledFreqPiezoTherm +from .simulations import * +from .mesh import * +from .io import * +from .postprocessing import * +from .materials import * +from .enums import * +from .coupled_sim import * diff --git a/plutho/coupled_sim.py b/plutho/coupled_sim.py index 14a5135..b8d7bd3 100644 --- a/plutho/coupled_sim.py +++ b/plutho/coupled_sim.py @@ -1,20 +1,36 @@ """Module for coupled simulations.""" # Python stanard libraries -from typing import Dict, Union +from dataclasses import dataclass # Third party libraries import numpy as np import numpy.typing as npt # Local libraries -# TODO Make relative -from plutho import SingleSimulation, SimulationData, Mesh, MaterialData, \ - FieldType -from plutho.simulation.base import get_avg_temp_field_per_element +from .simulations import ThermoPiezoTime, ThermoTime, PiezoFreq +from .simulations.helpers import get_avg_temp_field_per_element +from .mesh.mesh import Mesh +from .materials import MaterialData +from .enums import FieldType -class CoupledThermoPiezoThermoSim: +__all__ = [ + "SimulationData", + "CoupledThermoPiezoTime", + "CoupledPiezoThermoFreq" +] + + +@dataclass +class SimulationData: + delta_t: float + number_of_time_steps: int + gamma: float + beta: float + + +class CoupledThermoPiezoTime: """Couples a thermo-piezoelectric simulation with a heat conduction simulation where only a temperature field is calculated. The power losses from the thermo-piezoeletric simulation are assumed to be @@ -24,58 +40,35 @@ class CoupledThermoPiezoThermoSim: Also a 'Symaxis' for the symmetric axis property is needed. Attributes: - working_directory: Path to the directory where the simulation - folder is created. simulation_name: Name of the simulation and the simulation folder. piezo_sim: Thermo-piezoeletric simulation. heat_cond_sim: Heat conduction simulation. """ # Basic settings - working_directory: str simulation_name: str # Simulations - thermo_piezo_sim: SingleSimulation - thermo_sim: SingleSimulation + thermo_piezo_sim: ThermoPiezoTime + thermo_sim: ThermoTime def __init__( self, - working_directory: str, simulation_name: str, - thermo_piezo_simulation_data: SimulationData, - thermo_simulation_data: SimulationData, mesh: Mesh ): - self.working_directory = working_directory self.simulation_name = simulation_name # Init simulations - thermo_piezo_sim = SingleSimulation( - working_directory, + thermo_piezo_sim = ThermoPiezoTime( simulation_name, mesh ) - # TODO Unecessary to create a new simulation data object - # in the function - thermo_piezo_sim.setup_thermo_piezo_time_domain( - thermo_piezo_simulation_data.number_of_time_steps, - thermo_piezo_simulation_data.delta_t, - thermo_piezo_simulation_data.gamma, - thermo_piezo_simulation_data.beta - ) - - thermo_sim = SingleSimulation( - working_directory, + thermo_sim = ThermoTime( simulation_name, mesh ) - thermo_sim.setup_thermo_time_domain( - thermo_simulation_data.delta_t, - thermo_simulation_data.number_of_time_steps, - thermo_simulation_data.gamma - ) self.thermo_piezo_sim = thermo_piezo_sim self.thermo_sim = thermo_sim @@ -84,7 +77,7 @@ def add_material( self, material_name: str, material_data: MaterialData, - physical_group_name: Union[str, None] + physical_group_name: str ): """Add a material to the simulation. @@ -162,17 +155,22 @@ def set_convection_bc( ] ) - self.thermo_sim.solver.set_convection_bc( + self.thermo_sim.set_convection_bc( convective_bc_elements, alpha, outer_temperature ) + def assemble(self): + """Assembles both underlying simulations.""" + self.thermo_piezo_sim.assemble() + self.thermo_sim.assemble() + def simulate( self, - averaging_time_step_count: int = 100, - thermo_piezo_kwargs: Dict = {}, - thermo_kwargs: Dict = {} + thermo_piezo_simulation_data: SimulationData, + thermo_simulation_data: SimulationData, + averaging_time_step_count: int = 100 ): """Runs the coupled simulation. First the thermo-piezoelectric simulation is run. Then the power losses in the stationary state @@ -181,91 +179,103 @@ def simulate( is done. Parameters: + thermo_piezo_simulation_data: Simulation data object, which + contains the settings for the thermo-piezo simulation. + thermo_simulation_data: Simulation data object, which contains + the settings for the thermo simulation. averaging_time_step_count: Number of time steps starting from the last time step to calculate the average power loss. - thermo_piezo_kwargs: Parameters given to the - SingleSimulation.simulate() function for the thermo piezo - simulation. - thermo_kwargs: Parameters given to the SingleSimulation.simulate() - function for the thermo simulation. - """ # Run thermo_piezoelectric simulation first - self.thermo_piezo_sim.simulate(**thermo_piezo_kwargs) + self.thermo_piezo_sim.simulate( + thermo_piezo_simulation_data.delta_t, + thermo_piezo_simulation_data.number_of_time_steps, + thermo_piezo_simulation_data.gamma, + thermo_piezo_simulation_data.beta, + ) # Calculate mech losses which are sources in heat cond sim avg_mech_loss = np.mean( - self.thermo_piezo_sim.solver.mech_loss[ + self.thermo_piezo_sim.mech_loss[ :, -averaging_time_step_count: ], axis=1 ) - self.thermo_sim.solver.set_constant_volume_heat_source( + self.thermo_sim.set_constant_volume_heat_source( avg_mech_loss, - self.thermo_sim.solver.simulation_data.number_of_time_steps + thermo_simulation_data.number_of_time_steps ) # Set the result temp field of piezo sim to start field of heat cond # sim number_of_nodes = len(self.thermo_piezo_sim.mesh_data.nodes) - theta_start = self.thermo_piezo_sim.solver.u[3*number_of_nodes:, -1] + theta_start = self.thermo_piezo_sim.u[-1, 3*number_of_nodes:] self.thermo_sim.simulate( - initial_theta_field=theta_start, - **thermo_kwargs + thermo_simulation_data.delta_t, + thermo_simulation_data.number_of_time_steps, + thermo_simulation_data.gamma, + theta_start + ) + + def save_simulation_settings(self, directory: str): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation settings in it. + + Parameters: + directory: Path to the directory in which the simulation folder is + created. + """ + self.thermo_piezo_sim.save_simulation_settings( + directory, + "thermo_piezo" ) + self.thermo_sim.save_simulation_settings(directory, "heat_cond") - def save_simulation_results(self): - """Saves the simulation results to the simulation folder.""" - self.thermo_piezo_sim.save_simulation_results("thermo_piezo") - self.thermo_sim.save_simulation_results("heat_cond") + def save_simulation_results(self, directory: str): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation results in it. - def save_simulation_settings(self): - """Saves the simulation settings to the simulation folder.""" - self.thermo_piezo_sim.save_simulation_settings("thermo_piezo") - self.thermo_sim.save_simulation_settings("heat_cond") + Parameters: + directory: Path to the directory in which the simulation folder is + created. + """ + self.thermo_piezo_sim.save_simulation_results( + directory, + "thermo_piezo" + ) + self.thermo_sim.save_simulation_results( + directory, + "heat_cond" + ) -class CoupledFreqPiezoTherm: +class CoupledPiezoThermoFreq: - # Basic setup - working_directory: str simulation_name: str # Simulations - piezo_freq: SingleSimulation - thermo_time_sim: SingleSimulation + piezo_freq: PiezoFreq + thermo_time_sim: ThermoTime def __init__( self, - working_directory: str, simulation_name: str, - piezo_freq_frequency: float, - thermo_time_simulation_data: SimulationData, mesh: Mesh ): - self.working_directory = working_directory self.simulation_name = simulation_name - piezo_freq = SingleSimulation( - working_directory, + piezo_freq = PiezoFreq( simulation_name, mesh ) - piezo_freq.setup_piezo_freq_domain(np.array([piezo_freq_frequency])) - thermo_time_sim = SingleSimulation( - working_directory, + thermo_time_sim = ThermoTime( simulation_name, mesh ) - thermo_time_sim.setup_thermo_time_domain( - thermo_time_simulation_data.delta_t, - thermo_time_simulation_data.number_of_time_steps, - thermo_time_simulation_data.gamma - ) self.piezo_freq = piezo_freq self.thermo_time_sim = thermo_time_sim @@ -274,7 +284,7 @@ def add_material( self, material_name: str, material_data: MaterialData, - physical_group_name: Union[str, None] + physical_group_name: str ): """Add a material to the simulation. @@ -325,8 +335,15 @@ def set_excitation( np.zeros(1) ) + def assemble(self): + """Assembles the underlying simulations.""" + self.piezo_freq.assemble() + self.thermo_time_sim.assemble() + def simulate( self, + piezo_freq_frequency: float, + thermo_time_simulation_data: SimulationData, material_starting_temperature: float, temperature_dependent: bool ): @@ -334,26 +351,36 @@ def simulate( simulation classes. Parameters: - material_starting_temperature: Sets the temperature at time step 0. - Sets the starting temperature field as well as the material + piezo_freq_frequency: Frequency for which the frqeuency domain + simulation is solved. + thermo_time_simulation_data: Simulation data object for the + thermo time simulation. + material_starting_temperature: Sets the temperature at time step + 0. Sets the starting temperature field as well as the material parameters at the start. temperature_dependent: Set to true if the materials shall be - temperature dependent (Changing material parameters based on the - current field temperature). + temperature dependent (Changing material parameters based on + the current field temperature). """ time_index = 0 - number_of_time_steps = \ - self.thermo_time_sim.solver.simulation_data.number_of_time_steps - number_of_nodes = len(self.thermo_time_sim.solver.mesh_data.nodes) + number_of_time_steps = thermo_time_simulation_data \ + .number_of_time_steps + number_of_nodes = len(self.thermo_time_sim.mesh_data.nodes) + number_of_elements = len(self.thermo_time_sim.mesh_data.elements) if temperature_dependent: temp_field_per_element = None while time_index < number_of_time_steps: # Run piezo_freq simulation if time_index == 0: + self.piezo_freq.material_manager.update_temperature( + material_starting_temperature * np.ones( + number_of_elements + ) + ) self.piezo_freq.simulate( - calculate_mech_loss=True, - material_starting_temperature=material_starting_temperature + np.array([piezo_freq_frequency]), + True ) else: if temp_field_per_element is None: @@ -366,17 +393,18 @@ def simulate( temp_field_per_element ) self.piezo_freq.simulate( - calculate_mech_loss=True, + np.array([piezo_freq_frequency]), + True, ) # Get mech losses mech_loss = np.real( - self.piezo_freq.solver.mech_loss[:, 0] + self.piezo_freq.mech_loss[time_index, :] ) # Use mech losses for heat conduction simulation - self.thermo_time_sim.solver.f *= 0 # Reset load vector - self.thermo_time_sim.solver.set_constant_volume_heat_source( + self.thermo_time_sim.f *= 0 # Reset load vector + self.thermo_time_sim.set_constant_volume_heat_source( mech_loss, number_of_time_steps ) @@ -385,19 +413,33 @@ def simulate( # changed enough the simulation stops and a frequency domain # simulation is done. if time_index == 0: - time_index = self.thermo_time_sim.simulate( - initial_theta_time_step=0, - initial_theta_field=( - material_starting_temperature - * np.ones(number_of_nodes) - ), - material_starting_temperature=material_starting_temperature + self.thermo_time_sim.material_manager.update_temperature( + material_starting_temperature*np.ones( + number_of_elements) ) - else: - time_index = self.thermo_time_sim.simulate( - initial_theta_time_step=time_index, - material_starting_temperature=material_starting_temperature + time_index = self.thermo_time_sim \ + .simulate_until_material_parameters_change( + thermo_time_simulation_data.delta_t, + thermo_time_simulation_data.number_of_time_steps, + thermo_time_simulation_data.gamma, + ( + material_starting_temperature + * np.ones(number_of_nodes) + ), + 0 ) + else: + time_index = self.thermo_time_sim \ + .simulate_until_material_parameters_change( + thermo_time_simulation_data.delta_t, + thermo_time_simulation_data.number_of_time_steps, + thermo_time_simulation_data.gamma, + ( + material_starting_temperature + * np.ones(number_of_nodes) + ), + time_index, + ) if time_index == number_of_time_steps: # Simulation is done @@ -410,37 +452,52 @@ def simulate( f"Updating material parameters at time step {time_index}" ) temp_field_per_element = get_avg_temp_field_per_element( - self.thermo_time_sim.solver.theta[:, time_index], - self.thermo_time_sim.solver.mesh_data.elements + self.thermo_time_sim.theta[time_index, :], + self.thermo_time_sim.mesh_data.elements ) else: # Run piezofreq simulation self.piezo_freq.simulate( - calculate_mech_loss=True + np.array([piezo_freq_frequency]), + True ) # Get mech losses - mech_loss = np.real(self.piezo_freq.solver.mech_loss[:, 0]) + mech_loss = np.real(self.piezo_freq.mech_loss[time_index, :]) # Use mech losses for heat conduction simulation - self.thermo_time_sim.solver.set_constant_volume_heat_source( + self.thermo_time_sim.set_constant_volume_heat_source( mech_loss, number_of_time_steps ) # Run thermo time domain simulation - self.thermo_time_sim.simulate( - initial_theta_field=( - material_starting_temperature*np.ones(number_of_nodes) - ) + self.thermo_time_sim.simulate_until_material_parameters_change( + thermo_time_simulation_data.delta_t, + thermo_time_simulation_data.number_of_time_steps, + thermo_time_simulation_data.gamma, + material_starting_temperature*np.ones(number_of_nodes), + 0 ) - def save_simulation_results(self): - """Saves the simulation results to the simulation folder.""" - self.piezo_freq.save_simulation_results("piezo_freq") - self.thermo_time_sim.save_simulation_results("heat_cond") + def save_simulation_settings(self, directory: str): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation settings in it. - def save_simulation_settings(self): - """Saves the simulation settings to the simulation folder.""" - self.piezo_freq.save_simulation_settings("piezo_freq") - self.thermo_time_sim.save_simulation_settings("heat_cond") + Parameters: + directory: Path to the directory in which the simulation folder is + created. + """ + self.piezo_freq.save_simulation_settings(directory, "piezo_freq") + self.thermo_time_sim.save_simulation_settings(directory, "heat_cond") + + def save_simulation_results(self, directory: str): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation results in it. + + Parameters: + directory: Path to the directory in which the simulation folder is + created. + """ + self.piezo_freq.save_simulation_results(directory, "piezo_freq") + self.thermo_time_sim.save_simulation_results(directory, "heat_cond") diff --git a/plutho/enums.py b/plutho/enums.py new file mode 100644 index 0000000..3072d7a --- /dev/null +++ b/plutho/enums.py @@ -0,0 +1,37 @@ +"""Module for enums.""" + + +# Python standard library +from enum import Enum + + +__all__ = [ + "SolverType", + "NonlinearType", + "FieldType" +] + + +class SolverType(Enum): + """Contains the possible solver types.""" + PiezoFreq = "PiezoFreq", + PiezoTime = "PiezoTime", + ThermoPiezoTime = "ThermoPiezoTime", + ThermoTime = "ThermoTime" + + +class NonlinearType(Enum): + """Contains the various possible nonlinear simulation types.""" + QuadraticRayleigh = 0, + QuadraticCustom = 1, + CubicRayleigh = 2, + CubicCustom = 3 + + +class FieldType(Enum): + """Possible field types which are calculated using differnet simulations. + """ + U_R = 0 + U_Z = 1 + PHI = 2 + THETA = 3 diff --git a/plutho/io.py b/plutho/io.py index 27beb7b..73304a1 100644 --- a/plutho/io.py +++ b/plutho/io.py @@ -1,7 +1,6 @@ """Module for io functions for the resulting fields and values.""" # Python standard libraries -from typing import Tuple import os # Third party libraries @@ -10,83 +9,14 @@ import yaml # Local libraries -from .simulation.base import MaterialData +from .materials import MaterialData -def parse_charge_hist_file(file_path: str) -> Tuple[npt.NDArray, npt.NDArray]: - """Reads the charge file from an OpenCFS simulation. - - Parameters: - file_path: Path to the charge (*.hist) file. - - Returns: - Tuple containing the time list and charge list read from the - given file. - """ - lines = [] - with open(file_path, "r", encoding="UTF-8") as fd: - lines = fd.readlines()[3:] - - time = [] - charge = [] - for line in lines: - current_time, current_charge = line.split() - time.append(float(current_time)) - charge.append(float(current_charge)) - - return np.array(time), np.array(charge) - - -def parse_charge_freq_hist_file( - file_path: str -) -> Tuple[npt.NDArray, npt.NDArray]: - """Reads the charge file from an OpenCFS simulation. - - Parameters: - file_path: Path to the charge (*.hist) file. - - Returns: - Tuple containing the time list and charge list read from the - given file. - """ - lines = [] - with open(file_path, "r", encoding="UTF-8") as fd: - lines = fd.readlines()[3:] - - frequencies = [] - charge = [] - for line in lines: - current_frequency, current_charge, _ = line.split() - frequencies.append(float(current_frequency)) - charge.append(float(current_charge)) - - return np.array(frequencies), np.array(charge) - - -def parse_displacement_hist_file( - file_path: str -) -> Tuple[npt.NDArray, npt.NDArray, npt.NDArray]: - """Reads the displacement file from and OpenCFS simulation. - - Parameters: - file_path: Path to the displacement (*.hist) file. - - Returns: - Tuple containg the time steps, u_r and u_z. - """ - lines = [] - with open(file_path, "r", encoding="UTF-8") as fd: - lines = fd.readlines() - time_steps = [] - u_r = [] - u_z = [] - for _, line in enumerate(lines[3:]): - current_time_step, current_u_r, current_u_z = line.split() - time_steps.append(float(current_time_step)) - u_r.append(float(current_u_r)) - u_z.append(float(current_u_z)) - - return np.array(time_steps), np.array(u_r), np.array(u_z) +__all__ = [ + "create_scalar_field_as_csv", + "create_vector_field_as_csv", + "load_temperature_dependent_material_data_pic181" +] def create_scalar_field_as_csv( diff --git a/plutho/materials.py b/plutho/materials.py index d41fd70..021a076 100644 --- a/plutho/materials.py +++ b/plutho/materials.py @@ -1,20 +1,92 @@ """Contains different materials which can be used in the simulations.""" # Python standard libraries +import json +from dataclasses import dataclass, fields from typing import List, Union -import numpy as np -import numpy.typing as npt # Third party libraries +import numpy as np +import numpy.typing as npt from scipy import interpolate -# Local libraries -from .simulation.base import MaterialData + +__all__ = [ + "MaterialData", + "MaterialManager" +] + + +@dataclass +class MaterialData: + """Contains the plain material data. Some parameters can either be + a float or an array depending if they are temperature dependent. + If they are temperature dependent the index in the array corresponds to + the temperature value from the temperatures array at the same index. + """ + c11: Union[float, npt.NDArray] + c12: Union[float, npt.NDArray] + c13: Union[float, npt.NDArray] + c33: Union[float, npt.NDArray] + c44: Union[float, npt.NDArray] + e15: Union[float, npt.NDArray] + e31: Union[float, npt.NDArray] + e33: Union[float, npt.NDArray] + eps11: Union[float, npt.NDArray] + eps33: Union[float, npt.NDArray] + alpha_m: float + alpha_k: float + thermal_conductivity: float + heat_capacity: float + temperatures: Union[float, npt.NDArray] + density: float + + def to_dict(self): + """Convert the dataclass to dict for json serialization.""" + json_dict = {} + for attribute in fields(self.__class__): + value = getattr(self, attribute.name) + if isinstance(value, float) or isinstance(value, int): + json_dict[attribute.name] = value + elif isinstance(value, np.ndarray): + json_dict[attribute.name] = value.tolist() + else: + raise ValueError( + "Wrong type saved in MaterialData. Value is of type " + f"{type(value)}" + ) + return json_dict + + @staticmethod + def from_dict(contents): + """Convert given dict, e.g. from a json deserialization, to a + MaterialData object.""" + for key, value in contents.items(): + if isinstance(value, List): + contents[key] = np.array(value) + + return MaterialData(**contents) + + @staticmethod + def load_from_file(file_path: str): + """Load the data from given file. + + Parameters: + file_path: Path to the file + """ + if not os.path.exists(file_path): + raise IOError( + "Given file path {} does not exist. Cannot load " + "material data." + ) + + with open(file_path, "r", encoding="UTF-8") as fd: + return MaterialData.from_dict(json.load(fd)) class Material: - """Class to handle a single material. The material can have either constant - material properties or temperature variant properties. + """Class to handle a single material. The material can have either + constant material properties or temperature variant properties. The temperature properties are determined using linear interpolation. Attributes: diff --git a/plutho/mesh/mesh.py b/plutho/mesh/mesh.py index 6de14a0..86eaf16 100644 --- a/plutho/mesh/mesh.py +++ b/plutho/mesh/mesh.py @@ -2,6 +2,7 @@ # Python standard libraries import os +from dataclasses import dataclass from typing import Union, Tuple, Dict, List # Third party libraries @@ -14,6 +15,14 @@ from .custom_mesh_parser import CustomParser +@dataclass +class MeshData: + """Contains the mesh data is used in the simulation.""" + nodes: npt.NDArray + elements: npt.NDArray + element_order: int + + class Mesh: """Class to generate and read meshes. diff --git a/plutho/postprocessing.py b/plutho/postprocessing.py index 1b74d3e..ee6aa6b 100644 --- a/plutho/postprocessing.py +++ b/plutho/postprocessing.py @@ -6,8 +6,14 @@ import numpy.typing as npt # Local libraries -from plutho.simulation.base import energy_integral_theta, \ - gradient_local_shape_functions_2d +from .simulations.integrals import energy_integral_theta + + +__all__ = [ + "calculate_impedance", + "calculate_electrical_input_energy", + "calculate_stored_thermal_energy" +] def calculate_impedance( @@ -38,7 +44,7 @@ def calculate_electrical_input_energy( voltage_excitation: npt.NDArray, charge: npt.NDArray, delta_t: float -): +) -> float: """Calculates the energy of the elctric input. Parameters: @@ -85,7 +91,7 @@ def calculate_stored_thermal_energy( stored_energies = np.zeros(theta.shape[1]) for time_index in range(theta.shape[1]): - for element_index, element in enumerate(elements): + for _, element in enumerate(elements): node_points = np.zeros(shape=(2, points_per_element)) for node_index in range(points_per_element): node_points[:, node_index] = [ @@ -112,7 +118,7 @@ def calculate_stored_thermal_energy( # Only one time step stored_energy = 0 - for element_index, element in enumerate(elements): + for _, element in enumerate(elements): node_points = np.zeros(shape=(2, points_per_element)) for node_index in range(points_per_element): node_points[:, node_index] = [ diff --git a/plutho/simulation/__init__.py b/plutho/simulation/__init__.py deleted file mode 100644 index d404de5..0000000 --- a/plutho/simulation/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -from .thermo_piezo_time import ThermoPiezoSimTime -from .piezo_time import PiezoSimTime -from .base import MaterialData, SimulationData, MeshData, ModelType, \ - SimulationType -from .thermo_time import ThermoSimTime -from .piezo_freq import PiezoSimFreq -from .nonlinear import NLPiezoTime, NonlinearType, NLPiezoHB, Nonlinearity diff --git a/plutho/simulation/nonlinear/__init__.py b/plutho/simulation/nonlinear/__init__.py deleted file mode 100644 index 92d6ab1..0000000 --- a/plutho/simulation/nonlinear/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .base import assemble, Nonlinearity, NonlinearType -from .piezo_time import NLPiezoTime -from .piezo_hb import NLPiezoHB diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py deleted file mode 100644 index 3351c40..0000000 --- a/plutho/simulation/nonlinear/base.py +++ /dev/null @@ -1,410 +0,0 @@ -"""Module for base functions used in nonlinear simulations""" - -# Python standard libraries -from enum import Enum -from typing import Tuple - -# Third party libraries -import numpy as np -from scipy import sparse - -# Local libraries -from ..base import MeshData, local_to_global_coordinates, b_operator_global, \ - integral_m, integral_ku, integral_kuv, integral_kve, \ - create_node_points, quadratic_quadrature, gradient_local_shape_functions_2d -from plutho.materials import MaterialManager - -# -# -------- Enums -------- -# -class NonlinearType(Enum): - QuadraticRayleigh = 0, - QuadraticCustom = 1, - CubicRayleigh = 2, - CubicCustom = 3 - -# -# -------- Functions -------- -# -def assemble( - mesh_data: MeshData, - material_manager: MaterialManager -) -> Tuple[ - sparse.lil_array, - sparse.lil_array, - sparse.lil_array -]: - # TODO Assembly takes to long rework this algorithm - # Maybe the 2x2 matrix slicing is not very fast - # TODO Change to lil_array - nodes = mesh_data.nodes - elements = mesh_data.elements - element_order = mesh_data.element_order - node_points = create_node_points(nodes, elements, element_order) - - number_of_nodes = len(nodes) - mu = sparse.lil_matrix( - (2*number_of_nodes, 2*number_of_nodes), - dtype=np.float64 - ) - ku = sparse.lil_matrix( - (2*number_of_nodes, 2*number_of_nodes), - dtype=np.float64 - ) - kuv = sparse.lil_matrix( - (2*number_of_nodes, number_of_nodes), - dtype=np.float64 - ) - kv = sparse.lil_matrix( - (number_of_nodes, number_of_nodes), - dtype=np.float64 - ) - - for element_index, element in enumerate(elements): - current_node_points = node_points[element_index] - - # TODO Check if its necessary to calculate all integrals - # --> Dirichlet nodes could be leaved out? - # Multiply with jac_det because its integrated with respect to - # local coordinates. - # Mutiply with 2*pi because theta is integrated from 0 to 2*pi - mu_e = ( - material_manager.get_density(element_index) - * integral_m(current_node_points, element_order) - * 2 * np.pi - ) - ku_e = ( - integral_ku( - current_node_points, - material_manager.get_elasticity_matrix(element_index), - element_order - ) * 2 * np.pi - ) - kuv_e = ( - integral_kuv( - current_node_points, - material_manager.get_piezo_matrix(element_index), - element_order - ) * 2 * np.pi - ) - kve_e = ( - integral_kve( - current_node_points, - material_manager.get_permittivity_matrix( - element_index - ), - element_order - ) * 2 * np.pi - ) - - # Now assemble all element matrices - for local_p, global_p in enumerate(element): - for local_q, global_q in enumerate(element): - # Mu_e is returned as a 3x3 matrix (should be 6x6) - # because the values are the same. - # Only the diagonal elements have values. - mu[2*global_p, 2*global_q] += mu_e[local_p][local_q] - mu[2*global_p+1, 2*global_q+1] += mu_e[local_p][local_q] - - # Ku_e is a 6x6 matrix and 2x2 matrices are sliced out - # of it. - ku[2*global_p:2*global_p+2, 2*global_q:2*global_q+2] += \ - ku_e[2*local_p:2*local_p+2, 2*local_q:2*local_q+2] - - # KuV_e is a 6x3 matrix and 2x1 vectors are sliced out. - # [:, None] converts to a column vector. - kuv[2*global_p:2*global_p+2, global_q] += \ - kuv_e[2*local_p:2*local_p+2, local_q][:, None] - - # KVe_e is a 3x3 matrix and therefore the values can - # be set directly. - kv[global_p, global_q] += kve_e[local_p, local_q] - - # Calculate damping matrix - # Currently in the simulations only one material is used - # TODO Add algorithm to calculate cu for every element on its own - # in order to use multiple materials - cu = ( - material_manager.get_alpha_m(0) * mu - + material_manager.get_alpha_k(0) * ku - ) - - # Calculate block matrices - zeros1x1 = np.zeros((number_of_nodes, number_of_nodes)) - zeros2x1 = np.zeros((2*number_of_nodes, number_of_nodes)) - zeros1x2 = np.zeros((number_of_nodes, 2*number_of_nodes)) - - m = sparse.block_array([ - [mu, zeros2x1], - [zeros1x2, zeros1x1], - ]).tolil() - c = sparse.block_array([ - [cu, zeros2x1], - [zeros1x2, zeros1x1], - ]).tolil() - k = sparse.block_array([ - [ku, kuv], - [kuv.T, -1*kv] - ]).tolil() - - return m, c, k - - -def integral_u_nonlinear_quadratic( - node_points, - nonlinear_elasticity_tensor, - element_order -): - def inner(s, t): - dn = gradient_local_shape_functions_2d(s, t, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - jacobian_det = np.linalg.det(jacobian) - - b_op = b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - - B_TC = np.einsum("im,mjk->ijk", b_op.T, nonlinear_elasticity_tensor) - B_TCB = np.einsum("ink,nj->ijk", B_TC, b_op) - B_TCBB = np.einsum("ijp,pk->ijk", B_TCB, b_op) - - return B_TCBB * r * jacobian_det - - return quadratic_quadrature(inner, element_order) - -# -# -------- Classes -------- -# -class Nonlinearity: - - def __init__(self): - self.nonlinear_type = None - self.mesh_data = None - self.node_points = None - - def set_mesh_data(self, mesh_data, node_points): - self.mesh_data = mesh_data - self.node_points = node_points - - def set_quadratic_rayleigh(self, zeta): - self.nonlinear_type = NonlinearType.QuadraticRayleigh - self.zeta = zeta - - def set_quadratic_custom(self, nonlinear_data): - self.nonlinear_type = NonlinearType.QuadraticCustom - - # Reduce to axisymmetric 4x4x4 matrix - # TODO Update this for order 3 or make it n dimensional - voigt_map = [0, 1, 3, 2] - nm_reduced = np.zeros(shape=(4, 4, 4)) - for i_new, i_old in enumerate(voigt_map): - for j_new, j_old in enumerate(voigt_map): - for k_new, k_old in enumerate(voigt_map): - nm_reduced[ - i_new, - j_new, - k_new - ] = nonlinear_data[i_old, j_old, k_old] - - self.custom_data = nm_reduced - - def set_cubic_rayleigh(self, zeta): - self.nonlinear_type = NonlinearType.CubicRayleigh - self.zeta = zeta - - def set_cubic_custom(self, nonlinear_data): - self.nonlinear_type = NonlinearType.CubicCustom - - # Reduce to axisymmetric 4x4x4 matrix - # TODO Update this for order 3 or make it n dimensional - voigt_map = [0, 1, 3, 2] - nm_reduced = np.zeros(shape=(4, 4, 4)) - for i_new, i_old in enumerate(voigt_map): - for j_new, j_old in enumerate(voigt_map): - for k_new, k_old in enumerate(voigt_map): - for l_new, l_old in enumerate(voigt_map): - nm_reduced[ - i_new, - j_new, - k_new, - l_new - ] = nonlinear_data[i_old, j_old, k_old, l_old] - - self.custom_data = nm_reduced - - def evaluate_force_vector(self, u, m, c, k): - if self.nonlinear_type is None: - raise ValueError("Cannot evaluate force vector, since no \ - nonlinear type is set") - - match self.nonlinear_type: - case NonlinearType.QuadraticRayleigh: - return self.ln@(u**2) - case NonlinearType.QuadraticCustom: - # Use jit compiler --> wimp - """ - result = np.zeros(3*number_of_nodes) - - for index in range(3*number_of_nodes): - L_k = L[index*3*number_of_nodes:(index+1)*3*number_of_nodes, :] - - result[index] = u @ (L_k @ u) - - return result - """ - raise NotImplementedError() - case NonlinearType.CubicRayleigh: - return self.ln@(u**3) - case NonlinearType.CubicCustom: - raise NotImplementedError() - - def evaluate_jacobian(self, u, m, c, k): - if self.nonlinear_type is None: - raise ValueError("Cannot evaluate jacobian, since no \ - nonlinear type is set") - - match self.nonlinear_type: - case NonlinearType.QuadraticRayleigh: - return 2*self.ln@sparse.diags_array(u) - case NonlinearType.QuadraticCustom: - """ - k_tangent = sparse.lil_matrix(( - 3*number_of_nodes, 3*number_of_nodes - )) - - k_tangent += k - - m = np.arange(3*number_of_nodes) - - # TODO This calculation is very slow - for i in range(3*number_of_nodes): - l_k = ln[3*number_of_nodes*i:3*number_of_nodes*(i+1), :] - l_i = ln[3*number_of_nodes*m+i, :] - - k_tangent += (l_k*u[i]).T - k_tangent += (l_i*u[i]).T - """ - raise NotImplementedError() - case NonlinearType.CubicRayleigh: - return 3*self.ln@sparse.diags_array(u**2) - case NonlinearType.CubicCustom: - raise NotImplementedError() - - def apply_dirichlet_bc(self, dirichlet_nodes): - if self.nonlinear_type is None: - raise ValueError("Cannot apply dirichlet boundary \ - conditions, since no nonlinear type is set") - - match self.nonlinear_type: - case NonlinearType.QuadraticRayleigh | \ - NonlinearType.CubicRayleigh: - self.ln[dirichlet_nodes, :] = 0 - case NonlinearType.QuadraticCustom: - raise NotImplementedError() - case NonlinearType.CubicCustom: - raise NotImplementedError() - - def assemble(self, m, c, k): - if self.nonlinear_type is None: - raise ValueError("Cannot assemble matrices, since no \ - nonlinear type is set") - - match self.nonlinear_type: - case NonlinearType.QuadraticRayleigh | \ - NonlinearType.CubicRayleigh: - self.ln = k*self.zeta - case NonlinearType.QuadraticCustom: - self._assemble_quadratic_custom() - case NonlinearType.CubicCustom: - self._assemble_cubic_custom() - - def _assemble_quadratic_custom(self): - nodes = self.mesh_data.nodes - elements = self.mesh_data.elements - element_order = self.mesh_data.element_order - number_of_nodes = len(nodes) - node_points = create_node_points( - nodes, - elements, - element_order - ) - - lu = sparse.lil_array( - (4*number_of_nodes**2, 2*number_of_nodes), - dtype=np.float64 - ) - - for element_index, element in enumerate(elements): - current_node_points = node_points[element_index] - - lu_e = integral_u_nonlinear_quadratic( - current_node_points, - self.custom_data, - element_order - ) * 2 * np.pi - - for local_i, global_i in enumerate(element): - for local_j, global_j in enumerate(element): - for local_k, global_k in enumerate(element): - # lu_e is 6x6x6 but has to be assembled into an - # 4*N**2x2*N sparse matrix - # The k-th plane is therefore append using the i-th - # index - # Since lu_e is 6x6x6, 2x2x2 slices are taken out and - # used for assembling - k_0 = lu_e[ - 2*local_i:2*(local_i+1), - 2*local_j:2*(local_j+1), - 2*local_k - ] - k_1 = lu_e[ - 2*local_i:2*(local_i+1), - 2*local_j:2*(local_j+1), - 2*local_k+1 - ] - - lu[ - ( - 2*global_i - + 3*number_of_nodes*global_k - ):( - 2*(global_i+1) - + 3*number_of_nodes*global_k - ), - 2*global_j:2*(global_j+1), - ] += k_0 - lu[ - ( - 2*global_i - + 3*number_of_nodes*global_k+1 - ):( - 2*(global_i+1) - + 3*number_of_nodes*global_k+1 - ), - 2*global_j:2*(global_j+1), - ] += k_1 - - ln = sparse.lil_array( - (9*number_of_nodes**2, 3*number_of_nodes), - dtype=np.float64 - ) - - for k in range(3*number_of_nodes): - ln[ - k*3*number_of_nodes:k*3*number_of_nodes+2*number_of_nodes, - :2*number_of_nodes - ] = lu[ - k*2*number_of_nodes:(k+1)*2*number_of_nodes, : - ] - - self.ln = ln - - - def _assemble_cubic_custom(self): - pass diff --git a/plutho/simulation/nonlinear/piezo_hb.py b/plutho/simulation/nonlinear/piezo_hb.py deleted file mode 100644 index 6c098dd..0000000 --- a/plutho/simulation/nonlinear/piezo_hb.py +++ /dev/null @@ -1,505 +0,0 @@ -"""Module for the simulation of nonlinaer piezoelectric systems using the -harmonic balancing method.""" - -# Python standard libraries -from typing import Union, Tuple, Callable - -# Third party libraries -import numpy as np -import numpy.typing as npt -import scipy -from scipy import sparse -from scipy.sparse import linalg - -# Local libraries -from .base import assemble, NonlinearType -from ..base import MeshData -from plutho.materials import MaterialManager - - -class NLPiezoHB: - """Implementes a nonlinear FEM harmonic balancing simulation. - """ - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - - # Boundary conditions - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray - - # FEM matrices - k: sparse.lil_array - c: sparse.lil_array - m: sparse.lil_array - lu: sparse.lil_array - - # Resulting fields - u: npt.NDArray - - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - harmonic_order: int - ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.harmonic_order = harmonic_order - self.dirichlet_nodes = np.array([]) - self.dirichlet_valuse = np.array([]) - - def assemble( - self, - nonlinear_order: int, - nonlinear_type: NonlinearType, - **kwargs - ): - """Redirect to general nonlinear assembly function. - - Parameters: - nonlinear_order: Order of the nonlinearity. - nonlinear_type: Type of the nonlinear material. - **kwargs: Parameters for the nonlinear material. - """ - # Get the default FEM matrices - m, c, k, lu = assemble( - self.mesh_data, - self.material_manager, - nonlinear_order, - nonlinear_type, - **kwargs - ) - self.nonlinear_order = nonlinear_order - - self.m = m - self.c = c - self.k = k - self.lu = lu - - def solve_linear(self): - """Solves the linear problem Ku=f (ln is not used). - """ - # Apply boundary conditions - _, _, k = NLPiezoHB.apply_dirichlet_bc( - None, - None, - self.k.copy(), - self.dirichlet_nodes - ) - - f = self.get_load_vector( - self.dirichlet_nodes, - self.dirichlet_values - ) - - # Calculate u using phi as load vector - self.u = linalg.spsolve( - k, - f - ).todense() - - def solve_newton( - self, - angular_frequency: float, - tolerance: float = 1e-10, - max_iter: int = 100, - u_start: Union[npt.NDArray, None] = None, - alpha: float = 1, - load_factor: float = 1 - ): - """Solves the system of nonlinear equations using the Newton-Raphson - method. Saves the result in self.u - - Parameters: - angular_frequency: Angualr frequency at which the simulation is - done. - tolerance: Upper bound for the residuum. Iteration stops when the - residuum is lower than tolerance. - max_iter: Maximum number of iteration when the tolerance is not - met. Then the best solution (lowest residuum) is returned - u_start: Initial guess for the nonlinear problem. If it is None - a linear solution is calculated and used as a guess for - nonlinear u. - alpha: Damping parameter for updating the guess for u. Can be used - to improve convergence. Choose a value between 0 and 1. - load_factor: Multiplied with the load vector. Used to apply less - than than the full excitation. Could improve conversion of the - algorithm. Choose a value between 0 and 1. - """ - number_of_nodes = len(self.mesh_data.nodes) - - # Get FEM matrices - fem_m = self.m.copy() - fem_c = self.c.copy() - fem_k = self.k.copy() - fem_lu = self.lu.copy() - - # Create the harmonic balancing matrices from the FEM matrices - hb_m, hb_c, hb_k = self.get_harmonic_balancing_matrices( - fem_m, - fem_c, - fem_k, - angular_frequency - ) - - # Apply boundary conditions - hb_m, hb_c, hb_k = NLPiezoHB.apply_dirichlet_bc( - hb_m, - hb_c, - hb_k, - self.dirichlet_nodes - ) - - f = self.get_load_vector( - self.dirichlet_nodes, - self.dirichlet_values - ) - - # Set start value - if u_start is None: - # If no start value is given calculate one using a linear - # simulation - current_u = linalg.spsolve( - ( - hb_m - + hb_c - + hb_k - ), - load_factor*f - ) - else: - current_u = u_start - - # Residual of the Newton algorithm - def residual(u, nonlinear_force): - return ( - hb_m@u - + hb_c@u - + hb_k@u - + nonlinear_force - - load_factor*f - ) - - # Nonlinear forces based on starting field guess - nonlinear_force = self.calculate_nonlinear_force( - current_u, - fem_lu, - number_of_nodes - ) - - # best_field tracks the field with the best norm - # best norm tracks the value of the norm itself - # Calculate the first norms based on the starting field guess - best_field = current_u.copy() - best_norm = np.linalg.norm(residual(best_field, nonlinear_force)) - - # Check if the initial value already suffices the tolerance condition - #if best_norm < tolerance: - # print(f"Initial value is already best {best_norm}") - # self.u = best_field - # return - - # Run Newton method - for iteration_count in range(max_iter): - if iteration_count > 0: - # Does not need to be updated at step 0 - nonlinear_force = self.calculate_nonlinear_force( - current_u, - fem_lu, - number_of_nodes - ) - - # Calculate tangential stiffness matrix - k_tangent = self.calculate_tangent_matrix_analytical( - current_u, - hb_m, - hb_c, - hb_k - ) - - # Solve for displacement increment - delta_u = linalg.spsolve( - k_tangent, - residual(current_u, nonlinear_force) - ) - - # Update step - next_u = current_u - alpha * delta_u - - # Check for convergence - norm = scipy.linalg.norm(residual(next_u, nonlinear_force)) - if norm < tolerance: - print( - f"Newton found solution after {iteration_count+1} " - f"steps with residual {norm}" - ) - self.u = next_u - return - elif norm < best_norm: - best_norm = norm - best_field = next_u.copy() - - if iteration_count % 100 == 0 and iteration_count > 0: - print("Iteration:", iteration_count) - - # Update for next iteration - current_u = next_u - - print("Simulation did not converge") - print("Error from best iteration:", best_norm) - - self.u = best_field - - def calculate_nonlinear_force( - self, - u: npt.NDArray, - fem_lu: sparse.lil_array, - number_of_nodes: int - ) -> npt.NDArray: - """Calculates the forces due to the nonlinearity. - - Parameters: - u: Current solution vector. - fem_lu: FEM nonlinearity matrix. - number_of_nodes: Number of nodes of the mesh. - - Returns: - Vector of nonlinear forces. - """ - # TODO Currently only for third order nonlinearity - - u_c = u[::2] - u_s = u[1::2] - - lu_c = fem_lu.dot( - 3/4*u_c**3+3/4*np.multiply(u_c, u_s**2) - ) - lu_s = fem_lu.dot( - 3/4*u_s**3+3/4*np.multiply(u_s, u_c**2) - ) - - nonlinear_force = np.zeros(2*3*number_of_nodes) - nonlinear_force[::2] = lu_c - nonlinear_force[1::2] = lu_s - - """ - nonlinear_force = np.zeros(2*3*number_of_nodes) - for i in range(3*number_of_nodes): - a = u[2*i] - b = u[2*i+1] - nonlinear_force[2*i:2*i+2] = fem_lu[i,i]*np.array([ - 3/4*(a**3+a*b**2), - 3/4*(a**2*b+b**3) - ]) - """ - - # Add dirichlet boundary conditions - nonlinear_force[self.dirichlet_nodes] = 0 - - return nonlinear_force - - def get_harmonic_balancing_matrices( - self, - fem_m: sparse.lil_array, - fem_c: sparse.lil_array, - fem_k: sparse.lil_array, - angular_frequency: float - ) -> Tuple[sparse.lil_array, sparse.lil_array, sparse.lil_array]: - """Calculates the harmonic balancing matrices from the FEM matrices. - The matrices are getting bigger by a factor of (2*harmonic_order)**2 - - Parameters: - fem_m: FEM mass matrix. - fem_c: FEM damping matrix. - fem_k: FEM stiffness matrix. - - Returns: - Tuple of HB mass matrix, HB damping matrix and HB stiffness - matrix. - """ - harmonic_order = self.harmonic_order - - # Create harmonic balance derivative matrices - nabla_m = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - nabla_c = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - nabla_k = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - for i in range(harmonic_order): - k = i + 1 - nabla_m[2*i:2*i+2, 2*i:2*i+2] = np.array([ - [-(angular_frequency**2*k**2), 0], - [0, -(angular_frequency**2*k**2)] - ]) - nabla_c[2*i:2*i+2, 2*i:2*i+2] = np.array([ - [0, angular_frequency*k], - [-(angular_frequency*k), 0] - ]) - nabla_k[2*i:2*i+2, 2*i:2*i+2] = np.array([ - [1, 0], - [0, 1] - ]) - - # Apply HB derivative matrices on FEM matrices - fem_m = sparse.kron(fem_m, nabla_m, format='lil').tolil() - fem_c = sparse.kron(fem_c, nabla_c, format='lil').tolil() - fem_k = sparse.kron(fem_k, nabla_k, format='lil').tolil() - - return fem_m, fem_c, fem_k - - def get_load_vector( - self, - nodes: npt.NDArray, - values: npt.NDArray, - ) -> npt.NDArray: - """Calculates the load vector (right hand side) vector for the - simulation. - - Parameters: - nodes: Nodes at which the dirichlet value shall be set. - values: Dirichlet value which is set at the corresponding node. - - Returns: - Right hand side vector for the simulation. - """ - number_of_nodes = len(self.mesh_data.nodes) - - # Can be initialized to 0 because external load and volume - # charge density is 0. - f = np.zeros( - 3*number_of_nodes*2*self.harmonic_order, - dtype=np.float64 - ) - - for node, value in zip(nodes, values): - f[node] = value - - return f - - def calculate_tangent_matrix_analytical( - self, - u: npt.NDArray, - hb_m: sparse.lil_array, - hb_c: sparse.lil_array, - hb_k: sparse.lil_array - ) -> sparse.csc_array: - """Calculates the tangent matrix using a priori analytical - calculations. - - Parameters: - u: Solution vector containing mechanical and electrical field - values. - hb_m: Harmonic balancing mass matrix. - hb_c: Harmonic balancing damping matrix. - hb_k: Harmonic balancing stiffness matrix. - - Returns: - The tangent stiffness matrix. - """ - # TODO Only valid for harmonic_order = 1 and nonlinearity_order = 3 - fem_lu = self.lu - lu_diag = fem_lu.diagonal() - number_of_nodes = len(self.mesh_data.nodes) - size = number_of_nodes*3*self.harmonic_order*2 - - center_diag = np.zeros(size) - for i in range(len(lu_diag)): - diag_value_first = 3/4*lu_diag[i]*(3*u[2*i]**2+u[2*i+1]**2) - diag_value_second = 3/4*lu_diag[i]*(u[2*i]**2+3*u[2*i+1]**2) - - center_diag[2*i] = diag_value_first - center_diag[2*i+1] = diag_value_second - - left_diag = np.zeros(size-1) - left_diag[::2] = np.multiply(u[1::2], lu_diag) - - right_diag = np.zeros(size-1) - right_diag[::2] = np.multiply(u[:-1:2], lu_diag) - - return hb_m + hb_c + hb_k + sparse.diags( - diagonals=[ - left_diag, - center_diag, - right_diag, - ], - offsets=[-1, 0, 1], - format="csc" - ) - - - # -------- Static functions -------- - - @staticmethod - def calculate_tangent_matrix_numerical( - u: npt.NDArray, - nonlinear_force: npt.NDArray, - residual: Callable - ) -> sparse.csc_array: - """Calculates the tanget stiffness matrix numerically using a first - order finite difference approximation. - Important note: Calculation is very slow. - - Parameters: - u: Current solution vector. - nonlinear_force: Vector of nonlinear forces. - residual: Residual function. - - Returns: - Tangential stiffness matrix in sparse csc format. - """ - h = 1e-14 - k_tangent = sparse.lil_matrix( - (len(u), len(u)), - dtype=np.float64 - ) - for i in range(len(u)): - e_m = np.zeros(len(u)) - e_m[i] = h - k_tangent[:, i] = ( - residual(u+e_m, nonlinear_force) - - residual(u, nonlinear_force) - )/h - - return k_tangent.tocsc() - - - @staticmethod - def apply_dirichlet_bc( - m: sparse.lil_array, - c: sparse.lil_array, - k: sparse.lil_array, - nodes: npt.NDArray - ) -> Tuple[ - sparse.csc_array, - sparse.csc_array, - sparse.csc_array - ]: - """Applies dirichlet boundary conditions to the given matrix based on - the given nodes. Essentially the rows of the m c and k matrices at the - node indices are set to 0. The element of the k matrix at position - (node_index, node_index) is set to 1. - - Parameters: - m: Mass matrix. - c: Damping matrix. - k: Stiffness matrix. - nodes: Nodes at which the dirichlet boundary conditions shall be - applied. - - Returns: - Modified mass, damping and stiffness matrix. - """ - # Set rows of matrices to 0 and diagonal of K to 1 (at node points) - # Matrices for u_r component - for node in nodes: - # Set rows to 0 - if m is not None: - m[node, :] = 0 - if c is not None: - c[node, :] = 0 - if k is not None: - k[node, :] = 0 - - # Set diagonal values to 1 - k[node, node] = 1 - - return m.tocsc(), c.tocsc(), k.tocsc() diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py deleted file mode 100644 index a4703fc..0000000 --- a/plutho/simulation/nonlinear/piezo_time.py +++ /dev/null @@ -1,410 +0,0 @@ -"""Module for the simulation of nonlinear piezoelectric systems""" - -# Third party libraries -from typing import Tuple, Union -import numpy as np -import numpy.typing as npt -import scipy -from scipy import sparse -from scipy.sparse import linalg - -# Local libraries -from .base import assemble, Nonlinearity -from ..base import SimulationData, MeshData, MaterialData, FieldType, \ - create_node_points -from plutho.simulation.piezo_time import calculate_charge -from plutho.materials import MaterialManager -from plutho.mesh.mesh import Mesh - - -class NLPiezoTime: - - def __init__( - self, - simulation_name: str, - mesh: Mesh, - nonlinearity: Nonlinearity - ): - self.simulation_name = simulation_name - - # Setup simulation data - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - # Setup mesh and material data - nodes, elements = mesh.get_mesh_nodes_and_elements() - self.mesh = mesh - self.mesh_data = MeshData(nodes, elements, mesh.element_order) - self.node_points = create_node_points( - nodes, - elements, - mesh.element_order - ) - - self.material_manager = MaterialManager(len(elements)) - self.nonlinearity = nonlinearity - self.nonlinearity.set_mesh_data(self.mesh_data, self.node_points) - - def add_material( - self, - material_name: str, - material_data: MaterialData, - physical_group_name: str - ): - """Adds a material to the simulation - - Parameters: - material_name: Name of the material. - material_data: Material data. - physical_group_name: Name of the physical group for which the - material shall be set. If this is None or "" the material will - be set for the whole model. - """ - if physical_group_name is None or physical_group_name == "": - element_indices = np.arange(len(self.mesh_data.elements)) - else: - element_indices = self.mesh.get_elements_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - self.material_manager.add_material( - material_name, - material_data, - physical_group_name, - element_indices - ) - - def add_dirichlet_bc( - self, - field_type: FieldType, - physical_group_name: str, - values: npt.NDArray - ): - """Adds a dirichlet boundary condition (dirichlet bc) to the - simulation. The given values are for each time step. Therefore each - value is applied to every node from the given physical_group. - - Parameters: - field_type: The type of field for which the bc shall be set. - physical_group_name: Name of the physical group for which the - boundary condition shall be set. - values: List of values per time step. Value of the bc for each time - step. The value for one time step is applied to every element - equally. Length: number_of_time_step. - """ - # Save boundary condition for serialization - self.boundary_conditions.append({ - "field_type": field_type.name, - "physical_group": physical_group_name, - "values": values.tolist() - }) - - # Apply the given values for all nodes from the given physical_group - node_indices = self.mesh.get_nodes_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - number_of_nodes = len(self.mesh_data.nodes) - - for node_index in node_indices: - real_index = 0 - - # Depending on the variable type and the simulation types - # the corresponding field variable may be found at different - # positions of the solution vector - if field_type is FieldType.PHI: - real_index = 2*number_of_nodes+node_index - elif field_type is FieldType.U_R: - real_index = 2*node_index - elif field_type is FieldType.U_Z: - real_index = 2*node_index+1 - else: - raise ValueError(f"Unknown variable type {field_type}") - - self.dirichlet_nodes.append(real_index) - self.dirichlet_values.append(values) - - def clear_dirichlet_bcs(self): - """Resets the dirichlet boundary conditions.""" - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - def set_electrode(self, electrode_name: str): - self.electrode_elements = self.mesh.get_elements_by_physical_groups( - [electrode_name] - )[electrode_name] - self.electrode_normals = np.tile( - [0, 1], - (len(self.electrode_elements), 1) - ) - - def assemble(self): - """Redirect to general nonlinear assembly function""" - self.material_manager.initialize_materials() - m, c, k = assemble( - self.mesh_data, - self.material_manager - ) - self.m = m - self.c = c - self.k = k - - self.nonlinearity.assemble(m, c, k) - - def simulate( - self, - delta_t: float, - number_of_time_steps: int, - gamma: float, - beta: float, - tolerance: float = 1e-11, - max_iter: int = 300, - u_start: Union[npt.NDArray, None] = None - ): - dirichlet_nodes = np.array(self.dirichlet_nodes) - dirichlet_values = np.array(self.dirichlet_values) - - self.simulation_data = SimulationData( - delta_t, - number_of_time_steps, - gamma, - beta - ) - number_of_nodes = len(self.mesh_data.nodes) - - m = self.m.copy() - c = self.c.copy() - k = self.k.copy() - - # Time integration constants - a1 = 1/(beta*(delta_t**2)) - a2 = 1/(beta*delta_t) - a3 = (1-2*beta)/(2*beta) - a4 = gamma/(beta*delta_t) - a5 = 1-gamma/beta - a6 = (1-gamma/(2*beta))*delta_t - - # Init arrays - # Displacement u which is calculated - u = np.zeros( - (3*number_of_nodes, number_of_time_steps), - dtype=np.float64 - ) - # Displacement u derived after t (du/dt) - v = np.zeros( - (3*number_of_nodes, number_of_time_steps), - dtype=np.float64 - ) - # v derived after u (d^2u/dt^2) - a = np.zeros( - (3*number_of_nodes, number_of_time_steps), - dtype=np.float64 - ) - # Charge - q = np.zeros(number_of_time_steps, dtype=np.float64) - - # Apply dirichlet bc to matrices - m, c, k = self._apply_dirichlet_bc( - m, - c, - k, - dirichlet_nodes - ) - self.nonlinearity.apply_dirichlet_bc(dirichlet_nodes) - - # Residual for the newton iterations - def residual(next_u, current_u, v, a, f): - return ( - m@(a1*(next_u-current_u)-a2*v-a3*a) - + c@(a4*(next_u-current_u)+a5*v+a6*a) - + k@next_u+self.nonlinearity.evaluate_force_vector( - next_u, - m, - c, - k - )-f - ) - - if u_start is not None: - u[:, 0] = u_start - - print("Starting nonlinear time domain simulation") - for time_index in range(number_of_time_steps-1): - # Calculate load vector - f = self._get_load_vector( - dirichlet_nodes, - dirichlet_values[:, time_index+1] - ) - - # Values of the current time step - current_u = u[:, time_index] - current_v = v[:, time_index] - current_a = a[:, time_index] - - # Current iteration value of u of the next time step - # as a start value this is set to the last converged u of the last - # time step - u_i = current_u.copy() - - # Best values during newton are saved - best_u_i = u_i.copy() - best_norm = scipy.linalg.norm(residual( - u_i, - current_u, - current_v, - current_a, - f - )) - - # Placeholder for result of newton - next_u = np.zeros(3*number_of_nodes) - self.converged = False - - if best_norm > tolerance: - # Start newton iterations - for i in range(max_iter): - # Calculate tangential stiffness matrix - tangent_matrix = self._calculate_tangent_matrix( - u_i, - m, - c, - k - ) - delta_u = linalg.spsolve( - (a1*m+a4*c+tangent_matrix), - residual(u_i, current_u, current_v, current_a, f) - ) - u_i_next = u_i - delta_u - - # Check for convergence - norm = scipy.linalg.norm(residual( - u_i_next, current_u, current_v, current_a, f - )) - - if norm < tolerance: - print( - f"Newton converged at time step {time_index} " - f"after {i+1} iteration(s)" - ) - # print(u_i_next) - next_u = u_i_next - self.converged = True - break - elif norm < best_norm: - best_norm = norm - best_u_i = u_i_next - - if i % 100 == 0 and i > 0: - print("Iteration:", i) - - # Update for next iteration - u_i = u_i_next - if not self.converged: - print( - "Newton did not converge.. Choosing best value: " - f"{best_norm}" - ) - next_u = best_u_i - else: - print("Start value norm already below tolerance") - next_u = best_u_i - - # Calculate next v and a - a[:, time_index+1] = ( - a1*(next_u-current_u) - - a2*current_v - - a3*current_a - ) - v[:, time_index+1] = ( - a4*(next_u-current_u) - + a5*current_v - + a6*current_a - ) - - # Set u array value - u[:, time_index+1] = next_u - - # Calculate charge - if (self.electrode_elements is not None - and self.electrode_elements.shape[0] > 0): - q[time_index] = calculate_charge( - u[:, time_index+1], - self.material_manager, - self.electrode_elements, - self.electrode_normals, - self.mesh_data.nodes, - self.mesh_data.element_order - ) - - self.u = u - self.q = q - - def _get_load_vector( - self, - nodes: npt.NDArray, - values: npt.NDArray - ) -> npt.NDArray: - """Calculates the load vector (right hand side) vector for the - simulation. - - Parameters: - nodes: Nodes at which the dirichlet value shall be set. - values: Dirichlet value which is set at the corresponding node. - - Returns: - Right hand side vector for the simulation. - """ - number_of_nodes = len(self.mesh_data.nodes) - - # Can be initialized to 0 because external load and volume - # charge density is 0. - f = np.zeros(3*number_of_nodes, dtype=np.float64) - - for node, value in zip(nodes, values): - f[node] = value - - return f - - def _apply_dirichlet_bc( - self, - m: sparse.lil_array, - c: sparse.lil_array, - k: sparse.lil_array, - nodes: npt.NDArray - ) -> Tuple[ - sparse.csc_array, - sparse.csc_array, - sparse.csc_array - ]: - # Set rows of matrices to 0 and diagonal of K to 1 (at node points) - m[nodes, :] = 0 - c[nodes, :] = 0 - k[nodes, :] = 0 - k[nodes, nodes] = 1 - - return m.tocsc(), c.tocsc(), k.tocsc() - - def _calculate_tangent_matrix( - self, - u: npt.NDArray, - m: sparse.csc_array, - c: sparse.csc_array, - k: sparse.csc_array - ): - # TODO Duplicate function in piezo_stationary.py - """Calculates the tangent matrix based on an analytically - expression. - - Parameters: - u: Current mechanical displacement. - k: FEM stiffness matrix. - ln: FEM nonlinear stiffness matrix. - """ - linear = k - nonlinear = self.nonlinearity.evaluate_jacobian( - u, m, c, k - ) - - return (linear + nonlinear).tocsc() diff --git a/plutho/simulation/piezo_time.py b/plutho/simulation/piezo_time.py deleted file mode 100644 index b0f834a..0000000 --- a/plutho/simulation/piezo_time.py +++ /dev/null @@ -1,479 +0,0 @@ -"""Main module for the simulation of piezoelectric systems.""" - -# Python standard libraries -from typing import List -import numpy as np -import numpy.typing as npt -import scipy.sparse as sparse -import scipy.sparse.linalg as slin - -# Local libraries -from .base import SimulationData, MeshData, \ - gradient_local_shape_functions_2d, \ - local_to_global_coordinates, b_operator_global, integral_m, \ - integral_ku, integral_kuv, integral_kve, apply_dirichlet_bc, \ - line_quadrature, create_node_points -from ..materials import MaterialManager - - -def charge_integral_u( - node_points: npt.NDArray, - u_e: npt.NDArray, - piezo_matrix: npt.NDArray, - element_order: int -): - """Calculates the integral of eBu of the given element. - - Parameters: - node_points: List of node points [[x1, x2], [y1, y2]] of - one line. - u_e: List of u values at the nodes of the triangle - [u1_r, u1_z, u2_r, u2_z]. - piezo_matrix: Piezo matrix for the current element (e matrix). - element_order: Order of the shape functions. - - Returns: - Float: Integral of eBu of the current triangle. - """ - def inner(s): - dn = gradient_local_shape_functions_2d(s, 0, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - - b_opt_global = b_operator_global( - node_points, - jacobian_inverted_t, - s, - 0, - element_order - ) - r = local_to_global_coordinates(node_points, s, 0, element_order)[0] - - return -np.dot(np.dot(piezo_matrix, b_opt_global), u_e)*r - - return line_quadrature(inner, element_order) - - -def charge_integral_v( - node_points: npt.NDArray, - v_e: npt.NDArray, - permittivity_matrix: npt.NDArray, - element_order: int -): - """Calculates the integral of epsilonBVe of the given element. - - Parameters: - node_points: List of node points [[x1, x2], [y1, y2]] of - one line. - v_e: List of u values at the nodes of the line - [v1, v2]. - permittivity_matrix: Permittivity matrix for the current - element (e matrix). - element_order: Order of the shape functions. - - Returns: - Float: Integral of epsilonBVe of the current triangle. - """ - def inner(s): - dn = gradient_local_shape_functions_2d(s, 0, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - - global_dn = np.dot(jacobian_inverted_t, dn) - r = local_to_global_coordinates(node_points, s, 0, element_order)[0] - - return -np.dot(np.dot(permittivity_matrix, global_dn), v_e)*r - - return line_quadrature(inner, element_order) - - -def calculate_charge( - u: npt.NDArray, - material_manager: MaterialManager, - elements: npt.NDArray, - element_normals: npt.NDArray, - nodes: npt.NDArray, - element_order: int, - is_complex: bool = False -) -> float: - """Calculates the charge of the given elements. - - Parameters: - u: List of u values for every node. - material_manager: MaterialManager object. - elements: Elements for which the charge shall be calculated. - element_normals: List of element normal vectors (index corresponding to - the elements list). - nodes: All nodes used in the simulation. - element_order: Order of the shape functions. - """ - points_per_element = int(1/2*(element_order+1)*(element_order+2)) - number_of_nodes = len(nodes) - q = 0 - - for element_index, element in enumerate(elements): - node_points = np.zeros(shape=(2, points_per_element)) - for node_index in range(points_per_element): - node_points[:, node_index] = [ - nodes[element[node_index]][0], - nodes[element[node_index]][1] - ] - - # The integration is always done along the first two points of the - # triangle. - # TODO Not true for element_order > 1 - jacobian_det = np.sqrt( - np.square(nodes[element[1]][0]-nodes[element[0]][0]) - + np.square(nodes[element[1]][1]-nodes[element[0]][1]) - ) - - # Altough it is a line integral, it is necessary to take the field of - # all element points even if they are not on the line. - # This is due to the charge calculation needing the derivatives of the - # shape function (dN_i/dr and dN_i/dz) which are (+/-) 1 along the - # whole triangle for all 3 points. - # TODO Is this true? - if is_complex: - u_e = np.zeros(2*points_per_element, dtype=np.complex128) - ve_e = np.zeros(points_per_element, dtype=np.complex128) - else: - u_e = np.zeros(2*points_per_element) - ve_e = np.zeros(points_per_element) - - for i in range(points_per_element): - u_e[2*i] = u[2*element[i]] - u_e[2*i+1] = u[2*element[i]+1] - ve_e[i] = u[element[i]+2*number_of_nodes] - - q_u = charge_integral_u( - node_points, - u_e, - material_manager.get_piezo_matrix(element_index), - element_order - ) * 2 * np.pi * jacobian_det - q_v = charge_integral_v( - node_points, - ve_e, - material_manager.get_permittivity_matrix(element_index), - element_order - ) * 2 * np.pi * jacobian_det - - # Now take the component normal to the line (outer direction) - q_u_normal = np.dot(q_u, element_normals[element_index]) - q_v_normal = np.dot(q_v, element_normals[element_index]) - - q += q_u_normal - q_v_normal - - return q - - -class PiezoSimTime: - """Class for the simulation of mechanical-electric fields. - - Parameters: - mesh_data: MeshData format. - material_manager: MaterialManager object. - simulation_data: SimulationData format. - - Attributes: - mesh_data: Contains the mesh information. - material_manager: Contains information about the materials. - simulation_data: Contains the information about the simulation. - dirichlet_nodes: List of nodes for which dirichlet values are set. - dirichlet_values: List of values of the corresponding dirichlet nodes - per time step. - m: Mass matrix. - c: Damping matrix. - k: Stiffness matrix. - u: Resulting vector of size 4*number_of_nodes containing u_r, u_z - and v. u[node_index, time_index]. An offset needs to be - added to the node_index in order to access the different field - paramteters: - u_r: 0, - u_z: 1*number_of_nodes, - v: 2*number_of_nodes - q: Resulting charges for each time step q[time_index]. - """ - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - simulation_data: SimulationData - - # Dirichlet boundary condition - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray - - # FEM matrices - m: sparse.lil_array - c: sparse.lil_array - k: sparse.lil_array - - # Resulting fields - u: npt.NDArray - q: npt.NDArray - - # Internal simulation data - node_points: npt.NDArray - - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - simulation_data: SimulationData - ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.simulation_data = simulation_data - - def assemble(self): - """Assembles the FEM matrices based on the set mesh_data and - material_data. - The matrices are stored in self.m, self.c, self.k. - """ - # TODO Assembly takes to long rework this algorithm - # Maybe the 2x2 matrix slicing is not very fast - nodes = self.mesh_data.nodes - elements = self.mesh_data.elements - element_order = self.mesh_data.element_order - self.node_points = create_node_points(nodes, elements, element_order) - - number_of_nodes = len(nodes) - mu = sparse.lil_matrix( - (2*number_of_nodes, 2*number_of_nodes), - dtype=np.float64 - ) - ku = sparse.lil_matrix( - (2*number_of_nodes, 2*number_of_nodes), - dtype=np.float64 - ) - kuv = sparse.lil_matrix( - (2*number_of_nodes, number_of_nodes), - dtype=np.float64 - ) - kve = sparse.lil_matrix( - (number_of_nodes, number_of_nodes), - dtype=np.float64 - ) - - for element_index, element in enumerate(self.mesh_data.elements): - node_points = self.node_points[element_index] - - # TODO Check if its necessary to calculate all integrals - # --> Dirichlet nodes could be leaved out? - # Multiply with jac_det because its integrated with respect to - # local coordinates. - # Mutiply with 2*pi because theta is integrated from 0 to 2*pi - mu_e = ( - self.material_manager.get_density(element_index) - * integral_m(node_points, element_order) - * 2 * np.pi - ) - ku_e = ( - integral_ku( - node_points, - self.material_manager.get_elasticity_matrix(element_index), - element_order - ) - * 2 * np.pi - ) - kuv_e = ( - integral_kuv( - node_points, - self.material_manager.get_piezo_matrix(element_index), - element_order - ) - * 2 * np.pi - ) - kve_e = ( - integral_kve( - node_points, - self.material_manager.get_permittivity_matrix( - element_index - ), - element_order - ) - * 2 * np.pi - ) - - # Now assemble all element matrices - for local_p, global_p in enumerate(element): - for local_q, global_q in enumerate(element): - # Mu_e is returned as a 3x3 matrix (should be 6x6) - # because the values are the same. - # Only the diagonal elements have values. - mu[2*global_p, 2*global_q] += mu_e[local_p][local_q] - mu[2*global_p+1, 2*global_q+1] += mu_e[local_p][local_q] - - # Ku_e is a 6x6 matrix and 2x2 matrices are sliced out - # of it. - ku[2*global_p:2*global_p+2, 2*global_q:2*global_q+2] += \ - ku_e[2*local_p:2*local_p+2, 2*local_q:2*local_q+2] - - # KuV_e is a 6x3 matrix and 2x1 vectors are sliced out. - # [:, None] converts to a column vector. - kuv[2*global_p:2*global_p+2, global_q] += \ - kuv_e[2*local_p:2*local_p+2, local_q][:, None] - - # KVe_e is a 3x3 matrix and therefore the values can - # be set directly. - kve[global_p, global_q] += kve_e[local_p, local_q] - - # Calculate damping matrix - # Currently in the simulations only one material is used - # TODO Add algorithm to calculate cu for every element on its own - # in order to use multiple materials - cu = ( - self.material_manager.get_alpha_m(0) * mu - + self.material_manager.get_alpha_k(0) * ku - ) - - # Calculate block matrices - zeros1x1 = np.zeros((number_of_nodes, number_of_nodes)) - zeros2x1 = np.zeros((2*number_of_nodes, number_of_nodes)) - zeros1x2 = np.zeros((number_of_nodes, 2*number_of_nodes)) - - m = sparse.bmat([ - [mu, zeros2x1], - [zeros1x2, zeros1x1], - ]) - c = sparse.bmat([ - [cu, zeros2x1], - [zeros1x2, zeros1x1], - ]) - k = sparse.bmat([ - [ku, kuv], - [kuv.T, -1*kve] - ]) - - self.m = m.tolil() - self.c = c.tolil() - self.k = k.tolil() - - def solve_time( - self, - electrode_elements: npt.NDArray, - electrode_normals: npt.NDArray - ): - """Runs the simulation using the assembled m, c and k matrices as well - as the set excitation. - Calculates the displacement field and potential field of the given - model (all saved in u). - Calculates the electric charge in the electrode. - Saves u[node_index, time_index] and q[time_index]. - - Parameters: - electrode_elements: List of all elements which are in - the electrode. - electrode_normals: List of normal vectors corresponding to the - electrode_elements list. - """ - number_of_time_steps = self.simulation_data.number_of_time_steps - beta = self.simulation_data.beta - gamma = self.simulation_data.gamma - delta_t = self.simulation_data.delta_t - - m = self.m - c = self.c - k = self.k - - if m is None or c is None or k is None: - print("Cannot solve simulation since matrices are not assembled") - - # Init arrays - # Displacement u which is calculated - u = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) - # Displacement u derived after t (du/dt) - v = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) - # v derived after u (d^2u/dt^2) - a = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) - - # Charge calculated during simulation (for impedance) - q = np.zeros(number_of_time_steps, dtype=np.float64) - - m, c, k = apply_dirichlet_bc( - m, - c, - k, - self.dirichlet_nodes - ) - - k_star = (k+gamma/(beta*delta_t)*c+1/(beta*delta_t**2)*m).tocsr() - - # Set initial value of potential at electrode nodes to given excitation - # TODO Could be removed? - for index, node in enumerate(self.dirichlet_nodes): - u[node, 0] = self.dirichlet_values[index, 0] - - print("Starting simulation") - for time_index in range(number_of_time_steps-1): - # Calculate load vector and add dirichlet boundary conditions - f = self.get_load_vector( - self.dirichlet_nodes, - self.dirichlet_values[:, time_index+1] - ) - - # Perform Newmark method - # Predictor step - u_tilde = ( - u[:, time_index] - + delta_t*v[:, time_index] - + delta_t**2/2*(1-2*beta)*a[:, time_index] - ) - v_tilde = ( - v[:, time_index] - + (1-gamma)*delta_t*a[:, time_index] - ) - - # Solve for next time step - u[:, time_index+1] = slin.spsolve( - k_star, ( - f - - c*v_tilde - + (1/(beta*delta_t**2)*m - + gamma/(beta*delta_t)*c)*u_tilde)) - # Perform corrector step - a[:, time_index+1] = (u[:, time_index+1]-u_tilde)/(beta*delta_t**2) - v[:, time_index+1] = v_tilde + gamma*delta_t*a[:, time_index+1] - - # Calculate charge - if electrode_elements is not None: - q[time_index+1] = calculate_charge( - u[:, time_index+1], - self.material_manager, - electrode_elements, - electrode_normals, - self.mesh_data.nodes, - self.mesh_data.element_order - ) - - if (time_index + 1) % 100 == 0: - print(f"Finished time step {time_index+1}") - - self.q = q - self.u = u - - def get_load_vector( - self, - nodes: npt.NDArray, - values: npt.NDArray - ) -> npt.NDArray: - """Calculates the load vector (right hand side) vector for the - simulation. - - Parameters: - nodes: Nodes at which the dirichlet value shall be set. - values: Dirichlet value which is set at the corresponding node. - - Returns: - Right hand side vector for the simulation. - """ - number_of_nodes = len(self.mesh_data.nodes) - - # Can be initialized to 0 because external load and volume - # charge density is 0. - f = np.zeros(3*number_of_nodes, dtype=np.float64) - - for node, value in zip(nodes, values): - f[node] = value - - return f diff --git a/plutho/simulations/__init__.py b/plutho/simulations/__init__.py new file mode 100644 index 0000000..1990f22 --- /dev/null +++ b/plutho/simulations/__init__.py @@ -0,0 +1,5 @@ +from .piezo_freq import * +from .piezo_time import * +from .thermo_piezo_time import * +from .thermo_time import * +from .helpers import * diff --git a/plutho/simulations/helpers.py b/plutho/simulations/helpers.py new file mode 100644 index 0000000..45b532e --- /dev/null +++ b/plutho/simulations/helpers.py @@ -0,0 +1,224 @@ +"""This module contains various helper functions needed in some simulations. +""" + +# Python standard libraries +from typing import Union, Tuple + +# Third party libraries +import numpy as np +import numpy.typing as npt +from scipy import sparse + +# Local libraries +from .integrals import integral_charge_u, integral_charge_v, integral_volume +from ..materials import MaterialManager + + +__all__ = [ + "create_node_points", + "calculate_volumes", + "get_avg_temp_field_per_element", + "calculate_charge", +] + + +def create_node_points( + nodes: npt.NDArray, + elements: npt.NDArray, + element_order: int +) -> npt.NDArray: + """Create the local node data for every given element. + + Parameters: + nodes: Nodes of the mesh. + elements: Elements of the mesh. + element_order: Order of the elements. + + Returns: + List of nodes. + """ + points_per_element = int(1/2*(element_order+1)*(element_order+2)) + + node_points = np.zeros( + shape=(len(elements), 2, points_per_element) + ) + + for element_index, element in enumerate(elements): + # Get node points of element in format + # [x1 x2 x3 ... xn] + # [y1 y2 y3 ... yn] where (xi, yi) are the coordinates for Node i + for node_index in range(points_per_element): + node_points[element_index, :, node_index] = [ + nodes[element[node_index]][0], + nodes[element[node_index]][1] + ] + + return node_points + + +def calculate_volumes( + node_points: npt.NDArray, + element_order: int +) -> npt.NDArray: + """Calculates the volume of each element. The element information + is given by the local_element_data + + Parameters: + node_points: Node points of all elements for which the volume shall be + calculated. + element_order: Order of the elements. + + Returns: + List of volumes of the elements. + """ + volumes = [] + + number_of_elements = node_points.shape[0] + + for element_index in range(number_of_elements): + volumes.append( + integral_volume(node_points[element_index], element_order) + *2*np.pi + ) + + return np.array(volumes) + + +def get_avg_temp_field_per_element( + theta: npt.NDArray, + elements: npt.NDArray +) -> npt.NDArray: + """Returns the average temperature for each element. + + Parameters: + theta: Temperatures for each node. + elements: List of elements. + + Returns: + Mean temperature for each element. + """ + theta_elements = np.zeros(len(elements)) + + for element_index, element in enumerate(elements): + theta_elements[element_index] = np.mean(theta[element]) + + return np.array(theta_elements) + + +def calculate_charge( + u: npt.NDArray, + material_manager: MaterialManager, + elements: npt.NDArray, + element_normals: npt.NDArray, + nodes: npt.NDArray, + element_order: int, + is_complex: bool = False +) -> float: + """Calculates the charge of the given elements. + + Parameters: + u: List of u values for every node. + material_manager: MaterialManager object. + elements: Elements for which the charge shall be calculated. + element_normals: List of element normal vectors (index corresponding to + the elements list). + nodes: All nodes used in the simulation. + element_order: Order of the shape functions. + is_complex: Set to true of u can be complex typed. + """ + points_per_element = int(1/2*(element_order+1)*(element_order+2)) + number_of_nodes = len(nodes) + q = 0 + + for element_index, element in enumerate(elements): + node_points = np.zeros(shape=(2, points_per_element)) + for node_index in range(points_per_element): + node_points[:, node_index] = [ + nodes[element[node_index]][0], + nodes[element[node_index]][1] + ] + + # The integration is always done along the first two points of the + # triangle. + # TODO Not true for element_order > 1 + jacobian_det = np.sqrt( + np.square(nodes[element[1]][0]-nodes[element[0]][0]) + + np.square(nodes[element[1]][1]-nodes[element[0]][1]) + ) + + # Altough it is a line integral, it is necessary to take the field of + # all element points even if they are not on the line. + # This is due to the charge calculation needing the derivatives of the + # shape function (dN_i/dr and dN_i/dz) which are (+/-) 1 along the + # whole triangle for all 3 points. + # TODO Is this true? + if is_complex: + u_e = np.zeros(2*points_per_element, dtype=np.complex128) + ve_e = np.zeros(points_per_element, dtype=np.complex128) + else: + u_e = np.zeros(2*points_per_element) + ve_e = np.zeros(points_per_element) + + for i in range(points_per_element): + u_e[2*i] = u[2*element[i]] + u_e[2*i+1] = u[2*element[i]+1] + ve_e[i] = u[element[i]+2*number_of_nodes] + + q_u = integral_charge_u( + node_points, + u_e, + material_manager.get_piezo_matrix(element_index), + element_order + ) * 2 * np.pi * jacobian_det + q_v = integral_charge_v( + node_points, + ve_e, + material_manager.get_permittivity_matrix(element_index), + element_order + ) * 2 * np.pi * jacobian_det + + # Now take the component normal to the line (outer direction) + q_u_normal = np.dot(q_u, element_normals[element_index]) + q_v_normal = np.dot(q_v, element_normals[element_index]) + + q += q_u_normal - q_v_normal + + return q + +def mat_apply_dbcs( + m: Union[sparse.lil_array, None], + c: Union[sparse.lil_array, None], + k: sparse.lil_array, + nodes: npt.NDArray +) -> Tuple[ + Union[sparse.lil_array, None], + Union[sparse.lil_array, None], + sparse.lil_array +]: + """Prepares the given matrices m, c and k for the dirichlet boundary + conditions. This is done by setting the corresponding rows to 0 + excepct for the node which will contain the specific value (this is set + to 1). Right now only boundary conditions for v and u_r can be set, not + for u_z (not neede yet). + + Parameters: + m: Mass matrix M. + c: Damping matrix C. + k: Stiffness matrix K. + nodes: List of nodes at which a dirichlet boundary condition + shall be applied. + + Returns: + Modified mass, damping and stiffness matrix. + """ + # Set rows of matrices to 0 and diagonal of K to 1 (at node points) + if m is not None: + m[nodes, :] = 0 + + if c is not None: + c[nodes, :] = 0 + + k[nodes, :] = 0 + k[nodes, nodes] = 1 + + return m, c, k diff --git a/plutho/simulation/base.py b/plutho/simulations/integrals.py similarity index 56% rename from plutho/simulation/base.py rename to plutho/simulations/integrals.py index 7b4d796..4ad5bb9 100644 --- a/plutho/simulation/base.py +++ b/plutho/simulations/integrals.py @@ -1,159 +1,22 @@ -"""Module for base functionalities needed for the simulations.""" +"""General module for the implementation of various functions on local and +global shape functions as well as integration techniques""" # Python standard libraries -import os -import json -from typing import Tuple, Callable, List, Any, Union -from dataclasses import dataclass, fields -from enum import Enum +from typing import Callable + +# Third party libraries import numpy as np import numpy.typing as npt -from scipy import sparse - -# Local libraries -from ..mesh import Mesh - -# -------- ENUMS AND DATACLASSES -------- - -class FieldType(Enum): - """Possible field types which are calculated using differnet simulations. - """ - U_R = 0 - U_Z = 1 - PHI = 2 - THETA = 3 - - -class ModelType(Enum): - """Containts the model type. Since for different model types - different boundary conditions are needed it is necessary to make a - distinction. - Additionaly the Ring model type needs an appropiate mesh set separately - (using x_offset in the mesh generation). - """ - DISC = "disc" - RING = "ring" - - -class SimulationType(Enum): - """Contains the simulation type. Currently it is possible to have - a simulation with or without thermal field.""" - PIEZOELECTRIC = "piezoelectric" - THERMOPIEZOELECTRIC = "thermo-piezoelectric" - FREQPIEZOELECTRIC = "frequency-piezoelectric" - -@dataclass -class SimulationData: - """Contains data for the simulation itself.""" - delta_t: float - number_of_time_steps: int - gamma: float - beta: float - -@dataclass -class MeshData: - """Contains the mesh data is used in the simulation.""" - nodes: npt.NDArray - elements: npt.NDArray +# +# -------- Local shape functions -------- +# +def local_shape_functions_2d( + s: float, + t: float, element_order: int - - -class ExcitationType(Enum): - """Sets the excitation type of the simulation.""" - SINUSOIDAL = "sinusoidal" - TRIANGULAR_PULSE = "triangular_pulse" - - -@dataclass -class ExcitationInfo: - """Contains information about the excitation. Is used to save the - excitation data in the simulation config file.""" - amplitude: Union[float, npt.NDArray] - frequency: Union[float, npt.NDArray] - excitation_type: ExcitationType - - def asdict(self): - """Returns this object as a dictionary.""" - content = self.__dict__ - if self.frequency is None: - del content["frequency"] - content["excitation_type"] = self.excitation_type.value - return content - - -@dataclass -class MaterialData: - """Contains the plain material data. Some parameters can either be - a float or an array depending if they are temperature dependent. - If they are temperature dependent the index in the array corresponds to - the temperature value from the temperatures array at the same index. - """ - c11: Union[float, npt.NDArray] - c12: Union[float, npt.NDArray] - c13: Union[float, npt.NDArray] - c33: Union[float, npt.NDArray] - c44: Union[float, npt.NDArray] - e15: Union[float, npt.NDArray] - e31: Union[float, npt.NDArray] - e33: Union[float, npt.NDArray] - eps11: Union[float, npt.NDArray] - eps33: Union[float, npt.NDArray] - alpha_m: float - alpha_k: float - thermal_conductivity: float - heat_capacity: float - temperatures: Union[float, npt.NDArray] - density: float - - def to_dict(self): - """Convert the dataclass to dict for json serialization.""" - json_dict = {} - for attribute in fields(self.__class__): - value = getattr(self, attribute.name) - if isinstance(value, float) or isinstance(value, int): - json_dict[attribute.name] = value - elif isinstance(value, np.ndarray): - json_dict[attribute.name] = value.tolist() - else: - raise ValueError( - "Wrong type saved in MaterialData. Value is of type " - f"{type(value)}" - ) - return json_dict - - @staticmethod - def from_dict(contents): - """Convert given dict, e.g. from a json deserialization, to a - MaterialData object.""" - for key, value in contents.items(): - if isinstance(value, List): - contents[key] = np.array(value) - - return MaterialData(**contents) - - @staticmethod - def load_from_file(file_path: str): - """Load the data from given file. - - Parameters: - file_path: Path to the file - """ - if not os.path.exists(file_path): - raise IOError( - "Given file path {} does not exist. Cannot load " - "material data." - ) - - with open(file_path, "r", encoding="UTF-8") as fd: - return MaterialData.from_dict(json.load(fd)) - - -# -------- Local functions and integrals -------- - - -def local_shape_functions_2d(s, t, element_order): +) -> npt.NDArray: """Returns the local linear shape functions based on a reference triangle with corner points [(0,0), (1,0), (1,1)] for the given coordinates. @@ -207,7 +70,11 @@ def local_shape_functions_2d(s, t, element_order): ) -def gradient_local_shape_functions_2d(s, t, element_order) -> npt.NDArray: +def gradient_local_shape_functions_2d( + s: float, + t: float, + element_order: int +) -> npt.NDArray: """Returns the gradient of the local shape functions. Parameters: @@ -285,7 +152,7 @@ def local_to_global_coordinates( s: float, t: float, element_order: int -) -> Any: +) -> npt.NDArray: """Transforms the local coordinates given by dimensions using the node points to the global coordinates r, z. @@ -314,7 +181,7 @@ def b_operator_global( s: float, t: float, element_order: int -): +) -> npt.NDArray: """Calculates the B operator for the local coordinantes which is needed for voigt-notation. The derivates are with respect to the global coordinates r and z. @@ -363,8 +230,10 @@ def b_operator_global( return b - -def integral_m(node_points: npt.NDArray, element_order: int): +# +# -------- Integrals -------- +# +def integral_m(node_points: npt.NDArray, element_order: int) -> npt.NDArray: """Calculates the M integral. Parameters: @@ -375,7 +244,7 @@ def integral_m(node_points: npt.NDArray, element_order: int): Returns: npt.NDArray: 3x3 M matrix for the given element. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_det = np.linalg.det(jacobian) @@ -397,7 +266,7 @@ def integral_ku( node_points: npt.NDArray, elasticity_matrix: npt.NDArray, element_order: int -): +) -> npt.NDArray: """Calculates the Ku integral Parameters: @@ -410,7 +279,7 @@ def integral_ku( Returns: npt.NDArray: 6x6 Ku matrix for the given element. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_inverted_t = np.linalg.inv(jacobian).T @@ -437,7 +306,7 @@ def integral_kuv( node_points: npt.NDArray, piezo_matrix: npt.NDArray, element_order: int -): +) -> npt.NDArray: """Calculates the KuV integral. Parameters: @@ -449,7 +318,7 @@ def integral_kuv( Returns: npt.NDArray: 6x3 KuV matrix for the given element. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_inverted_t = np.linalg.inv(jacobian).T @@ -477,7 +346,7 @@ def integral_kve( node_points: npt.NDArray, permittivity_matrix: npt.NDArray, element_order: int -): +) -> npt.NDArray: """Calculates the KVe integral. Parameters: @@ -490,7 +359,7 @@ def integral_kve( Returns: npt.NDArray: 3x3 KVe matrix for the given element. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_inverted_t = np.linalg.inv(jacobian).T @@ -514,7 +383,7 @@ def energy_integral_theta( node_points: npt.NDArray, theta: npt.NDArray, element_order: int -): +) -> npt.NDArray: """Integrates the given element over the given theta field. Parameters: @@ -524,7 +393,7 @@ def energy_integral_theta( [theta1, theta2, theta3]. element_order: Order of the shape functions. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_det = np.linalg.det(jacobian) @@ -537,7 +406,10 @@ def inner(s, t): return quadratic_quadrature(inner, element_order) -def integral_volume(node_points: npt.NDArray, element_order: int): +def integral_volume( + node_points: npt.NDArray, + element_order: int +) -> npt.NDArray: """Calculates the volume of the triangle given by the node points. HINT: Must be multiplied with 2*np.pi and the jacobian determinant in order to give the correct volume of any rotationsymmetric triangle. @@ -550,7 +422,7 @@ def integral_volume(node_points: npt.NDArray, element_order: int): Returns: Float. Volume of the triangle. """ - def inner(s, t): + def inner(s: float, t: float) -> npt.NDArray: dn = gradient_local_shape_functions_2d(s, t, element_order) jacobian = np.dot(node_points, dn.T) jacobian_det = np.linalg.det(jacobian) @@ -562,7 +434,270 @@ def inner(s, t): return quadratic_quadrature(inner, element_order) -def quadratic_quadrature(func: Callable, element_order: int): +def integral_loss_scs( + node_points: npt.NDArray, + u_e: npt.NDArray, + elasticity_matrix: npt.NDArray, + element_order: int +): + """Calculate sthe integral of dS/dt*c*dS/dt over one triangle in the + frequency domain for the given frequency. + + Parameters: + node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of one + triangle. + u_e: Displacement at this element [u1_r, u1_z, u_2r, u_2z, u_3r, u_3z]. + angular_frequency: Angular frequency at which the loss shall be + calculated. + elasticity_matrix: Elasticity matrix for the current element + (c matrix). + element_order: Order of the shape functions. + """ + def inner(s, t): + dn = gradient_local_shape_functions_2d(s, t, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + jacobian_det = np.linalg.det(jacobian) + + b_opt = b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + s_e = np.dot(b_opt, u_e) + + # return np.dot(dt_s.T, np.dot(elasticity_matrix.T, dt_s))*r + return np.dot( + np.conjugate(s_e).T, + np.dot( + elasticity_matrix.T, + s_e + ) + ) * r * jacobian_det + + return quadratic_quadrature(inner, element_order) + + +def integral_charge_u( + node_points: npt.NDArray, + u_e: npt.NDArray, + piezo_matrix: npt.NDArray, + element_order: int +): + """Calculates the integral of eBu of the given element. + + Parameters: + node_points: List of node points [[x1, x2], [y1, y2]] of + one line. + u_e: List of u values at the nodes of the triangle + [u1_r, u1_z, u2_r, u2_z]. + piezo_matrix: Piezo matrix for the current element (e matrix). + element_order: Order of the shape functions. + + Returns: + Float: Integral of eBu of the current triangle. + """ + def inner(s): + dn = gradient_local_shape_functions_2d(s, 0, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + + b_opt_global = b_operator_global( + node_points, + jacobian_inverted_t, + s, + 0, + element_order + ) + r = local_to_global_coordinates(node_points, s, 0, element_order)[0] + + return -np.dot(np.dot(piezo_matrix, b_opt_global), u_e)*r + + return line_quadrature(inner, element_order) + + +def integral_charge_v( + node_points: npt.NDArray, + v_e: npt.NDArray, + permittivity_matrix: npt.NDArray, + element_order: int +): + """Calculates the integral of epsilonBVe of the given element. + + Parameters: + node_points: List of node points [[x1, x2], [y1, y2]] of + one line. + v_e: List of u values at the nodes of the line + [v1, v2]. + permittivity_matrix: Permittivity matrix for the current + element (e matrix). + element_order: Order of the shape functions. + + Returns: + Float: Integral of epsilonBVe of the current triangle. + """ + def inner(s): + dn = gradient_local_shape_functions_2d(s, 0, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + + global_dn = np.dot(jacobian_inverted_t, dn) + r = local_to_global_coordinates(node_points, s, 0, element_order)[0] + + return -np.dot(np.dot(permittivity_matrix, global_dn), v_e)*r + + return line_quadrature(inner, element_order) + + +def integral_heat_flux( + node_points: npt.NDArray, + heat_flux: npt.NDArray, + element_order: int +) -> npt.NDArray: + """Integrates the heat flux using the shape functions. + + Parameters: + node_points: List of node points [[x1, x2, ..], [y1, y2, ..]]. + heat_flux: Heat flux at the points. + element_order: Order of the shape functions. + + Returns: + npt.NDArray heat flux integral on each point. + """ + def inner(s: float) -> npt.NDArray: + n = local_shape_functions_2d(s, 0, element_order) + r = local_to_global_coordinates(node_points, s, 0, element_order)[0] + + return n*heat_flux*r + + return line_quadrature(inner, element_order) + + +def integral_ktheta( + node_points: npt.NDArray, + element_order: int +) -> npt.NDArray: + """Calculates the Ktheta integral. + + Parameters: + node_points: List of node points [[x1, x2, ..], [y1, y2, ..]]. + element_order: Order of the shape functions. + + Returns: + npt.NDArray: 3x3 Ktheta matrix for the given element. + """ + def inner(s: float, t: float) -> npt.NDArray: + dn = gradient_local_shape_functions_2d(s, t, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + jacobian_det = np.linalg.det(jacobian) + + global_dn = np.dot(jacobian_inverted_t, dn) + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + return np.dot(global_dn.T, global_dn)*r*jacobian_det + + return quadratic_quadrature(inner, element_order) + + +def integral_theta_load( + node_points: npt.NDArray, + mech_loss: float, + element_order: int +) -> npt.NDArray: + """Returns the load value for the temperature field (f) for the specific + element. + + Parameters: + node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of + one triangle. + point_loss: Loss power on each node (heat source). + element_order: Order of the shape functions. + + Returns: + npt.NDArray: f vector value at the specific ndoe + """ + def inner(s: float, t: float) -> npt.NDArray: + dn = gradient_local_shape_functions_2d(s, t, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_det = np.linalg.det(jacobian) + + n = local_shape_functions_2d(s, t, element_order) + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + return n*mech_loss*r*jacobian_det + + return quadratic_quadrature(inner, element_order) + + +def integral_loss_scs_time( + node_points: npt.NDArray, + u_e_t: npt.NDArray, + u_e_t_minus_1: npt.NDArray, + u_e_t_minus_2: npt.NDArray, + delta_t: float, + elasticity_matrix: npt.NDArray, + element_order: int +) -> npt.NDArray: + """Calculates the integral of dS/dt*c*dS/dt over one triangle. Since foward + difference quotient of second oder is used the last 2 values of e_u are + needed. + + Parameters: + node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of + one triangle. + u_e_t: Values of u_e at the current time point. + Format: [u1_r, u1_z, u2_r, u2_z, u3_r, u3_z]. + u_e_t_minus_1: Values of u_e at one time point earlier. Same format + as u_e_t. + u_e_t_minus_2: Values of u_e at two time points earlier. Same format + as u_e_t. + delta_t: Difference between the time steps. + jacobian_inverted_t: Jacobian matrix inverted and transposed, needed + for calculation of global derivatives. + elasticity_matrix: Elasticity matrix for the current element + (c matrix). + element_order: Order of the shape functions. + """ + def inner(s: float, t: float) -> npt.NDArray: + dn = gradient_local_shape_functions_2d(s, t, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + jacobian_det = np.linalg.det(jacobian) + + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + b_opt = b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + + s_e = np.dot(b_opt, u_e_t) + s_e_t_minus_1 = np.dot(b_opt, u_e_t_minus_1) + s_e_t_minus_2 = np.dot(b_opt, u_e_t_minus_2) + dt_s = (3*s_e-4*s_e_t_minus_1+s_e_t_minus_2)/(2*delta_t) + + return np.dot( + dt_s.T, + np.dot( + elasticity_matrix.T, + dt_s + ) + ) * r * jacobian_det + + return quadratic_quadrature(inner, element_order) + + +# +# -------- Numerical calculations -------- +# +def quadratic_quadrature(func: Callable, element_order: int) -> npt.NDArray: """Integrates the given function of 2 variables using gaussian quadrature along 2 variables for a reference triangle. This gives exact results for linear shape functions. @@ -628,7 +763,7 @@ def quadratic_quadrature(func: Callable, element_order: int): return sum -def line_quadrature(func: Callable, element_order: int): +def line_quadrature(func: Callable, element_order: int) -> npt.NDArray: """Integrates the given function of 2 variables along one variable for a reference triangle. This gives exact results for linear shape functions. @@ -669,235 +804,3 @@ def line_quadrature(func: Callable, element_order: int): return sum - -# -------- Boundary condition functions -------- - - -def apply_dirichlet_bc( - m: sparse.lil_array, - c: sparse.lil_array, - k: sparse.lil_array, - nodes: npt.NDArray -) -> Tuple[sparse.lil_array, sparse.lil_array, sparse.lil_array]: - """Prepares the given matrices m, c and k for the dirichlet boundary - conditions. This is done by setting the corresponding rows to 0 - excepct for the node which will contain the specific value (this is set - to 1). Right now only boundary conditions for v and u_r can be set, not - for u_z (not neede yet). - - Parameters: - m: Mass matrix M. - c: Damping matrix C. - k: Stiffness matrix K. - nodes: List of nodes at which a dirichlet boundary condition - shall be applied. - - Returns: - Modified mass, damping and stiffness matrix. - """ - # Set rows of matrices to 0 and diagonal of K to 1 (at node points) - - # Matrices for u_r component - for node in nodes: - # Set rows to 0 - m[node, :] = 0 - c[node, :] = 0 - k[node, :] = 0 - - # Set diagonal values to 1 - k[node, node] = 1 - - return m, c, k - - -def create_dirichlet_bc_nodes_freq( - mesh: Mesh, - amplitudes: npt.NDArray, - number_of_frequencies: int, - set_symmetric_bc: bool = True -) -> Tuple: - """Create the dirichlet boundary condition nodes for a simulation in the - frequency domain. - - Parameters: - mesh: Mesh class for accessing nodes. - amplitudes: List of amplitude per frequency step. - number_of_frequencies: Total number of frequency steps. - set_symmetric_bc: Set to True of the left border of the model is on - the ordinate axis. - - Returns: - The dirichlet nodes and values for u and v respectively. - """ - - # Get nodes from mesh - pg_nodes = mesh.get_nodes_by_physical_groups( - ["Electrode", "Symaxis", "Ground"]) - - electrode_nodes = pg_nodes["Electrode"] - symaxis_nodes = pg_nodes["Symaxis"] - ground_nodes = pg_nodes["Ground"] - - if set_symmetric_bc: - dirichlet_nodes_u = symaxis_nodes - dirichlet_values_u = np.zeros( - (number_of_frequencies, len(dirichlet_nodes_u), 2)) - else: - dirichlet_nodes_u = np.array([]) - dirichlet_values_u = np.array([]) - - # For potential v set electrode to excitation and ground to 0 - dirichlet_nodes_v = np.concatenate((electrode_nodes, ground_nodes)) - dirichlet_values_v = np.zeros( - (number_of_frequencies, len(dirichlet_nodes_v))) - - # Set excitation value for electrode nodes points - for freq_index, amplitude in enumerate(amplitudes): - dirichlet_values_v[freq_index, :len(electrode_nodes)] = \ - amplitude - - return [dirichlet_nodes_u, dirichlet_nodes_v], \ - [dirichlet_values_u, dirichlet_values_v] - - -def create_dirichlet_bc_nodes( - mesh: Mesh, - electrode_excitation: npt.NDArray, - number_of_time_steps: int, - set_symmetric_bc: bool = True -): - """Creates lists of nodes and values used in the simulation to - apply the boundary conditions. - There are 4 lists returning: nodes_u, values_u, nodes_v, values_v. - The values are corresponding to the nodes with the same index. - The lists are set using the given electrode excitation as well - as the information if a symmetric bc is needed - Disc -> Symmetric bc - Ring -> No symmtric bc - The symmetric bc sets the u_r values at the symmetric axis nodes to 0. - The ground electrode nodes are implicitly set to 0. - - Parameters: - electrode_nodes: Nodes in the electrode region. - symaxis_nodes: Nodes on the symmetrical axis (r=0). - ground_nodes: Nodes in the ground region. - electrode_excitation: Excitation values for each time step. - number_of_time_steps: Total number of time steps of the simulation. - - Returns: - Tuple of 2 tuples. The first inner tuple is a tuple containing - the nodes and values for u and the second inner tuple contains - the nodes and values for v. - """ - # Get nodes from gmsh handler - pg_nodes = mesh.get_nodes_by_physical_groups( - ["Electrode", "Symaxis", "Ground"]) - - electrode_nodes = pg_nodes["Electrode"] - symaxis_nodes = pg_nodes["Symaxis"] - ground_nodes = pg_nodes["Ground"] - - # For Potential set "Electrode" to excitation function - # "Ground" is set to 0 - # For displacement set symaxis values to 0. - # Zeros are set for u_r and u_z but the u_z component is not used. - # Others are implicit natrual bcs (Neumann B.C. -> 0). - if set_symmetric_bc: - dirichlet_nodes_u = symaxis_nodes - dirichlet_values_u = np.zeros( - (number_of_time_steps, len(dirichlet_nodes_u), 2)) - else: - dirichlet_nodes_u = np.array([]) - dirichlet_values_u = np.array([]) - - # For potential v set electrode to excitation and ground to 0 - dirichlet_nodes_v = np.concatenate((electrode_nodes, ground_nodes)) - dirichlet_values_v = np.zeros( - (number_of_time_steps, len(dirichlet_nodes_v))) - - # Set excitation value for electrode nodes points - for time_index, excitation_value in enumerate(electrode_excitation): - dirichlet_values_v[time_index, :len(electrode_nodes)] = \ - excitation_value - - return [dirichlet_nodes_u, dirichlet_nodes_v], \ - [dirichlet_values_u, dirichlet_values_v] - - -def create_node_points( - nodes: npt.NDArray, - elements: npt.NDArray, - element_order: int -) -> npt.NDArray: - """Create the local node data for every given element. - - Parameters: - nodes: Nodes of the mesh. - elements: Elements of the mesh. - element_order: Order of the elements. - - Returns: - List of nodes. - """ - points_per_element = int(1/2*(element_order+1)*(element_order+2)) - - node_points = np.zeros( - shape=(len(elements), 2, points_per_element) - ) - - for element_index, element in enumerate(elements): - # Get node points of element in format - # [x1 x2 x3 ... xn] - # [y1 y2 y3 ... yn] where (xi, yi) are the coordinates for Node i - for node_index in range(points_per_element): - node_points[element_index, :, node_index] = [ - nodes[element[node_index]][0], - nodes[element[node_index]][1] - ] - - return node_points - - -def calculate_volumes(node_points: npt.NDArray, element_order): - """Calculates the volume of each element. The element information - is given by the local_element_data - - Parameters: - node_points: Node points of all elements for which the volume shall be - calculated. - - Returns: - List of volumes of the elements. - """ - volumes = [] - - number_of_elements = node_points.shape[0] - - for element_index in range(number_of_elements): - volumes.append( - integral_volume(node_points[element_index], element_order) - *2*np.pi - ) - - return volumes - - -def get_avg_temp_field_per_element( - theta: npt.NDArray, - elements: npt.NDArray -): - """Returns the average temperature for each element. - - Parameters: - theta: Temperatures for each node. - elements: List of elements. - - Returns: - Mean temperature for each element. - """ - theta_elements = np.zeros(len(elements)) - - for element_index, element in enumerate(elements): - theta_elements[element_index] = np.mean(theta[element]) - - return np.array(theta_elements) diff --git a/plutho/simulation/piezo_freq.py b/plutho/simulations/piezo_freq.py similarity index 50% rename from plutho/simulation/piezo_freq.py rename to plutho/simulations/piezo_freq.py index deeacfc..d9df1ff 100644 --- a/plutho/simulation/piezo_freq.py +++ b/plutho/simulations/piezo_freq.py @@ -1,104 +1,39 @@ -"""Main module for the simulation of piezoelectric systems.""" +"""Module for the simulation of freqwuency domain piezoelectric systems.""" # Python standard libraries -from typing import List +from typing import Union import numpy as np import numpy.typing as npt import scipy.sparse as sparse import scipy.sparse.linalg as slin # Local libraries -from .base import MeshData, local_to_global_coordinates, b_operator_global, \ - integral_m, integral_ku, integral_kuv, integral_kve, apply_dirichlet_bc, \ - create_node_points, quadratic_quadrature, \ - calculate_volumes, gradient_local_shape_functions_2d -from .piezo_time import calculate_charge -from ..materials import MaterialManager - - -def loss_integral_scs( - node_points: npt.NDArray, - u_e: npt.NDArray, - elasticity_matrix: npt.NDArray, - element_order: int -): - """Calculate sthe integral of dS/dt*c*dS/dt over one triangle in the - frequency domain for the given frequency. - - Parameters: - node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of one - triangle. - u_e: Displacement at this element [u1_r, u1_z, u_2r, u_2z, u_3r, u_3z]. - angular_frequency: Angular frequency at which the loss shall be - calculated. - elasticity_matrix: Elasticity matrix for the current element - (c matrix). - element_order: Order of the shape functions. - """ - def inner(s, t): - dn = gradient_local_shape_functions_2d(s, t, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - jacobian_det = np.linalg.det(jacobian) - - b_opt = b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] +from .solver import FEMSolver +from .helpers import create_node_points, calculate_volumes, \ + mat_apply_dbcs +from .integrals import integral_m, integral_ku, integral_kuv, integral_kve, \ + integral_loss_scs +from ..enums import SolverType +from ..mesh import Mesh - s_e = np.dot(b_opt, u_e) - - # return np.dot(dt_s.T, np.dot(elasticity_matrix.T, dt_s))*r - return np.dot( - np.conjugate(s_e).T, - np.dot( - elasticity_matrix.T, - s_e - ) - ) * r * jacobian_det - return quadratic_quadrature(inner, element_order) +__all__ = [ + "PiezoFreq" +] -class PiezoSimFreq: - """Class for the simulation of mechanical-electric fields. - - Parameters: - mesh_data: MeshData format. - material_manager: MaterialManager object. - simulation_data: SimulationData format. +class PiezoFreq(FEMSolver): + """Simulation class for frequency-domain piezoelectric simulations. Attributes: - mesh_data: Contains the mesh information. - material_manager: Contains the information about the materials. - simulation_data: Contains the information about the simulation. - dirichlet_nodes: List of nodes for which dirichlet values are set. - dirichlet_values: List of values of the corresponding dirichlet nodes - per time step. - m: Mass matrix. - c: Damping matrix. - k: Stiffness matrix. - u: Resulting vector of size 4*number_of_nodes containing u_r, u_z - and v. u[node_index, time_index]. An offset needs to be - added to the node_index in order to access the different field - paramteters: - u_r: 0, - u_z: 1*number_of_nodes, - v: 2*number_of_nodes - q: Resulting charges for each time step q[time_index]. + node_points: List of node points per elements. + m: Sparse mass matrix. + c: Sparse damping matrix. + k: Sparse stiffness matrix. + mech_loss: Mechanical loss field. """ - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - frequencies: npt.NDArray - - # Dirichlet boundary condition - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray + # Internal simulation data + node_points: npt.NDArray # FEM matrices m: sparse.lil_array @@ -106,22 +41,13 @@ class PiezoSimFreq: k: sparse.lil_array # Resulting fields - u: npt.NDArray - q: npt.NDArray - mech_loss: npt.NDArray + mech_loss: Union[npt.NDArray, None] - # Internal simulation data - node_points: npt.NDArray + def __init__(self, simulation_name: str, mesh: Mesh): + super().__init__(simulation_name, mesh) - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - frequencies: npt.NDArray - ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.frequencies = frequencies + self.solver_type = SolverType.PiezoFreq + self.mech_loss = None def assemble(self): """Assembles the FEM matrices based on the set mesh_data and @@ -130,6 +56,12 @@ def assemble(self): """ # TODO Assembly takes to long rework this algorithm # Maybe the 2x2 matrix slicing is not very fast + self.material_manager.initialize_materials() + + # Check if some materials are set + if len(self.material_manager.materials) == 0: + raise ValueError("Before assembly some materials must be added.") + nodes = self.mesh_data.nodes elements = self.mesh_data.elements element_order = self.mesh_data.element_order @@ -244,48 +176,44 @@ def assemble(self): self.c = c.tolil() self.k = k.tolil() - def solve_frequency( + def simulate( self, - electrode_elements: npt.NDArray, - electrode_normals: npt.NDArray, - calculate_mech_loss: bool + frequencies: npt.NDArray, + calculate_mech_loss: bool = False ): - """Run the frequency simulation using the given frequencies. - If electrode elements are given the charge at those elements is - calculated. + """Run the frequency domain simulation for each given frequency. + The resulting displacement field is saved in self.u. Parameters: frequencies: Array of frequencies at which the simulation is done. - electrode_elements: Array of element indices. At those indices - the charge is calculated, summed up and saved in q. - electrode_normals: List of normal vectors corresponding to the - electrode_elements list. + calculate_mech_loss: Set to true if the mechanical losses shall be + calculated. They are saved in self.mech_loss. """ m = self.m c = self.c k = self.k - frequencies = self.frequencies number_of_nodes = len(self.mesh_data.nodes) number_of_elements = len(self.mesh_data.elements) element_order = self.mesh_data.element_order points_per_element = int(1/2*(element_order+1)*(element_order+2)) + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) u = np.zeros( - (3*number_of_nodes, len(frequencies)), + (len(frequencies), 3*number_of_nodes), dtype=np.complex128 ) - q = np.zeros((len(frequencies)), dtype=np.complex128) mech_loss = np.zeros( - (number_of_elements, len(frequencies)), + (len(frequencies), number_of_elements), dtype=np.complex128 ) - m, c, k = apply_dirichlet_bc( + m, c, k = mat_apply_dbcs( m, c, k, - self.dirichlet_nodes + dirichlet_nodes ) volumes = calculate_volumes( @@ -305,39 +233,28 @@ def solve_frequency( + 1j*angular_frequency*c + k ) - - f = self.get_load_vector( - self.dirichlet_nodes, - self.dirichlet_values, - frequency_index - ) - u[:, frequency_index] = slin.spsolve(a, f) - - if electrode_elements is not None: - q[frequency_index] = calculate_charge( - u[:, frequency_index], - self.material_manager, - electrode_elements, - electrode_normals, - self.mesh_data.nodes, - self.mesh_data.element_order, - True - ) + # Apply dirichlet bc for f + f = np.zeros(3*number_of_nodes) + f[dirichlet_nodes] = dirichlet_values[:, frequency_index] + + u[frequency_index, :] = slin.spsolve(a, f) # Calculate mech_loss for every element - for element_index, element in enumerate(self.mesh_data.elements): - node_points = self.node_points[element_index] - - # Get field values - u_e = np.zeros(2*points_per_element, dtype=np.complex128) - for i in range(points_per_element): - u_e[2*i] = u[2*element[i], frequency_index] - u_e[2*i+1] = u[2*element[i]+1, frequency_index] - - if calculate_mech_loss: - mech_loss[element_index, frequency_index] = ( - loss_integral_scs( + if calculate_mech_loss: + for element_index, element in enumerate( + self.mesh_data.elements + ): + node_points = self.node_points[element_index] + + # Get field values + u_e = np.zeros(2*points_per_element, dtype=np.complex128) + for i in range(points_per_element): + u_e[2*i] = u[frequency_index, 2*element[i]] + u_e[2*i+1] = u[frequency_index, 2*element[i]+1] + + mech_loss[frequency_index, element_index] = ( + integral_loss_scs( node_points, u_e, self.material_manager.get_elasticity_matrix( @@ -356,33 +273,6 @@ def solve_frequency( print(f"Frequency step {frequency_index} finished") self.u = u - self.q = q - self.mech_loss = mech_loss - - def get_load_vector( - self, - dirichlet_nodes: npt.NDArray, - dirichlet_values: npt.NDArray, - frequency_index: int - ) -> npt.NDArray: - """Calculates the load vector (right hand side) vector for the - simulation. - - Parameters: - dirichlet_nodes: Nodes at which the dirichlet value shall be set. - dirichlet_values: Dirichlet value which is set at the corresponding - node. - - Returns: - Right hand side vector for the simulation. - """ - number_of_nodes = len(self.mesh_data.nodes) - - # Can be initialized to 0 because external load and volume - # charge density is 0. - f = np.zeros(3*number_of_nodes, dtype=np.float64) - - for node, value in zip(dirichlet_nodes, dirichlet_values): - f[node] = value[frequency_index] - return f + if calculate_mech_loss: + self.mech_loss = mech_loss diff --git a/plutho/simulations/piezo_time.py b/plutho/simulations/piezo_time.py new file mode 100644 index 0000000..735c196 --- /dev/null +++ b/plutho/simulations/piezo_time.py @@ -0,0 +1,247 @@ +"""Module for the simulation of time domain piezoelectric systems.""" + + +# Third party libraries +import numpy as np +import numpy.typing as npt +import scipy.sparse as sparse +import scipy.sparse.linalg as slin + +# Local libraries +from .helpers import create_node_points, mat_apply_dbcs +from .integrals import integral_m, integral_ku, integral_kuv, \ + integral_kve +from ..enums import SolverType +from .solver import FEMSolver +from ..mesh.mesh import Mesh + + +__all__ = [ + "PiezoTime" +] + + +class PiezoTime(FEMSolver): + """Class for the simulation of time domain piezoelectric systems. + + Attributes: + node_points: List of node points per elements. + m: Sparse mass matrix. + c: Sparse damping matrix. + k: Sparse stiffness matrix. + """ + # Internal simulation data + node_points: npt.NDArray + + # FEM matrices + m: sparse.lil_array + c: sparse.lil_array + k: sparse.lil_array + + def __init__(self, simulation_name: str, mesh: Mesh): + super().__init__(simulation_name, mesh) + + self.solver_type = SolverType.PiezoTime + + def assemble(self): + """Assembles the FEM matrices based on the set mesh_data and + material_data. + The matrices are stored in self.m, self.c, self.k. + """ + # TODO Assembly takes to long rework this algorithm + # Maybe the 2x2 matrix slicing is not very fast + self.material_manager.initialize_materials() + nodes = self.mesh_data.nodes + elements = self.mesh_data.elements + element_order = self.mesh_data.element_order + self.node_points = create_node_points(nodes, elements, element_order) + + number_of_nodes = len(nodes) + mu = sparse.lil_matrix( + (2*number_of_nodes, 2*number_of_nodes), + dtype=np.float64 + ) + ku = sparse.lil_matrix( + (2*number_of_nodes, 2*number_of_nodes), + dtype=np.float64 + ) + kuv = sparse.lil_matrix( + (2*number_of_nodes, number_of_nodes), + dtype=np.float64 + ) + kve = sparse.lil_matrix( + (number_of_nodes, number_of_nodes), + dtype=np.float64 + ) + + for element_index, element in enumerate(self.mesh_data.elements): + node_points = self.node_points[element_index] + + # TODO Check if its necessary to calculate all integrals + # --> Dirichlet nodes could be leaved out? + # Multiply with jac_det because its integrated with respect to + # local coordinates. + # Mutiply with 2*pi because theta is integrated from 0 to 2*pi + mu_e = ( + self.material_manager.get_density(element_index) + * integral_m(node_points, element_order) + * 2 * np.pi + ) + ku_e = ( + integral_ku( + node_points, + self.material_manager.get_elasticity_matrix(element_index), + element_order + ) + * 2 * np.pi + ) + kuv_e = ( + integral_kuv( + node_points, + self.material_manager.get_piezo_matrix(element_index), + element_order + ) + * 2 * np.pi + ) + kve_e = ( + integral_kve( + node_points, + self.material_manager.get_permittivity_matrix( + element_index + ), + element_order + ) + * 2 * np.pi + ) + + # Now assemble all element matrices + for local_p, global_p in enumerate(element): + for local_q, global_q in enumerate(element): + # Mu_e is returned as a 3x3 matrix (should be 6x6) + # because the values are the same. + # Only the diagonal elements have values. + mu[2*global_p, 2*global_q] += mu_e[local_p][local_q] + mu[2*global_p+1, 2*global_q+1] += mu_e[local_p][local_q] + + # Ku_e is a 6x6 matrix and 2x2 matrices are sliced out + # of it. + ku[2*global_p:2*global_p+2, 2*global_q:2*global_q+2] += \ + ku_e[2*local_p:2*local_p+2, 2*local_q:2*local_q+2] + + # KuV_e is a 6x3 matrix and 2x1 vectors are sliced out. + # [:, None] converts to a column vector. + kuv[2*global_p:2*global_p+2, global_q] += \ + kuv_e[2*local_p:2*local_p+2, local_q][:, None] + + # KVe_e is a 3x3 matrix and therefore the values can + # be set directly. + kve[global_p, global_q] += kve_e[local_p, local_q] + + # Calculate damping matrix + # Currently in the simulations only one material is used + # TODO Add algorithm to calculate cu for every element on its own + # in order to use multiple materials + cu = ( + self.material_manager.get_alpha_m(0) * mu + + self.material_manager.get_alpha_k(0) * ku + ) + + # Calculate block matrices + zeros1x1 = np.zeros((number_of_nodes, number_of_nodes)) + zeros2x1 = np.zeros((2*number_of_nodes, number_of_nodes)) + zeros1x2 = np.zeros((number_of_nodes, 2*number_of_nodes)) + + m = sparse.bmat([ + [mu, zeros2x1], + [zeros1x2, zeros1x1], + ]) + c = sparse.bmat([ + [cu, zeros2x1], + [zeros1x2, zeros1x1], + ]) + k = sparse.bmat([ + [ku, kuv], + [kuv.T, -1*kve] + ]) + + self.m = m.tolil() + self.c = c.tolil() + self.k = k.tolil() + + def simulate( + self, + delta_t: float, + number_of_time_steps: int, + gamma: float, + beta: float + ): + """Runs the time domain simulation. Dirichlet boundary conditions must + be set beforehand. + The resulting field is stored in self.u, which contains the fields u_r, + u_z and phi. + + Parameters: + delta_t: Distance between two time steps. + number_of_time_steps: Total number of time steps for the + simulation. + gamma: Newmark integration parameter. + beta: Newmark integration parameter. + """ + m = self.m + c = self.c + k = self.k + + node_count = len(self.mesh_data.nodes) + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) + + # Init arrays + # Displacement u which is calculated + u = np.zeros((number_of_time_steps, 3*node_count), dtype=np.float64) + # Displacement u derived after t (du/dt) + v = np.zeros((number_of_time_steps, 3*node_count), dtype=np.float64) + # v derived after u (d^2u/dt^2) + a = np.zeros((number_of_time_steps, 3*node_count), dtype=np.float64) + + m, c, k = mat_apply_dbcs( + m, + c, + k, + dirichlet_nodes + ) + + k_star = (k+gamma/(beta*delta_t)*c+1/(beta*delta_t**2)*m).tocsr() + + print("Starting simulation") + for time_index in range(number_of_time_steps-1): + # Calculate load vector and add dirichlet boundary conditions + f = np.zeros(3*node_count) + f[dirichlet_nodes] = dirichlet_values[:, time_index+1] + + # Perform Newmark method + # Predictor step + u_tilde = ( + u[time_index, :] + + delta_t*v[time_index, :] + + delta_t**2/2*(1-2*beta)*a[time_index, :] + ) + v_tilde = ( + v[time_index, :] + + (1-gamma)*delta_t*a[time_index, :] + ) + + # Solve for next time step + u[time_index+1, :] = slin.spsolve( + k_star, ( + f + - c*v_tilde + + (1/(beta*delta_t**2)*m + + gamma/(beta*delta_t)*c)*u_tilde)) + # Perform corrector step + a[time_index+1, :] = (u[time_index+1, :]-u_tilde)/(beta*delta_t**2) + v[time_index+1, :] = v_tilde + gamma*delta_t*a[time_index+1, :] + + if (time_index + 1) % 100 == 0: + print(f"Finished time step {time_index+1}") + + self.u = u diff --git a/plutho/simulations/solver.py b/plutho/simulations/solver.py new file mode 100644 index 0000000..d180f67 --- /dev/null +++ b/plutho/simulations/solver.py @@ -0,0 +1,394 @@ +"""Base module for all linear simulations.""" + +# Python standard libraries +import os +import configparser +from abc import ABC, abstractmethod +from typing import Dict, Union, List + +# Third party libraries +import numpy as np +import numpy.typing as npt + + +# Local libraries +from ..enums import SolverType, FieldType +from ..materials import MaterialManager, MaterialData +from ..mesh.mesh import Mesh, MeshData +from .helpers import calculate_charge + + +__all__ = [ + "FEMSolver" +] + + +class FEMSolver(ABC): + """Base class for all linear solvers. The assemble() and simulation() + functions must be implemented in the child class + + Attributes: + simulation_name: Name of the simulation. + mesh: Mesh object, stores all the mesh information. + material_manager: MaterialManager object, stores all the material + informations. + dirichlet_nodes: List of nodes for which a dirichlet bc is applied. + dirichlet_values: Values to which the fields are set at the respective + nodes + boundary_conditions: List of all implemented boundary conditions. Is + used when the simulation settings are saved. + solver_type: Type of the current solver. + + Parameters: + simulation_name: Name of the simulation. + mesh: mesh object, stores all the mesh information. + """ + simulation_name: str + mesh: Mesh + material_manager: MaterialManager + dirichlet_nodes: List + dirichlet_values: List + boundary_conditions: List + solver_type: Union[SolverType, None] + + def __init__(self, simulation_name: str, mesh: Mesh): + self.simulation_name = simulation_name + self.mesh = mesh + self.solver_type = None + + nodes, elements = mesh.get_mesh_nodes_and_elements() + self.mesh_data = MeshData(nodes, elements, mesh.element_order) + self.material_manager = MaterialManager(len(elements)) + + self.dirichlet_nodes = [] + self.dirichlet_values = [] + self.boundary_conditions = [] + + self.u = None + self.q = None + + def add_material( + self, + material_name: str, + material_data: MaterialData, + physical_group_name: str + ): + """Adds a material to the simulation + + Parameters: + material_name: Name of the material. + material_data: Material data. + physical_group_name: Name of the physical group for which the + material shall be set. If this is None or "" the material will + be set for the whole model. + """ + if physical_group_name is None or physical_group_name == "": + element_indices = np.arange(len(self.mesh_data.elements)) + else: + element_indices = self.mesh.get_elements_by_physical_groups( + [physical_group_name] + )[physical_group_name] + + self.material_manager.add_material( + material_name, + material_data, + physical_group_name, + element_indices + ) + + def add_dirichlet_bc( + self, + field_type: FieldType, + physical_group_name: str, + values: npt.NDArray + ): + """Adds a dirichlet boundary condition (dirichlet bc) to the + simulation. The given values are for each time step. Therefore each + value is applied to every node from the given physical_group. + + Parameters: + field_type: The type of field for which the bc shall be set. + physical_group_name: Name of the physical group for which the + boundary condition shall be set. + values: List of values per time step. Value of the bc for each time + step. The value for one time step is applied to every element + equally. Length: number_of_time_step. + """ + if self.solver_type is None: + raise ValueError("Solver type is not set") + + match self.solver_type: + case SolverType.ThermoTime: + if (field_type is FieldType.U_R or + field_type is FieldType.U_Z or + field_type is FieldType.PHI): + raise ValueError( + f"Unknown variable type {field_type} given for" + f"simulation type {self.solver_type}" + ) + case SolverType.PiezoFreq: + if field_type is FieldType.THETA: + raise ValueError( + f"Unknown variable type {field_type} given for" + f"simulation type {self.solver_type}" + ) + + # Save boundary condition for serialization + self.boundary_conditions.append({ + "field_type": field_type.name, + "physical_group_name": physical_group_name, + "values": values.tolist() + }) + + # Apply the given values for all nodes from the given physical_group + node_indices = self.mesh.get_nodes_by_physical_groups( + [physical_group_name] + )[physical_group_name] + + number_of_nodes = len(self.mesh_data.nodes) + for node_index in node_indices: + real_index = 0 + + # Depending on the variable type and the simulation types + # the corresponding field variable may be found at different + # positions of the solution vector + if field_type is FieldType.PHI: + real_index = 2*number_of_nodes+node_index + elif field_type is FieldType.THETA: + if self.solver_type is SolverType.ThermoTime: + real_index = node_index + else: + real_index = 3*number_of_nodes+node_index + elif field_type is FieldType.U_R: + real_index = 2*node_index + elif field_type is FieldType.U_Z: + real_index = 2*node_index+1 + + self.dirichlet_nodes.append(real_index) + self.dirichlet_values.append(values) + + def clear_dirichlet_bcs(self): + """Resets the dirichlet boundary conditions.""" + self.boundary_conditions = [] + self.dirichlet_nodes = [] + self.dirichlet_values = [] + + def save_simulation_settings(self, directory: str, prefix: str = ""): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation settings in it. + + Parameters: + directory: Path to the directory in which the simulation folder is + created. + prefix: Possible prefix, which is prepend to the file names. + """ + if self.solver_type is None: + raise ValueError("Simulation settings can not be saved, since no \ + simulation is configured yet.") + + wd = os.path.join(directory, self.simulation_name) + + if not os.path.isdir(wd): + os.makedirs(wd) + + settings = configparser.ConfigParser() + settings["general"] = { + "simulation_name": self.simulation_name, + "solver_type": self.solver_type.value, + "mesh_file_path": self.mesh.mesh_file_path, + "mesh_element_order": str(self.mesh_data.element_order) + } + settings["materials"] = self._get_materials_settings() + settings["boundary_conditions"] = self._get_bc_settings() + + # Save config file + if prefix != "": + prefix += "_" + + config_file = os.path.join(wd, f"{prefix}settings.cfg") + with open(config_file, "w", encoding="UTF-8") as fd: + settings.write(fd) + + def save_simulation_results(self, directory: str, prefix: str = ""): + """Creates a folder 'simulation_name' in the given directory and saves + the current simulation results in it. + + Parameters: + directory: Path to the directory in which the simulation folder is + created. + prefix: Possible prefix, which is prepend to the file names. + """ + wd = os.path.join(directory, self.simulation_name) + + if not os.path.isdir(wd): + os.makedirs(wd) + + if prefix != "": + prefix += "_" + + if self.u is not None: + u_path = os.path.join(wd, f"{prefix}u.npy") + np.save(u_path, self.u) + + if self.q is not None: + q_path = os.path.join(wd, f"{prefix}q.npy") + np.save(q_path, self.q) + + match self.solver_type: + case SolverType.ThermoPiezoTime: + mech_loss_path = os.path.join(wd, f"{prefix}mech_loss.npy") + np.save(mech_loss_path, self.mech_loss) + case SolverType.ThermoTime: + theta_path = os.path.join(wd, f"{prefix}theta.npy") + np.save(theta_path, self.theta) + case _: + pass + + @classmethod + def load_simulation_settings(cls, simulation_folder: str): + """Loads the simulation settings from the given folder. + + Parameters: + simulation_folder: Path to the simulation folder, which shall be + loaded (must contain a *.cfg file). + + Returns: + The loaded simulation class object (some child of the FEMSolver + class). + """ + simulation_name = os.path.basename(simulation_folder) + config_path = os.path.join(simulation_folder, f"{simulation_name}.cfg") + settings = configparser.ConfigParser() + settings.read(config_path) + simulation_name = settings["general"]["simulation_name"] + mesh_file_path = settings["general"]["mesh_file_path"] + mesh_element_order = int(settings["general"]["mesh_element_order"]) + mesh = Mesh(mesh_file_path, mesh_element_order) + + # Check for solver type + solver_type = None + try: + solver_type = SolverType(settings["general"]["solver_type"]) + except ValueError: + raise ValueError( + f"Unknown solver type: {settings['general']['solver_type']}" + ) + + # Create simulation + sim = cls(simulation_name, mesh) + + if solver_type is not None and sim.solver_type is not solver_type: + raise ValueError( + "Wrong class used to instantiate a solver of type " + f"{solver_type}" + ) + + # Load materials + material_settings = settings["materials"] + for material_name, material in material_settings.items(): + material_data = MaterialData.from_dict(material["material_data"]) + physical_group_name = material["physical_group_name"] + + sim.add_material( + material_name, + material_data, + physical_group_name + ) + + # Load boundary conditions + bc_settings = settings["boundary_conditions"] + for bc in bc_settings: + field_type = FieldType(bc["field_type"]) + physical_group_name = bc["physical_group_name"] + values = np.array(list(bc["values"])) + + sim.add_dirichlet_bc(field_type, physical_group_name,values) + + return sim + + def _get_materials_settings(self) -> Dict: + """Local function. Is used to parse the material parameters to a + dictionary. + """ + material_settings = {} + + for material in self.material_manager.materials: + material_settings[material.material_name] = { + "material_data": material.to_dict(), + "physical_group_name": material.physical_group_name + } + + return material_settings + + def _get_bc_settings(self) -> Dict: + """Local function. Is used to parse the boundary conditions to a + dirctionary + """ + bc_settings = {} + + for index, bc in enumerate(self.boundary_conditions): + bc_settings[index] = bc + + return bc_settings + + def calculate_charge( + self, + electrode_name: str, + is_complex: bool = False + ): + """Calculates the charge for the displacement u stored in the current + object -> self.u. Therefore the simulate() function must be called + beforehand or a self.u is set manually. + + Parameters: + electrode_name: Name of the boundary on which the charge is + calculated, typically the electrode. + is_complex: Set to true if the displacement and thus the charge + can be of a complex type. E.g. in a frequency domain + simulation. + """ + if self.u is None: + raise ValueError("Cannot calculate charge since no simulation \ + has been done") + + # Get electrode elements and electrode normals + electrode_elements = self.mesh.get_elements_by_physical_groups( + [electrode_name] + )[electrode_name] + # TODO This must be calculated dynamically depending on the direction + # of the normal vector + electrode_normals = np.tile([0, 1], (len(electrode_elements), 1)) + + iteration_count, _ = self.u.shape + + if is_complex: + q = np.zeros(iteration_count, dtype=np.complex128) + else: + q = np.zeros(iteration_count) + + for index in range(iteration_count): + q[index] = calculate_charge( + self.u[index, :], + self.material_manager, + electrode_elements, + electrode_normals, + self.mesh_data.nodes, + self.mesh_data.element_order, + is_complex + ) + + self.q = q + + @abstractmethod + def assemble(self, *args, **kwargs): + """Needs to be implemented by any child class. Implements the matrix + assembly. + """ + raise NotImplementedError("'assembly' function not implemented.") + + @abstractmethod + def simulate(self, *args, **kwargs): + """Needs to be implemented by any child class. Implements the + simulation process + """ + raise NotImplementedError("'simulate' function not implemented.") diff --git a/plutho/simulation/thermo_piezo_time.py b/plutho/simulations/thermo_piezo_time.py similarity index 57% rename from plutho/simulation/thermo_piezo_time.py rename to plutho/simulations/thermo_piezo_time.py index 6473271..5884c21 100644 --- a/plutho/simulation/thermo_piezo_time.py +++ b/plutho/simulations/thermo_piezo_time.py @@ -1,124 +1,42 @@ -"""Main module for the simulation of piezoelectric systems with -thermal field.""" +"""Module for the simulation of thermo-piezoelectric systems in the time +domain.""" # Python standard libraries -from typing import List +from typing import Union + +# Third party libraries import numpy as np import numpy.typing as npt import scipy.sparse.linalg as slin from scipy import sparse # Local libraries -from .base import SimulationData, MeshData, \ - create_node_points, \ - local_to_global_coordinates, b_operator_global, integral_m, \ - integral_ku, integral_kuv, integral_kve, apply_dirichlet_bc, \ - quadratic_quadrature, calculate_volumes, \ - gradient_local_shape_functions_2d, get_avg_temp_field_per_element -from .piezo_time import calculate_charge -from .thermo_time import integral_ktheta, integral_theta_load -from ..materials import MaterialManager - - -def loss_integral_scs( - node_points: npt.NDArray, - u_e_t: npt.NDArray, - u_e_t_minus_1: npt.NDArray, - u_e_t_minus_2: npt.NDArray, - delta_t: float, - elasticity_matrix: npt.NDArray, - element_order: int -): - """Calculates the integral of dS/dt*c*dS/dt over one triangle. Since foward - difference quotient of second oder is used the last 2 values of e_u are - needed. - - Parameters: - node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of - one triangle. - u_e_t: Values of u_e at the current time point. - Format: [u1_r, u1_z, u2_r, u2_z, u3_r, u3_z]. - u_e_t_minus_1: Values of u_e at one time point earlier. Same format - as u_e_t. - u_e_t_minus_2: Values of u_e at two time points earlier. Same format - as u_e_t. - delta_t: Difference between the time steps. - jacobian_inverted_t: Jacobian matrix inverted and transposed, needed - for calculation of global derivatives. - elasticity_matrix: Elasticity matrix for the current element - (c matrix). - element_order: Order of the shape functions. - """ - def inner(s, t): - dn = gradient_local_shape_functions_2d(s, t, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - jacobian_det = np.linalg.det(jacobian) - - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - b_opt = b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - - s_e = np.dot(b_opt, u_e_t) - s_e_t_minus_1 = np.dot(b_opt, u_e_t_minus_1) - s_e_t_minus_2 = np.dot(b_opt, u_e_t_minus_2) - dt_s = (3*s_e-4*s_e_t_minus_1+s_e_t_minus_2)/(2*delta_t) - - return np.dot( - dt_s.T, - np.dot( - elasticity_matrix.T, - dt_s - ) - ) * r * jacobian_det +from .solver import FEMSolver +from ..enums import SolverType +from ..mesh import Mesh +from .helpers import create_node_points, mat_apply_dbcs, calculate_volumes, \ + get_avg_temp_field_per_element +from .integrals import integral_m, integral_ku, integral_kuv, \ + integral_kve, integral_ktheta, integral_theta_load, integral_loss_scs_time - return quadratic_quadrature(inner, element_order) +__all__ = [ + "ThermoPiezoTime" +] -class ThermoPiezoSimTime: - """Class for the coupled simulation of mechanical-electric and - thermal field. - Parameters: - mesh_data: MeshData format. - material_data: MaterialData format. - simulation_data: SimulationData format. +class ThermoPiezoTime(FEMSolver): + """Class for solving time domain themo-piezoelectric systems. Attributes: - mesh_data: Contains the mesh information. - material_data: Contains the information about the materials. - simulation_data: Contains the information about the simulation. - dirichlet_nodes: List of nodes for which dirichlet values are set. - dirichlet_values: List of values of the corresponding dirichlet nodes - per time step. - m: Mass matrix. - c: Damping matrix. - k: Stiffness matrix. - u: Resulting vector of size 4*number_of_nodes containing u_r, u_z, v - and theta. u[node_index, time_index]. An offset needs to be - added to the node_index in order to access the different field - paramteters: - u_r: 0, - u_z: 1*number_of_nodes, - v: 2*number_of_nodes, - theta: 3*number_of_nodes - q: Resulting charges for each time step q[time_index]. - mech_loss: Mechanical loss for each element per time step - mech_loss[element_index, time_index]. + node_points: List of node points per elements. + m: Sparse mass matrix. + c: Sparse damping matrix. + k: Sparse stiffness matrix. + mech_loss: Mechanical loss field. """ - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - simulation_data: SimulationData - - # Dirichlet boundary condition - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray + # Internal simulation data + node_points: npt.NDArray # FEM matrices m: sparse.lil_array @@ -126,22 +44,12 @@ class ThermoPiezoSimTime: k: sparse.lil_array # Resulting fields - u: npt.NDArray - q: npt.NDArray mech_loss: npt.NDArray - # Internal simulation data - node_points: npt.NDArray + def __init__(self, simulation_name: str, mesh: Mesh): + super().__init__(simulation_name, mesh) - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - simulation_data: SimulationData - ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.simulation_data = simulation_data + self.solver_type = SolverType.ThermoPiezoTime def assemble(self): """Assembles the FEM matrices based on the set mesh_data and @@ -150,6 +58,7 @@ def assemble(self): """ # TODO Assembly takes long rework this algorithm? # Maybe the 2x2 matrix slicing is not very fast + self.material_manager.initialize_materials() nodes = self.mesh_data.nodes elements = self.mesh_data.elements element_order = self.mesh_data.element_order @@ -284,39 +193,37 @@ def assemble(self): self.c = c.tolil() self.k = k.tolil() - def solve_time( + def simulate( self, - electrode_elements: npt.NDArray, - electrode_normals: npt.NDArray, - theta_start = None + delta_t: float, + number_of_time_steps: int, + gamma: float, + beta: float, + theta_start: Union[None, npt.NDArray] = None ): - """Runs the simulation using the assembled m, c and k matrices as well - as the set excitation. - Calculates the displacement field, potential field and the thermal - field of the given model (all saved in u). - Calculates the electric charge in the elctrode and the - mechanical losses of the whole body. - Saves u[node_index, time_index], q[time_index] and - mechanical_losses[element_index, time_index]. + """Runs the time domain simulation. Dirichlet boundary conditions must + be set beforehand. + The resulting displacement field is stored in self.u, which contains + the fields u_r, u_z, phi and theta. + The mechanical losses are calculated in each time step and used as a + source for the thermal field. Parameters: - electrode_elements: List of all elements which are in - the electrode. - electrode_normals: List of normal vectors corresponding to the - electrode_elements list. - set_symmetric_bc: True if the symmetric boundary condition for u - shall be set. False otherwise. + delta_t: Distance between two time steps. + number_of_time_steps: Total number of time steps for the + simulation. + gamma: Newmark integration parameter. + beta: Newmark integration parameter. + theta_start: Possible initial field for theta. """ - number_of_time_steps = self.simulation_data.number_of_time_steps - beta = self.simulation_data.beta - gamma = self.simulation_data.gamma - delta_t = self.simulation_data.delta_t elements = self.mesh_data.elements nodes = self.mesh_data.nodes element_order = self.mesh_data.element_order points_per_element = int(1/2*(element_order+1)*(element_order+2)) - number_of_elements = len(elements) - number_of_nodes = len(nodes) + node_count = len(nodes) + element_count = len(elements) + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) m = self.m c = self.c @@ -324,16 +231,13 @@ def solve_time( # Init arrays # Displacement u which is calculated - u = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) + u = np.zeros((number_of_time_steps, 4*node_count), dtype=np.float64) # Displacement u derived after t (du/dt) - v = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) + v = np.zeros((number_of_time_steps, 4*node_count), dtype=np.float64) # v derived after u (d^2u/dt^2) - a = np.zeros((m.shape[0], number_of_time_steps), dtype=np.float64) - - # Charge calculated during simulation (for impedance) - q = np.zeros(number_of_time_steps, dtype=np.float64) + a = np.zeros((number_of_time_steps, 4*node_count), dtype=np.float64) - m, c, k = apply_dirichlet_bc( + m, c, k = mat_apply_dbcs( m, c, k, @@ -344,7 +248,7 @@ def solve_time( # Mechanical loss calculated during simulation (for thermal field) mech_loss = np.zeros( - shape=(number_of_elements, number_of_time_steps), + shape=(number_of_time_steps, element_count), dtype=np.float64 ) @@ -354,11 +258,11 @@ def solve_time( ) if theta_start is not None: - if len(theta_start) != number_of_nodes: + if len(theta_start) != node_count: raise ValueError( "theta start must have the size of number of nodes" ) - u[3*number_of_nodes:, 0] = theta_start + u[3*node_count:, 0] = theta_start print("Starting simulation") for time_index in range(number_of_time_steps-1): @@ -366,7 +270,7 @@ def solve_time( # material parameters are used if self.material_manager.is_temperature_dependent: temp_field_per_elements = get_avg_temp_field_per_element( - u[3*number_of_nodes:, time_index], + u[time_index, 3*node_count:], self.mesh_data.elements ) update = self.material_manager.update_temperature( @@ -376,7 +280,7 @@ def solve_time( # Assemble the matrices with new parameters print(f"Assemble new (time step: {time_index})") self.assemble() - m, c, k = apply_dirichlet_bc( + m, c, k = mat_apply_dbcs( self.m, self.c, self.k, @@ -389,45 +293,32 @@ def solve_time( ).tocsr() # Calculate load vector and add dirichlet boundary conditions - f = self.get_load_vector( - self.dirichlet_nodes, - self.dirichlet_values[:, time_index+1], - mech_loss[:, time_index] - ) + f = self.calculate_load_vector(mech_loss[time_index, :]) + f[dirichlet_nodes] = dirichlet_values[:, time_index+1] # Perform Newmark method # Predictor step u_tilde = ( - u[:, time_index] - + delta_t*v[:, time_index] - + delta_t**2/2*(1-2*beta)*a[:, time_index] + u[time_index, :] + + delta_t*v[time_index, :] + + delta_t**2/2*(1-2*beta)*a[time_index, :] ) v_tilde = ( - v[:, time_index] - + (1-gamma)*delta_t*a[:, time_index] + v[time_index, :] + + (1-gamma)*delta_t*a[time_index, :] ) # Solve for next time step - u[:, time_index+1] = slin.spsolve( + u[time_index+1, :] = slin.spsolve( k_star, ( f - c*v_tilde + (1/(beta*delta_t**2)*m - + gamma/(beta*delta_t)*c)*u_tilde) + + gamma/(beta*delta_t)*c)*u_tilde) ) # Perform corrector step - a[:, time_index+1] = (u[:, time_index+1]-u_tilde)/(beta*delta_t**2) - v[:, time_index+1] = v_tilde + gamma*delta_t*a[:, time_index+1] - - if electrode_elements is not None: - q[time_index+1] = calculate_charge( - u[:, time_index+1], - self.material_manager, - electrode_elements, - electrode_normals, - nodes, - element_order - ) + a[time_index+1, :] = (u[time_index+1, :]-u_tilde)/(beta*delta_t**2) + v[time_index+1, :] = v_tilde + gamma*delta_t*a[time_index+1, :] # Calculate power_loss for element_index, element in enumerate(elements): @@ -438,18 +329,18 @@ def solve_time( u_e_t_minus_1 = np.zeros(2*points_per_element) u_e_t_minus_2 = np.zeros(2*points_per_element) for i in range(points_per_element): - u_e[2*i] = u[2*element[i], time_index+1] - u_e[2*i+1] = u[2*element[i]+1, time_index+1] - u_e_t_minus_1[2*i] = u[2*element[i], time_index] - u_e_t_minus_1[2*i+1] = u[2*element[i]+1, time_index] - u_e_t_minus_2[2*i] = u[2*element[i], time_index-1] - u_e_t_minus_2[2*i+1] = u[2*element[i]+1, time_index-1] + u_e[2*i] = u[time_index+1, 2*element[i]] + u_e[2*i+1] = u[time_index+1, 2*element[i]+1] + u_e_t_minus_1[2*i] = u[time_index, 2*element[i]] + u_e_t_minus_1[2*i+1] = u[time_index, 2*element[i]+1] + u_e_t_minus_2[2*i] = u[time_index-1, 2*element[i]] + u_e_t_minus_2[2*i+1] = u[time_index-1, 2*element[i]+1] # The mech loss of the element is divided by the volume # because it must be a power density. if time_index > 0: - mech_loss[element_index, time_index+1] = ( - loss_integral_scs( + mech_loss[time_index+1, element_index] = ( + integral_loss_scs_time( node_points, u_e, u_e_t_minus_1, @@ -468,34 +359,23 @@ def solve_time( if (time_index + 1) % 100 == 0: print(f"Finished time step {time_index+1}") - self.q = q self.u = u self.mech_loss = mech_loss - def get_load_vector( + def calculate_load_vector( self, - dirichlet_nodes: npt.NDArray, - values: npt.NDArray, mech_loss_density: npt.NDArray ) -> npt.NDArray: """Calculates the load vector (right hand side) vector for the - simulation. + simulation based on the mechanical loss density. Parameters: - nodes: Nodes where the dirichlet bc is set. - values: Values at which the corresponding nodes are set. mech_loss_density: Power loss density of every element of the current time step. Returns: Right hand side vector for the simulation. """ - # For u and v the load vector is set to the corresponding dirichlet - # values given by values_u and values_v - # at the nodes nodes_u and nodes_v since there is no inner charge and - # no forces given. - # For the temperature field the load vector represents the temperature - # sources given by the mechanical losses. nodes = self.mesh_data.nodes element_order = self.mesh_data.element_order points_per_element = int(1/2*(element_order+1)*(element_order+2)) @@ -505,10 +385,6 @@ def get_load_vector( # density is 0. f = np.zeros(4*number_of_nodes, dtype=np.float64) - # Set dirichlet values for given ndoes and values - for node, value in zip(dirichlet_nodes, values): - f[node] = value - # Calculation for theta load. # It needs to be assembled every step and every element since the # power is time- and position-dependent diff --git a/plutho/simulation/thermo_time.py b/plutho/simulations/thermo_time.py similarity index 51% rename from plutho/simulation/thermo_time.py rename to plutho/simulations/thermo_time.py index 67343b2..882aa82 100644 --- a/plutho/simulation/thermo_time.py +++ b/plutho/simulations/thermo_time.py @@ -1,7 +1,9 @@ -"""Extra module for the simulation of heat conduction over time.""" +"""Module for the simulation of a heat conduction problem over time. +Additonaly heat convection can be set.""" + # Python standard libraries -from typing import List +from typing import List, Union # Third party libraries import numpy as np @@ -10,148 +12,57 @@ from scipy import sparse # Local libraries -from .base import SimulationData, MeshData, \ - local_shape_functions_2d, gradient_local_shape_functions_2d, \ - local_to_global_coordinates, integral_m, \ - quadratic_quadrature, line_quadrature, create_node_points, \ - apply_dirichlet_bc, get_avg_temp_field_per_element -from ..materials import MaterialManager - - -def integral_heat_flux( - node_points: npt.NDArray, - heat_flux: npt.NDArray, - element_order: int -): - """Integrates the heat flux using the shape functions. - - Parameters: - node_points: List of node points [[x1, x2, ..], [y1, y2, ..]]. - heat_flux: Heat flux at the points. - element_order: Order of the shape functions. - - Returns: - npt.NDArray heat flux integral on each point. - """ - def inner(s): - n = local_shape_functions_2d(s, 0, element_order) - r = local_to_global_coordinates(node_points, s, 0, element_order)[0] - - return n*heat_flux*r - - return line_quadrature(inner, element_order) - - -def integral_ktheta( - node_points: npt.NDArray, - element_order: int -): - """Calculates the Ktheta integral. - - Parameters: - node_points: List of node points [[x1, x2, ..], [y1, y2, ..]]. - element_order: Order of the shape functions. - - Returns: - npt.NDArray: 3x3 Ktheta matrix for the given element. - """ - def inner(s, t): - dn = gradient_local_shape_functions_2d(s, t, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_inverted_t = np.linalg.inv(jacobian).T - jacobian_det = np.linalg.det(jacobian) - - global_dn = np.dot(jacobian_inverted_t, dn) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - - return np.dot(global_dn.T, global_dn)*r*jacobian_det - - return quadratic_quadrature(inner, element_order) - - -def integral_theta_load( - node_points: npt.NDArray, - mech_loss: float, - element_order: int -): - """Returns the load value for the temperature field (f) for the specific - element. - - Parameters: - node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of - one triangle. - point_loss: Loss power on each node (heat source). - element_order: Order of the shape functions. - - Returns: - npt.NDArray: f vector value at the specific ndoe - """ - def inner(s, t): - dn = gradient_local_shape_functions_2d(s, t, element_order) - jacobian = np.dot(node_points, dn.T) - jacobian_det = np.linalg.det(jacobian) - - n = local_shape_functions_2d(s, t, element_order) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] +from .integrals import integral_heat_flux, integral_ktheta, \ + integral_theta_load, integral_m +from .solver import FEMSolver +from .helpers import create_node_points, mat_apply_dbcs, \ + get_avg_temp_field_per_element +from ..enums import SolverType +from ..mesh import Mesh - return n*mech_loss*r*jacobian_det - return quadratic_quadrature(inner, element_order) +__all__ = [ + "ThermoTime" +] -class ThermoSimTime: - """Simulates a heat conduction problem for the given mesh, material and - simulation information. +class ThermoTime(FEMSolver): + """Class for the simulation of time domain heat conduction problems. Attributes: - mesh_data: Contains the mesh information. - material_manager: Organizes the material parameters. - simulation_data: Contains the simulation data. - theta: Temperature field defined for every node. - f: Load vector for the temperature field contains the information - from the given losses. - """ - # Simulation settings - mesh_data: MeshData - material_manager: MaterialManager - simulation_data: SimulationData - - # FEM Vectors - theta: npt.NDArray - f: npt.NDArray + node_points: List of node points per elements. + c: Sparse damping matrix. + k: Sparse stiffness matrix. + theta: Resulting thermal field. + enable_convection_bc: True if a convection boundary condition is set. + convective_b_e: List of elements for which the convection boundary + condition is applied. + convective_alpha: Heat transfer coefficient of the convection. + convective_outer_temp: Temperature of the material outside of the + simulated material. + """ + # Internal simulation data + node_points: npt.NDArray # FEM Matrices - k: sparse.lil_array c: sparse.lil_array + k: sparse.lil_array + + # Resulting arrays + theta: npt.NDArray # Convective bc enable_convection_bc: bool - convective_b_e: List[int] + convective_b_e: npt.NDArray convective_alpha: float convective_outer_temp: float - # Internal simulation data - node_points: npt.NDArray + def __init__(self, simulation_name: str, mesh: Mesh): + super().__init__(simulation_name, mesh) - # Dirichlet data - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray - - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - simulation_data: SimulationData): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.simulation_data = simulation_data - self.f = np.zeros( - ( - len(self.mesh_data.nodes), - self.simulation_data.number_of_time_steps - ) - ) + self.solver_type = SolverType.ThermoTime self.enable_convection_bc = False + self.f = None def assemble(self): """Assembles the FEM matrices based on the set mesh_data and @@ -160,16 +71,17 @@ def assemble(self): """ # TODO Assembly takes long rework this algorithm? # Maybe the 2x2 matrix slicing is not very fast + self.material_manager.initialize_materials() nodes = self.mesh_data.nodes elements = self.mesh_data.elements element_order = self.mesh_data.element_order self.node_points = create_node_points(nodes, elements, element_order) number_of_nodes = len(nodes) - c = sparse.lil_matrix( + c = sparse.lil_array( (number_of_nodes, number_of_nodes), dtype=np.float64) - k = sparse.lil_matrix( + k = sparse.lil_array( (number_of_nodes, number_of_nodes), dtype=np.float64) @@ -206,7 +118,7 @@ def set_constant_volume_heat_source( number_of_time_steps: float ): """Sets the excitation for the heat conduction simulation. - The mech_loss_density are set for every time step. + The mech_loss_density is set for every time step. Parameters: mech_loss_density: The mechanical losses for each element. They @@ -237,7 +149,10 @@ def set_constant_volume_heat_source( f[global_p] += f_theta_e[local_p] # Repeat it for n time steps - self.f += np.tile(f.reshape(-1, 1), (1, number_of_time_steps)) + if self.f is None: + self.f = np.tile(f.reshape(-1, 1), (1, number_of_time_steps)) + else: + self.f += np.tile(f.reshape(-1, 1), (1, number_of_time_steps)) def _calculate_convection_bc( self, @@ -256,19 +171,21 @@ def _calculate_convection_bc( Parameters: nodes: All nodes in the mesh. - elements: Boundary elements at which the convection boundary is - added. + boundary_elements: Boundary elements at which the convection + boundary is added. theta: Current temperature field at all nodes. alpha: Heat transfer coefficient. outer_temperature: Temperature outside of the model. + element_order: Element order of the mesh. """ + # This function is not yet tested points_per_element = int(1/2*(element_order+1)*(element_order+2)) f = np.zeros(len(nodes)) - for element_index, element in enumerate(boundary_elements): + for _, element in enumerate(boundary_elements): node_points = np.zeros(shape=(2, points_per_element)) for node_index in range(points_per_element): - node_points[element_index, :, node_index] = [ + node_points[:, node_index] = [ nodes[element[node_index]][0], nodes[element[node_index]][1] ] @@ -293,12 +210,14 @@ def _calculate_convection_bc( " not implemented for element order > 3" ) - theta_e = np.zeros(boundary_nodes.shape[0]) - for index, bnode in enumerate(boundary_nodes): - theta_e[index] = element[bnode] + theta_e = np.zeros(points_per_element) + theta_e[boundary_nodes] = theta[element[boundary_nodes]] + + outer_temps = np.zeros(points_per_element) + outer_temps[boundary_nodes] = outer_temperature heat_flux = alpha*( - theta_e - np.ones(boundary_nodes.shape[0])*outer_temperature + theta_e - outer_temps ) f_e = ( integral_heat_flux(node_points, heat_flux, element_order) @@ -312,9 +231,9 @@ def _calculate_convection_bc( def set_convection_bc( self, - convective_boundary_elements, - alpha, - outer_temperature + convective_boundary_elements: npt.NDArray, + alpha: float, + outer_temperature: float ): """Adds a convection boundary condition to the current simulation. @@ -329,62 +248,47 @@ def set_convection_bc( self.convective_outer_temp = outer_temperature self.enable_convection_bc = True - def set_dirichlet_bc_load_vector( + def simulate( self, - f: npt.NDArray, - nodes: npt.NDArray, - values: npt.NDArray, - time_step: int - ) -> npt.NDArray: - """Adds the dirichlet values for the given load vector. - - Parameters: - f: Existing load vector. Could be empty or contain already set - boundary conditions, which will be overwritten at - dirichlet nodes. - nodes: Nodes at which the dirichlet bc shall be set. - values: Values which shall be set at the corresponding node. - - Returns: - Right hand side vector for the simulation. - """ - for node, value in zip(nodes, values): - f[node] = value[time_step] - - return f - - def solve_time( - self, - theta_start = None + delta_t: float, + number_of_time_steps: int, + gamma: float, + theta_start: Union[npt.NDArray, int, None] = None ): - """Runs the simulation using the assembled c and k matrices as well - as the set excitation. - Calculates the temperature field and saves it in theta. + """Runs the time domain heat conduction simulation. + The resulting temperature field is stored in self.theta. Parameters: - theta_start: Sets a field for time step 0. + delta_t: Distance between two time steps. + number_of_time_steps: Total number of time steps of the simulation. + gamma: Newmark integration parameter. + theta_start: Sets an initial temperature field. """ - number_of_time_steps = self.simulation_data.number_of_time_steps - gamma = self.simulation_data.gamma - delta_t = self.simulation_data.delta_t - c = self.c k = self.k + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) + node_count = len(self.mesh_data.nodes) + # Init arrays - theta = np.zeros((k.shape[0], number_of_time_steps), dtype=np.float64) - if theta_start is not None: - theta[:, 0] = theta_start - # theta derived after t (du/dt) + theta = np.zeros( + (number_of_time_steps, node_count), + dtype=np.float64 + ) theta_dt = np.zeros( - (k.shape[0], number_of_time_steps), - dtype=np.float64) + (number_of_time_steps, node_count), + dtype=np.float64 + ) + + if theta_start is not None: + theta[0, :] = theta_start - _, c, k = apply_dirichlet_bc( - np.zeros(c.shape), + _, c, k = mat_apply_dbcs( + None, c, k, - self.dirichlet_nodes + dirichlet_nodes ) k_star = (k+1/(gamma*delta_t)*c).tocsr() @@ -394,115 +298,115 @@ def solve_time( print("Starting heat conduction simulation") for time_index in range(number_of_time_steps-1): - # Set dirichlet nodes for the load vector - # The load vector could already contain e.g. energy sources - # which are overwritten at the points where a dirichlet bc is set. - f = self.set_dirichlet_bc_load_vector( - self.f[:, time_index], - self.dirichlet_nodes, - self.dirichlet_values, - time_index+1 - ) + current_f = self.f[:, time_index+1] # Add a convection boundary condition if it set if self.enable_convection_bc: - f += self._calculate_convection_bc( + current_f += self._calculate_convection_bc( self.mesh_data.nodes, self.convective_b_e, - theta[:, time_index], + theta[time_index, :], self.convective_alpha, self.convective_outer_temp, self.mesh_data.element_order ) + if len(dirichlet_nodes) > 0: + current_f[dirichlet_nodes] = dirichlet_values + # Perform Newmark method # Predictor step theta_tilde = ( - theta[:, time_index] - + (1-gamma)*delta_t*theta_dt[:, time_index] + theta[time_index, :] + + (1-gamma)*delta_t*theta_dt[time_index, :] ) # Solve for next time step - theta[:, time_index+1] = slin.spsolve( + theta[time_index+1, :] = slin.spsolve( k_star, ( - f - + 1/(gamma*delta_t)*c*theta_tilde + current_f + + 1/(gamma*delta_t)*c@theta_tilde ) ) # Perform corrector step - theta_dt[:, time_index+1] = \ - (theta[:, time_index+1]-theta_tilde)/(gamma*delta_t) + theta_dt[time_index+1, :] = \ + (theta[time_index+1, :]-theta_tilde)/(gamma*delta_t) if (time_index + 1) % 100 == 0: print(f"Finished time step {time_index+1}") self.theta = theta - def solve_until_material_parameters_change( + def simulate_until_material_parameters_change( self, + delta_t: float, + number_of_time_steps: int, + gamma: float, initial_theta_field: npt.NDArray, initial_time_step: int ) -> int: - """Runs the simulation using the assembled c and k matrices as well - as the set excitation. - Calculates the temperature field and saves it in theta. + """Runs the time domain heat conduction simulation until the material + parameters of the set material model are changed sufficiently. + When the simulation is stopped due to a change of material parameters, + the last simulated time step is returned. Parameters: - theta_start: Field at time step start_index. - start_index: Sets the starting time step. + delta_t: Difference btweeen two time steps. + number_of_time_steps: Total number of time steps of the simulation. + gamma: Newmark integration parameter. + initial_theta_field: Initial temperature field at the first time + step. + initial_time_step: Sets the starting time step. This is needed, if + the simulation was interruped due to changed materials and is + now continued. The initial_thea_field must also be set + accordingly. """ - number_of_time_steps = self.simulation_data.number_of_time_steps - gamma = self.simulation_data.gamma - delta_t = self.simulation_data.delta_t - c = self.c k = self.k + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) + node_count = len(self.mesh_data.nodes) + # Init arrays if initial_time_step == 0: self.theta = np.zeros( - (k.shape[0], number_of_time_steps), + (number_of_time_steps, node_count), + dtype=np.float64 + ) + self.theta_dt = np.zeros( + (number_of_time_steps, node_count), dtype=np.float64 ) if initial_theta_field is not None: self.theta[:, initial_time_step] = initial_theta_field - # theta derived after t (du/dt) - self.theta_dt = np.zeros( - (k.shape[0], number_of_time_steps), - dtype=np.float64) - _, c, k = apply_dirichlet_bc( - np.zeros(c.shape), + _, c, k = mat_apply_dbcs( + None, c, k, - self.dirichlet_nodes + dirichlet_nodes ) k_star = (k+1/(gamma*delta_t)*c).tocsr() print("Starting heat conduction simulation") for time_index in range(initial_time_step, number_of_time_steps-1): - # Set dirichlet nodes for the load vector - # The load vector could already contain e.g. energy sources - # which are overwritten at the points where a dirichlet bc is set. - f = self.set_dirichlet_bc_load_vector( - self.f[:, time_index], - self.dirichlet_nodes, - self.dirichlet_values, - time_index+1 - ) + current_f = self.f[:, time_index+1] # Add a convection boundary condition if it set if self.enable_convection_bc: - f += self._calculate_convection_bc( - self.mesh_data.nodes, - self.convective_b_e, - self.theta[:, time_index], - self.convective_alpha, - self.convective_outer_temp, - self.mesh_data.element_order - ) + current_f += self._calculate_convection_bc( + self.mesh_data.nodes, + self.convective_b_e, + self.theta[time_index, :], + self.convective_alpha, + self.convective_outer_temp, + self.mesh_data.element_order + ) + + current_f[dirichlet_nodes] = dirichlet_values # Perform Newmark method # Predictor step @@ -512,23 +416,23 @@ def solve_until_material_parameters_change( ) # Solve for next time step - self.theta[:, time_index+1] = slin.spsolve( + self.theta[time_index+1, :] = slin.spsolve( k_star, ( - f - + 1/(gamma*delta_t)*c*theta_tilde + current_f + + 1/(gamma*delta_t)*c@theta_tilde ) ) # Perform corrector step - self.theta_dt[:, time_index+1] = \ - (self.theta[:, time_index+1]-theta_tilde)/(gamma*delta_t) + self.theta_dt[time_index+1, :] = \ + (self.theta[time_index+1, :]-theta_tilde)/(gamma*delta_t) if (time_index + 1) % 100 == 0: print(f"Finished time step {time_index+1}") # Check if material parameters would change temp_field_per_element = get_avg_temp_field_per_element( - self.theta[:, -1], + self.theta[-1, :], self.mesh_data.elements ) if self.material_manager.update_temperature( diff --git a/plutho/single_sim.py b/plutho/single_sim.py deleted file mode 100644 index efc0222..0000000 --- a/plutho/single_sim.py +++ /dev/null @@ -1,818 +0,0 @@ -"""Basic module to handle single simulations.""" -# TODO Remove this module and make an abstract class for the normal simulations -# Python standard libraries -import os -from enum import Enum -from typing import Union, List, Dict -import configparser -import json - -# Third party libraries -import numpy as np -import numpy.typing as npt - -# Local libraries -from .simulation.base import MeshData, SimulationData -from .simulation import ThermoSimTime, PiezoSimTime, ThermoPiezoSimTime, \ - PiezoSimFreq -from .materials import MaterialData, MaterialManager -from .mesh import Mesh -from .simulation.base import FieldType - -class SimulationException(Exception): - """Custom exception to simplify errors.""" - - - -class SingleSimulation: - """Base class to handle single simulations. Multiple different simulation - types can be chosen. - - Arguments: - simulation_name: Name of the current simulation. - simulation_directory: Directory where the results files are saved. - Combination of working_directory given to constructor and the - simulation name. - material_manager: MaterialManager object used to handle different - materials as well as temperature dependent material parameters. - mesh: Contains the mesh. - mesh_data: Contains the nodes and element arrays of the mesh. - dirichlet_nodes: Array of node indices at which a dirichlet boundary - shall be set. - dirichlet_values: Array of time dependent values which are set - at the corresponding dirichlet_node: [node_index, time_index]. - solver: Contains the simulation solver class. - charge_calculated: Is True if a charge is calculated in the simulation. - """ - # Basic settings - simulation_name: str - simulation_directory: str - - # Materials - material_manager: MaterialManager - - # Mesh - mesh: Mesh - mesh_data: MeshData - - # Dirichlet bc - boundary_conditions: List[Dict] # Used for serialization - dirichlet_nodes: List[int] - dirichlet_values: List[npt.NDArray] - - # Simulation - solver: Union[ThermoSimTime, PiezoSimFreq, - PiezoSimTime, ThermoPiezoSimTime] - charge_calculated: bool - mech_loss_calculated: bool - - def __init__( - self, - working_directory: str, - simulation_name: str, - mesh: Mesh - ): - simulation_directory = os.path.join(working_directory, simulation_name) - - self.simulation_directory = simulation_directory - self.simulation_name = simulation_name - self.mesh = mesh - - nodes, elements = mesh.get_mesh_nodes_and_elements() - self.mesh_data = MeshData(nodes, elements, mesh.element_order) - self.material_manager = MaterialManager(len(elements)) - self.charge_calculated = False - self.mech_loss_calculated = False - self.boundary_conditions = [] - - # Initialize dirichlet bc arrays - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - def add_material( - self, - material_name: str, - material_data: MaterialData, - physical_group_name: str - ): - """Adds a material to the simulation - - Parameters: - material_name: Name of the material. - material_data: Material data. - physical_group_name: Name of the physical group for which the - material shall be set. If this is None or "" the material will - be set for the whole model. - """ - if physical_group_name is None or physical_group_name == "": - element_indices = np.arange(len(self.mesh_data.elements)) - else: - element_indices = self.mesh.get_elements_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - self.material_manager.add_material( - material_name, - material_data, - physical_group_name, - element_indices - ) - - def add_dirichlet_bc( - self, - field_type: FieldType, - physical_group_name: str, - values: npt.NDArray - ): - """Adds a dirichlet boundary condition (dirichlet bc) to the - simulation. The given values are for each time step. Therefore each - value is applied to every node from the given physical_group. - - Parameters: - field_type: The type of field for which the bc shall be set. - physical_group_name: Name of the physical group for which the - boundary condition shall be set. - values: List of values per time step. Value of the bc for each time - step. The value for one time step is applied to every element - equally. Length: number_of_time_step. - """ - if isinstance(self.solver, ThermoSimTime): - if (field_type is FieldType.U_R or - field_type is FieldType.U_Z or - field_type is FieldType.PHI): - raise ValueError( - f"Unknown variable type {field_type} given for" - f"simulation type {type(self.solver)}") - elif (isinstance(self.solver, PiezoSimFreq) or - isinstance(self.solver, PiezoSimTime)): - if field_type is FieldType.THETA: - raise ValueError( - f"Unknown variable type {field_type} given for" - f"simulation type {type(self.solver)}") - - # Save boundary condition for serialization - self.boundary_conditions.append({ - "field_type": field_type.name, - "physical_group": physical_group_name, - "values": values.tolist() - }) - - # Apply the given values for all nodes from the given physical_group - node_indices = self.mesh.get_nodes_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - number_of_nodes = len(self.mesh_data.nodes) - for node_index in node_indices: - real_index = 0 - - # Depending on the variable type and the simulation types - # the corresponding field variable may be found at different - # positions of the solution vector - if field_type is FieldType.PHI: - real_index = 2*number_of_nodes+node_index - elif field_type is FieldType.THETA: - if isinstance(self.solver, ThermoSimTime): - real_index = node_index - else: - real_index = 3*number_of_nodes+node_index - elif field_type is FieldType.U_R: - real_index = 2*node_index - elif field_type is FieldType.U_Z: - real_index = 2*node_index+1 - else: - raise ValueError(f"Unknown variable type {field_type}") - - self.dirichlet_nodes.append(real_index) - self.dirichlet_values.append(values) - - def clear_dirichlet_bcs(self): - """Resets the dirichlet boundary conditions.""" - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - def setup_thermo_time_domain( - self, - delta_t: float, - number_of_time_steps: int, - gamma: float - ): - """Sets the simulation to a thermal simulation in the time domain. - - Parameters: - delta_t: Difference between the time steps. - number_of_time_steps: Total number of time steps. - gamma: Time integration parameter (Set to 0.5 for unconditionally - stable simulation). - """ - self.solver = ThermoSimTime( - self.mesh_data, - self.material_manager, - SimulationData( - delta_t, - number_of_time_steps, - gamma, - 0.0 - ) - ) - - def setup_piezo_time_domain( - self, - number_of_time_steps: int, - delta_t: float, - gamma: float, - beta: float - ): - """Sets a electro-mechanical (piezo) simulation in the time domain. - - Parameters: - number_of_time_steps: Number of time steps for which the simulation - is done - delta_t: Distance between two time steps. - gamma: Time integration parameter. - beta: Time integrationp parameter. - """ - self.solver = PiezoSimTime( - self.mesh_data, - self.material_manager, - SimulationData( - delta_t, - number_of_time_steps, - gamma, - beta - ) - ) - - def setup_piezo_freq_domain(self, frequencies: npt.NDArray): - """Sets a electro-mechanical (piezo) simulation in the frequency - domain - - Parameters: - frequencies: Frequencies which shall be simulated. - """ - self.solver = PiezoSimFreq( - self.mesh_data, - self.material_manager, - frequencies - ) - - def setup_thermo_piezo_time_domain( - self, - number_of_time_steps: int, - delta_t: float, - gamma: float, - beta: float - ): - """Sets a thermo-piezoelectrical simulation. - - Parameters: - number_of_time_steps: Number of time steps for which the simulation - is done. - delta_t: Distance between two time steps. - gamma: Time integration parameter. - beta: Time integrationp parameter. - """ - self.solver = ThermoPiezoSimTime( - self.mesh_data, - self.material_manager, - SimulationData( - delta_t, - number_of_time_steps, - gamma, - beta - ) - ) - - def simulate( - self, - material_starting_temperature: Union[ - npt.NDArray, - float, - None - ] = None, - temperature_dependent: bool = False, - initial_theta_field: Union[npt.NDArray, float, None] = None, - initial_theta_time_step: Union[npt.NDArray, None] = None, - electrode_elements: Union[npt.NDArray, None] = None, - electrode_normals: Union[npt.NDArray, None] = None, - calculate_mech_loss: bool = False - ): - # TODO Move parameters to **kwargs? - """Runs a simulation. Note that a simulation must be setup - beforehand. - - Parameters: - material_starting_temperature: Sets a starting temperature field. - initial_theta_field: Sets the starting field for the pure thermal - simulation. - initial_theta_time_step: Sets a starting time index for the pure - thermal simulation. - electrode_elements: List of element indices for which the - charge is calculated. - electrode_normals: List of normal vectors corresponding to each - electrode element (same index). - calculate_mech_loss: Set to true of the mech loss shall - be calculated. - """ - - # Check if solver is set - if self.solver is None: - raise SimulationException( - "Please setup a solver first using any " - "setup function" - ) - - # Check if there is a material defined for every point - if np.any(self.material_manager.element_index_to_material_data == -1): - raise SimulationException( - "Not every node has a material defined." - "Cannot start simulation." - ) - - # Check if mesh data is set - if self.mesh_data.nodes is None or self.mesh_data.elements is None: - raise SimulationException( - "Please setup mesh before starting simulation." - ) - - # Set material starting temperature if given - if material_starting_temperature is not None: - # Check if material temp is given as float - if isinstance(material_starting_temperature, float): - material_starting_temperature = ( - material_starting_temperature - * np.ones(len(self.mesh_data.elements)) - ) - - self.solver.material_manager.initialize_materials( - material_starting_temperature - ) - else: - self.solver.material_manager.initialize_materials() - - if electrode_elements is not None and electrode_normals is None: - raise SimulationException( - "If electrode elements are given, electrode normals " - "must be set too" - ) - - # Get boundary condition lists - self.solver.dirichlet_nodes = np.array(self.dirichlet_nodes) - self.solver.dirichlet_values = np.array(self.dirichlet_values) - - # Assemble - self.solver.assemble() - - # Check for the different solver types and run the simulation - if isinstance(self.solver, ThermoSimTime): - # Check if an initial theta field is given - if initial_theta_field is not None: - # Check if initial_theta_field is given as float - if isinstance(initial_theta_field, float): - initial_theta_field = ( - initial_theta_field - * np.ones(len(self.mesh_data.elements)) - ) - - # If the simulation is temperature dependent and there is an - # initial time step given the theta field will be solved until - # the temperature dependent material parameters change enough. - # If its enough is determined by the material manager. - if temperature_dependent and initial_theta_time_step is not None: - return self.solver.solve_until_material_parameters_change( - initial_theta_field, - initial_theta_time_step - ) - - self.solver.solve_time( - initial_theta_field - ) - elif isinstance(self.solver, PiezoSimFreq): - # Check if electrode_elements are given - if electrode_elements is not None: - self.charge_calculated = True - - # Check if mech loss shall be calculated - if calculate_mech_loss: - self.mech_loss_calculated = True - - self.solver.solve_frequency( - electrode_elements, - electrode_normals, - calculate_mech_loss - ) - elif isinstance(self.solver, PiezoSimTime): - # Check if electrode_elements are given - if electrode_elements is not None: - self.charge_calculated = True - - self.solver.solve_time( - electrode_elements, - electrode_normals - ) - elif isinstance(self.solver, ThermoPiezoSimTime): - # Check if electrode_elements are given - if electrode_elements is not None: - self.charge_calculated = True - - # Check if an initial theta field is given - if initial_theta_field is not None: - # Check if initial_theta_field is given as float - if isinstance(initial_theta_field, float): - initial_theta_field = ( - initial_theta_field - * np.ones(len(self.mesh_data.elements)) - ) - - # Mesh loss is always calculated - self.mech_loss_calculated = True - - self.solver.solve_time( - electrode_elements, - electrode_normals, - initial_theta_field - ) - else: - return - - def save_simulation_results(self, file_prefix: str = ""): - """Save the simulation results to the simulation folder. If a prefix - is given this is prepend to the name of the output file. - - Parameters: - file_prefix: String which is prepend to the output file names. - """ - if not os.path.exists(self.simulation_directory): - os.makedirs(self.simulation_directory) - - if file_prefix != "": - file_prefix += "_" - - if isinstance(self.solver, ThermoSimTime): - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}theta.npy" - ), - self.solver.theta - ) - elif isinstance(self.solver, PiezoSimFreq): - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}u.npy" - ), - self.solver.u - ) - if self.mech_loss_calculated: - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}mech_loss.npy" - ), - self.solver.mech_loss - ) - if self.charge_calculated: - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}q.npy" - ), - self.solver.q - ) - elif isinstance(self.solver, PiezoSimTime): - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}u.npy" - ), - self.solver.u - ) - if self.charge_calculated: - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}q.npy" - ), - self.solver.q - ) - elif isinstance(self.solver, ThermoPiezoSimTime): - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}u.npy" - ), - self.solver.u - ) - # Mech loss is always calculated in PiezoSimTherm - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}mech_loss.npy" - ), - self.solver.mech_loss - ) - if self.charge_calculated: - np.save( - os.path.join( - self.simulation_directory, - f"{file_prefix}q.npy" - ), - self.solver.q - ) - - def save_simulation_settings(self, prefix: str = ""): - """Save the simulation settings to the simulation folder. If a prefix - is given this is prepend to the name of the output file. - - Parameters: - file_prefix: String which is prepend to the output file names. - """ - if prefix != "": - prefix += "_" - - if not os.path.exists(self.simulation_directory): - os.makedirs(self.simulation_directory) - - settings = configparser.ConfigParser() - - if isinstance(self.solver, ThermoSimTime): - simulation_type = "ThermoTime" - elif isinstance(self.solver, PiezoSimTime): - simulation_type = "PiezoTime" - elif isinstance(self.solver, PiezoSimFreq): - simulation_type = "PiezoFreq" - elif isinstance(self.solver, ThermoPiezoSimTime): - simulation_type = "ThermoPiezoTime" - else: - raise ValueError( - f"Cannot save simulation type {type(self.solver)}" - ) - - # General simulation settings - general_settings = { - "name": self.simulation_name, - "mesh_file": self.mesh.mesh_file_path, - "simulation_type": simulation_type - } - settings["general"] = general_settings - - # Material settings and data - material_settings = {} - for material in self.material_manager.materials: - material_settings[material.material_name] = material.to_dict() - - if isinstance(self.solver, PiezoSimFreq): - # Save frequencies - np.savetxt( - os.path.join( - self.simulation_directory, - f"{prefix}frequencies.txt" - ), - self.solver.frequencies - ) - else: - # Simulation data - settings["simulation"] = self.solver.simulation_data.__dict__ - - # Boundary conditions - boundary_conditions = {} - for index, bc in enumerate(self.boundary_conditions): - boundary_conditions[index] = bc - - # Save simulation data to config file - with open( - os.path.join( - self.simulation_directory, - f"{prefix}{self.simulation_name}.cfg" - ), - "w", - encoding="UTF-8" - ) as fd: - settings.write(fd) - - # Save materials to json - with open( - os.path.join( - self.simulation_directory, - f"{prefix}materials.json" - ), - "w", - encoding="UTF-8" - ) as fd: - json.dump(material_settings, fd, indent=2) - - # Save boundary conditions to json - with open( - os.path.join( - self.simulation_directory, - f"{prefix}boundary_conditions.json" - ), - "w", - encoding="UTF-8" - ) as fd: - json.dump(boundary_conditions, fd, indent=2) - - @staticmethod - def load_simulation_settings( - simulation_folder: str, - prefix: str = "" - ): - """Loads the simulation settings from the given folder and return - a SingleSimulation object. - - Parameters: - simulation_folder: Path to a folder containing simulation files. - """ - # TODO Right now in a coupled sim the simulation settings - # are overwriding itself - # --> Add additional possible label to prevent this - - if prefix != "": - prefix += "_" - - simulation_name = os.path.basename(simulation_folder) - - # Check if folder has all necessary files - necessary_files = [ - os.path.join(simulation_folder, f"{prefix}{simulation_name}.cfg"), - os.path.join(simulation_folder, f"{prefix}materials.json"), - os.path.join( - simulation_folder, - f"{prefix}boundary_conditions.json" - ) - ] - for file_path in necessary_files: - if not os.path.isfile(file_path): - raise IOError(f"{file_path} does not exist.") - - settings = configparser.ConfigParser() - settings.read(necessary_files[0]) - - # Read general simulation settings - mesh_file = settings["general"]["mesh_file"] - simulation_type = settings["general"]["simulation_type"] - simulation = SingleSimulation( - simulation_folder, - "", - Mesh(mesh_file, True) - ) - # Workaround since simulation_folder and simulation_name are - # combined in the constructor, see empty string above for the - # constructor - simulation.simulation_name = simulation_name - - # Read simulation data - if simulation_type == "PiezoFreq": - # Check if file exists - frequencies_file = os.path.join( - simulation_folder, - "frequencies.txt" - ) - if os.path.isfile(frequencies_file): - frequencies = np.loadtxt(frequencies_file) - simulation.setup_piezo_freq_domain(frequencies) - else: - raise IOError( - "Frequencies file not found for frequency type sim" - ) - else: - simulation_data = SimulationData( - float(settings["simulation"]["delta_t"]), - int(settings["simulation"]["number_of_time_steps"]), - float(settings["simulation"]["gamma"]), - float(settings["simulation"]["beta"]), - ) - if simulation_type == "ThermoTime": - simulation.setup_thermo_time_domain( - simulation_data.delta_t, - simulation_data.number_of_time_steps, - simulation_data.gamma - ) - elif simulation_type == "PiezoTime": - simulation.setup_piezo_time_domain( - simulation_data.number_of_time_steps, - simulation_data.delta_t, - simulation_data.gamma, - simulation_data.beta - ) - elif simulation_type == "ThermoPiezoTime": - simulation.setup_thermo_piezo_time_domain( - simulation_data.number_of_time_steps, - simulation_data.delta_t, - simulation_data.gamma, - simulation_data.beta - ) - else: - raise IOError( - f"Cannot deserialize {simulation_type} simulation type" - ) - - # Read materials - with open(necessary_files[1], "r", encoding="UTF-8") as fd: - material_settings = json.load(fd) - for material_name in material_settings.keys(): - simulation.add_material( - material_name, - MaterialData( - **material_settings[material_name]["material_data"] - ), - material_settings[material_name]["physical_group_name"] - ) - simulation.material_manager.initialize_materials() - - # Read boundary conditions - with open(necessary_files[2], "r", encoding="UTF-8") as fd: - bc_settings = json.load(fd) - for index in bc_settings.keys(): - field_type = getattr( - FieldType, - bc_settings[index]["field_type"] - ) - physical_group = bc_settings[index]["physical_group"] - values = np.array(bc_settings[index]["values"]) - - simulation.add_dirichlet_bc( - field_type, - physical_group, - values - ) - - return simulation - - def load_simulation_results(self): - """Load the simulation results which are in the current simulation - folder to the current SingleSimulation object.""" - if isinstance(self.solver, ThermoSimTime): - theta_file = os.path.join( - self.simulation_directory, - "theta.npy" - ) - if os.path.isfile(theta_file): - self.solver.theta = np.load(theta_file) - else: - raise IOError("Couldn't find theta file.") - elif isinstance(self.solver, PiezoSimTime): - u_file = os.path.join( - self.simulation_directory, - "u.npy" - ) - q_file = os.path.join( - self.simulation_directory, - "q.npy" - ) - if os.path.isfile(u_file): - self.solver.u = np.load(u_file) - else: - raise IOError("Couldn't find u file.") - if os.path.isfile(q_file): - self.solver.q = np.load(q_file) - elif isinstance(self.solver, PiezoSimFreq): - u_file = os.path.join( - self.simulation_directory, - "u.npy" - ) - q_file = os.path.join( - self.simulation_directory, - "q.npy" - ) - mech_loss_file = os.path.join( - self.simulation_directory, - "mech_loss.npy" - ) - print(u_file) - if os.path.isfile(u_file): - self.solver.u = np.load(u_file) - else: - raise IOError("Couldn't find u file.") - if os.path.isfile(q_file): - self.solver.q = np.load(q_file) - if os.path.isfile(mech_loss_file): - self.solver.mech_loss = np.load(mech_loss_file) - elif isinstance(self.solver, ThermoPiezoSimTime): - u_file = os.path.join( - self.simulation_directory, - "u.npy" - ) - q_file = os.path.join( - self.simulation_directory, - "q.npy" - ) - mech_loss_file = os.path.join( - self.simulation_directory, - "mech_loss.npy" - ) - if os.path.isfile(u_file): - self.solver.u = np.load(u_file) - else: - raise IOError("Couldn't find u file.") - if os.path.isfile(q_file): - self.solver.q = np.load(q_file) - if os.path.isfile(mech_loss_file): - self.solver.mech_loss = np.load(mech_loss_file) - else: - raise ValueError( - "Cannot load simulation settings of simulation type", - type(self.solver) - ) diff --git a/scripts/basic_example.py b/scripts/basic_example.py index e600e61..402508e 100644 --- a/scripts/basic_example.py +++ b/scripts/basic_example.py @@ -87,21 +87,12 @@ def simulate_piezo_impedance(base_directory, show_results): # directory. mesh = load_mesh(os.path.join(base_directory, "disc_mesh.msh")) - # Create the simulation - sim = plutho.SingleSimulation( - # Folder IN which the simulation folder is created - working_directory=base_directory, - # Name of the simulation & sim folder name + # Create the frequency domain simulation + sim = plutho.PiezoFreq( simulation_name="pic255_impedance_example", - # Mesh which is used in the simulation mesh=mesh ) - # Frequencies for which the simulation is done - frequencies = np.linspace(0, 1e7, 1000)[1:] - # Set a frequency domain piezo simulation - sim.setup_piezo_freq_domain(frequencies) - # Add a material to the model and specify for which elements this element # is used. sim.add_material( @@ -113,6 +104,9 @@ def simulate_piezo_impedance(base_directory, show_results): physical_group_name="" ) + # Frequencies for which the simulation shall be done + frequencies = np.linspace(0, 1e7, 1000)[1:] + # Setup boundary conditions # First two boundary conditions for PHI (electrical potenzial) are added # those implement ground (0 V potential) on the bottom of the ceramic @@ -141,34 +135,21 @@ def simulate_piezo_impedance(base_directory, show_results): values=np.zeros(len(frequencies)) ) - # Since the charge shall be simulated in order to calculate the impedance - # it is necesarry to determine the elements for which a charge is - # calculated - electrode_elements = mesh.get_elements_by_physical_groups( - ["Electrode"] - )["Electrode"] - # Additionally the orientation of the electrode elements is needed to - # calculate the charge. Since they all lie on the uppder electrode, the - # (outer) normal vector only points in z direction. - electrode_normals = np.tile([0, 1], (len(electrode_elements), 1)) - # Now the simulation can be done - sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals - ) + sim.assemble() + sim.simulate(frequencies) - # The simulation settings and results are saved in the simulation folder - sim.save_simulation_settings() - sim.save_simulation_results() + # Charge can be calculated when given the name of the electrode (physical + # group in gmsh) + # Resulting charge is saved in sim.q + sim.calculate_charge("Electrode", is_complex=True) - # Calculate impedance - # The charge can be accessed using sim.solver - charge = sim.solver.q + sim.save_simulation_settings(base_directory) + sim.save_simulation_results(base_directory) # Since the simulation is already in the frequency domain the impedance # can be calculated directly - impedance = np.abs(1/(1j*2*np.pi*frequencies*charge)) + impedance = np.abs(1/(1j*2*np.pi*frequencies*sim.q)) if show_results: # Plot results @@ -193,27 +174,19 @@ def simulate_thermo_piezo(base_directory): """ mesh = load_mesh(os.path.join(base_directory, "disc_mesh.msh")) - sim = plutho.SingleSimulation( - base_directory, + sim = plutho.ThermoPiezoTime( "thermo_piezo_sim_20k", mesh ) + # Simulation parameters NUMBER_OF_TIME_STEPS = 20000 DELTA_T = 1e-8 - # Time integration parameters: Those values make sure that the simulation # is unconditionally stable. GAMMA = 0.5 BETA = 0.25 - sim.setup_thermo_piezo_time_domain( - NUMBER_OF_TIME_STEPS, - DELTA_T, - GAMMA, - BETA - ) - sim.add_material( "pic255", pic255, @@ -242,9 +215,18 @@ def simulate_thermo_piezo(base_directory): np.zeros(NUMBER_OF_TIME_STEPS) ) - sim.simulate(calculate_mech_loss=True) - sim.save_simulation_settings() - sim.save_simulation_results() + sim.assemble() + sim.simulate( + DELTA_T, + NUMBER_OF_TIME_STEPS, + GAMMA, + BETA + ) + + sim.calculate_charge("Electrode") + + sim.save_simulation_settings(base_directory) + sim.save_simulation_results(base_directory) def simulate_coupled_thermo_time(base_directory): @@ -278,11 +260,8 @@ def simulate_coupled_thermo_time(base_directory): 0 # Not needed in thermo simulation ) - coupled_sim = plutho.CoupledThermoPiezoThermoSim( - base_directory, + coupled_sim = plutho.CoupledThermoPiezoTime( "CoupledThermoPiezoelectricSim", - piezo_sim_data, - thermo_sim_data, mesh ) @@ -301,9 +280,14 @@ def simulate_coupled_thermo_time(base_directory): is_disc=True ) - coupled_sim.simulate() - coupled_sim.save_simulation_settings() - coupled_sim.save_simulation_results() + coupled_sim.assemble() + coupled_sim.simulate( + piezo_sim_data, + thermo_sim_data + ) + + coupled_sim.save_simulation_settings(base_directory) + coupled_sim.save_simulation_results(base_directory) def simulate_coupled_thermo_freq(base_directory): @@ -328,11 +312,8 @@ def simulate_coupled_thermo_freq(base_directory): 0 # Not needed in thermo simulation ) - coupled_sim = plutho.CoupledFreqPiezoTherm( - base_directory, + coupled_sim = plutho.CoupledPiezoThermoFreq( "CoupledThermopiezoelectricFreqSim", - FREQUENCY, - thermo_sim_data, mesh ) @@ -347,31 +328,30 @@ def simulate_coupled_thermo_freq(base_directory): is_disc=True ) + coupled_sim.assemble() coupled_sim.simulate( + FREQUENCY, + thermo_sim_data, material_starting_temperature=25, temperature_dependent=False ) - coupled_sim.save_simulation_settings() - coupled_sim.save_simulation_results() + + coupled_sim.save_simulation_settings(base_directory) + coupled_sim.save_simulation_results(base_directory) def simulate_thermo_time(base_directory): mesh = load_mesh(os.path.join(base_directory, "disc_mesh.msh")) - sim = plutho.SingleSimulation(base_directory, "ThermoTimeSim", mesh) + sim = plutho.ThermoTime("ThermoTimeSim", mesh) DELTA_T = 0.001 NUMBER_OF_TIME_STEPS = 1000 + GAMMA = 0.5 _, elements = mesh.get_mesh_nodes_and_elements() number_of_elements = len(elements) - sim.setup_thermo_time_domain( - DELTA_T, - NUMBER_OF_TIME_STEPS, - 0.5 - ) - sim.add_material( "pic255", pic255, @@ -379,7 +359,7 @@ def simulate_thermo_time(base_directory): ) # As an example set a constant volume heat source - sim.solver.set_constant_volume_heat_source( + sim.set_constant_volume_heat_source( np.ones(number_of_elements), NUMBER_OF_TIME_STEPS ) @@ -395,15 +375,22 @@ def simulate_thermo_time(base_directory): elements["RightBoundary"] ] ) - sim.solver.set_convection_bc( + + sim.set_convection_bc( convective_boundary_elements, 80, # Heat transfer coefficient 20 # Outer temperature ) - sim.simulate(initial_theta_field=25) - sim.save_simulation_settings() - sim.save_simulation_results() + sim.assemble() + sim.simulate( + DELTA_T, + NUMBER_OF_TIME_STEPS, + GAMMA + ) + + sim.save_simulation_settings(base_directory) + sim.save_simulation_results(base_directory) if __name__ == "__main__": @@ -415,7 +402,8 @@ def simulate_thermo_time(base_directory): if not os.path.isdir(CWD): os.makedirs(CWD) - simulate_piezo_impedance(CWD, False) - simulate_thermo_piezo(CWD) + # simulate_piezo_impedance(CWD, True) + # simulate_thermo_piezo(CWD) simulate_coupled_thermo_time(CWD) - simulate_coupled_thermo_freq(CWD) + # simulate_coupled_thermo_freq(CWD) + # simulate_thermo_time(CWD) diff --git a/tests/test_fields.py b/tests/test_fields.py index bf17e6d..09ea36c 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -9,7 +9,6 @@ # -------- Global variables -------- -NUMBER_OF_TIME_STEPS = 1000 ATOL = 1e-14 RTOL = 1e-4 @@ -51,27 +50,22 @@ def test_thermo_time(tmp_path): ) mesh = plutho.Mesh(mesh_path, element_order) - sim = plutho.SingleSimulation( - tmp_path, + sim = plutho.ThermoTime( "thermo_time_test", mesh ) # Simulation parameters INPUT_POWER_DENSITY = 1 - delta_t = 0.001 + NUMBER_OF_TIME_STEPS = 1000 + DELTA_T = 0.001 + GAMMA = 0.5 nodes, elements = mesh.get_mesh_nodes_and_elements() number_of_elements = len(elements) print(f"Nodes {nodes}") print(f"Elements {elements}") - sim.setup_thermo_time_domain( - delta_t, - NUMBER_OF_TIME_STEPS, - 0.5 - ) - sim.add_material( "pic255", pic255, @@ -79,16 +73,21 @@ def test_thermo_time(tmp_path): ) # Set constant volume loss density - sim.solver.set_constant_volume_heat_source( + sim.set_constant_volume_heat_source( INPUT_POWER_DENSITY*np.ones(number_of_elements), NUMBER_OF_TIME_STEPS ) - sim.simulate() + sim.assemble() + sim.simulate( + DELTA_T, + NUMBER_OF_TIME_STEPS, + GAMMA + ) # Calculated thermal stored energy temp_field_energy = plutho.calculate_stored_thermal_energy( - sim.solver.theta[:, -1], + sim.theta[-1, :], nodes, elements, element_order, @@ -97,8 +96,8 @@ def test_thermo_time(tmp_path): ) volume = np.sum( - plutho.simulation.base.calculate_volumes( - sim.solver.node_points, + plutho.calculate_volumes( + sim.node_points, element_order ) ) @@ -123,8 +122,7 @@ def test_piezo_time(tmp_path, test=True): ) mesh = plutho.Mesh(mesh_path, element_order) - sim = plutho.SingleSimulation( - tmp_path, + sim = plutho.PiezoTime( "piezo_time_test", mesh ) @@ -135,12 +133,10 @@ def test_piezo_time(tmp_path, test=True): "" ) - sim.setup_piezo_time_domain( - number_of_time_steps=NUMBER_OF_TIME_STEPS, - delta_t=1e-8, - gamma=0.5, - beta=0.25 - ) + DELTA_T = 1e-8 + NUMBER_OF_TIME_STEPS = 1000 + GAMMA = 0.5 + BETA = 0.25 # Set triangular excitation excitation = np.zeros(NUMBER_OF_TIME_STEPS) @@ -165,15 +161,16 @@ def test_piezo_time(tmp_path, test=True): np.zeros(NUMBER_OF_TIME_STEPS) ) - electrode_elements = mesh.get_elements_by_physical_groups( - ["Electrode"] - )["Electrode"] - electrode_normals = np.tile([0, 1], (len(electrode_elements), 1)) - + sim.assemble() sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals + DELTA_T, + NUMBER_OF_TIME_STEPS, + GAMMA, + BETA, ) + + sim.calculate_charge("Electrode") + test_folder_name = "piezo_time" if test: @@ -187,8 +184,8 @@ def test_piezo_time(tmp_path, test=True): test_q = np.load(os.path.join(test_folder, "q.npy")) test_u = np.load(os.path.join(test_folder, "u.npy")) - uut_q = sim.solver.q - uut_u = sim.solver.u[:, -1] + uut_q = sim.q + uut_u = sim.u[-1, :] # Compare arrays # For displacement just take last time step. TODO Is this sufficient? @@ -232,8 +229,7 @@ def test_piezo_freq(tmp_path, test=True): ) mesh = plutho.Mesh(mesh_path, element_order) - sim = plutho.SingleSimulation( - tmp_path, + sim = plutho.PiezoFreq( "piezo_freq_test", mesh ) @@ -244,8 +240,6 @@ def test_piezo_freq(tmp_path, test=True): "" ) - sim.setup_piezo_freq_domain(np.array([2e6])) - # Set boundary conditions sim.add_dirichlet_bc( plutho.FieldType.PHI, @@ -263,16 +257,13 @@ def test_piezo_freq(tmp_path, test=True): np.zeros(1) ) - electrode_elements = mesh.get_elements_by_physical_groups( - ["Electrode"] - )["Electrode"] - electrode_normals = np.array([[0, 1]] * len(electrode_elements)) - + sim.assemble() sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals + frequencies=np.array([2e6]) ) + sim.calculate_charge("Electrode", is_complex=True) + test_folder_name = "piezo_freq" if test: @@ -287,8 +278,8 @@ def test_piezo_freq(tmp_path, test=True): test_q = np.load(os.path.join(test_folder, "q.npy")) test_u = np.load(os.path.join(test_folder, "u.npy")) - uut_q = sim.solver.q - uut_u = sim.solver.u[:, -1] + uut_q = sim.q + uut_u = sim.u[-1, :] # Compare arrays # For displacement just take last time step. TODO Is this sufficient? @@ -334,8 +325,7 @@ def test_thermo_piezo_time(tmp_path, test=True): ) mesh = plutho.Mesh(mesh_path, element_order) - sim = plutho.SingleSimulation( - tmp_path, + sim = plutho.ThermoPiezoTime( "thermo_piezo_time", mesh ) @@ -351,13 +341,8 @@ def test_thermo_piezo_time(tmp_path, test=True): number_of_nodes = len(nodes) DELTA_T = 1e-8 NUMBER_OF_TIME_STEPS = 5000 - - sim.setup_thermo_piezo_time_domain( - number_of_time_steps=NUMBER_OF_TIME_STEPS, - delta_t=DELTA_T, - gamma=0.5, - beta=0.25 - ) + GAMMA = 0.5 + BETA = 0.25 # Set sinusoidal excitation AMPLITUDE = 1 @@ -388,30 +373,30 @@ def test_thermo_piezo_time(tmp_path, test=True): np.zeros(NUMBER_OF_TIME_STEPS) ) - electrode_elements = mesh.get_elements_by_physical_groups( - ["Electrode"] - )["Electrode"] - electrode_normals = np.array([[0, 1]] * len(electrode_elements)) - + sim.assemble() sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals + DELTA_T, + NUMBER_OF_TIME_STEPS, + GAMMA, + BETA ) + sim.calculate_charge("Electrode") + # Calculate input energy input_energy = plutho.calculate_electrical_input_energy( excitation, - sim.solver.q, + sim.q, DELTA_T ) # Calculate total loss energy - volumes = plutho.simulation.base.calculate_volumes( - sim.solver.node_points, element_order + volumes = plutho.calculate_volumes( + sim.node_points, element_order ) power = np.zeros(NUMBER_OF_TIME_STEPS) for time_step in range(NUMBER_OF_TIME_STEPS): - power[time_step] = np.dot(sim.solver.mech_loss[:, time_step], volumes) + power[time_step] = np.dot(sim.mech_loss[:, time_step], volumes) total_loss_energy = np.trapezoid( power, @@ -420,9 +405,9 @@ def test_thermo_piezo_time(tmp_path, test=True): ) # Calculate stored thermal energy - theta = sim.solver.u[3*number_of_nodes:] + theta = sim.u[3*number_of_nodes:] stored_thermal_energy = plutho.calculate_stored_thermal_energy( - theta[:, -1], + theta[-1, :], nodes, elements, element_order, @@ -452,9 +437,9 @@ def test_thermo_piezo_time(tmp_path, test=True): test_u = np.load(os.path.join(test_folder, "u.npy")) test_mech_loss = np.load(os.path.join(test_folder, "mech_loss.npy")) - uut_q = sim.solver.q - uut_u = sim.solver.u[:, -1] - uut_mech_loss = sim.solver.mech_loss[:, -1] + uut_q = sim.q + uut_u = sim.u[-1, :] + uut_mech_loss = sim.mech_loss[-1, :] # Compare arrays np.testing.assert_allclose( @@ -482,15 +467,15 @@ def test_thermo_piezo_time(tmp_path, test=True): # Save data np.save( os.path.join(tmp_path, test_folder_name, "q.npy"), - sim.solver.q + sim.q ) np.save( os.path.join(tmp_path, test_folder_name, "u.npy"), - sim.solver.u[:, -1] + sim.u[-1, :] ) np.save( os.path.join(tmp_path, test_folder_name, "mech_loss.npy"), - sim.solver.mech_loss[:, -1] + sim.mech_loss[-1, :] ) From 5b9fcd7e8bb00f72ff27f16d3f080527dba26f86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Mon, 27 Oct 2025 14:36:45 +0100 Subject: [PATCH 2/4] Readd nonlinear simulations scrips --- plutho/simulations/nonlinear/__init__.py | 3 + plutho/simulations/nonlinear/base.py | 519 +++++++++++++++++++++ plutho/simulations/nonlinear/piezo_hb.py | 503 ++++++++++++++++++++ plutho/simulations/nonlinear/piezo_time.py | 415 ++++++++++++++++ scripts/example_nonlinear.py | 8 +- 5 files changed, 1444 insertions(+), 4 deletions(-) create mode 100644 plutho/simulations/nonlinear/__init__.py create mode 100644 plutho/simulations/nonlinear/base.py create mode 100644 plutho/simulations/nonlinear/piezo_hb.py create mode 100644 plutho/simulations/nonlinear/piezo_time.py diff --git a/plutho/simulations/nonlinear/__init__.py b/plutho/simulations/nonlinear/__init__.py new file mode 100644 index 0000000..92d6ab1 --- /dev/null +++ b/plutho/simulations/nonlinear/__init__.py @@ -0,0 +1,3 @@ +from .base import assemble, Nonlinearity, NonlinearType +from .piezo_time import NLPiezoTime +from .piezo_hb import NLPiezoHB diff --git a/plutho/simulations/nonlinear/base.py b/plutho/simulations/nonlinear/base.py new file mode 100644 index 0000000..19eab78 --- /dev/null +++ b/plutho/simulations/nonlinear/base.py @@ -0,0 +1,519 @@ +"""Base module for various functions neede for the nonlinear simulation.""" + +# Python standard libraries +from typing import Tuple, Union + +# Third party libraries +import numpy as np +import numpy.typing as npt +from scipy import sparse + +# Local libraries +from ...mesh.mesh import MeshData +from ..integrals import local_to_global_coordinates, b_operator_global, \ + integral_m, integral_ku, integral_kuv, integral_kve, \ + quadratic_quadrature, gradient_local_shape_functions_2d +from ..helpers import create_node_points +from plutho.materials import MaterialManager +from ...enums import NonlinearType + +# +# -------- Functions -------- +# + +def assemble( + mesh_data: MeshData, + material_manager: MaterialManager +) -> Tuple[ + sparse.lil_array, + sparse.lil_array, + sparse.lil_array +]: + """Creates m, c, k FEM matrices based on the given mesh and material + data + + Parameters: + mesh_data: MeshData object containing the mesh nodes and elements. + material_manager: MaterialManager object containing the material data + + Returns: + The assembled mass matrix m, damping matrix c and stiffness matrix k. + """ + # TODO Assembly takes to long rework this algorithm + # Maybe the 2x2 matrix slicing is not very fast + # TODO Change to lil_array + nodes = mesh_data.nodes + elements = mesh_data.elements + element_order = mesh_data.element_order + node_points = create_node_points(nodes, elements, element_order) + + number_of_nodes = len(nodes) + mu = sparse.lil_matrix( + (2*number_of_nodes, 2*number_of_nodes), + dtype=np.float64 + ) + ku = sparse.lil_matrix( + (2*number_of_nodes, 2*number_of_nodes), + dtype=np.float64 + ) + kuv = sparse.lil_matrix( + (2*number_of_nodes, number_of_nodes), + dtype=np.float64 + ) + kv = sparse.lil_matrix( + (number_of_nodes, number_of_nodes), + dtype=np.float64 + ) + + for element_index, element in enumerate(elements): + current_node_points = node_points[element_index] + + # TODO Check if its necessary to calculate all integrals + # --> Dirichlet nodes could be leaved out? + # Multiply with jac_det because its integrated with respect to + # local coordinates. + # Mutiply with 2*pi because theta is integrated from 0 to 2*pi + mu_e = ( + material_manager.get_density(element_index) + * integral_m(current_node_points, element_order) + * 2 * np.pi + ) + ku_e = ( + integral_ku( + current_node_points, + material_manager.get_elasticity_matrix(element_index), + element_order + ) * 2 * np.pi + ) + kuv_e = ( + integral_kuv( + current_node_points, + material_manager.get_piezo_matrix(element_index), + element_order + ) * 2 * np.pi + ) + kve_e = ( + integral_kve( + current_node_points, + material_manager.get_permittivity_matrix( + element_index + ), + element_order + ) * 2 * np.pi + ) + + # Now assemble all element matrices + for local_p, global_p in enumerate(element): + for local_q, global_q in enumerate(element): + # Mu_e is returned as a 3x3 matrix (should be 6x6) + # because the values are the same. + # Only the diagonal elements have values. + mu[2*global_p, 2*global_q] += mu_e[local_p][local_q] + mu[2*global_p+1, 2*global_q+1] += mu_e[local_p][local_q] + + # Ku_e is a 6x6 matrix and 2x2 matrices are sliced out + # of it. + ku[2*global_p:2*global_p+2, 2*global_q:2*global_q+2] += \ + ku_e[2*local_p:2*local_p+2, 2*local_q:2*local_q+2] + + # KuV_e is a 6x3 matrix and 2x1 vectors are sliced out. + # [:, None] converts to a column vector. + kuv[2*global_p:2*global_p+2, global_q] += \ + kuv_e[2*local_p:2*local_p+2, local_q][:, None] + + # KVe_e is a 3x3 matrix and therefore the values can + # be set directly. + kv[global_p, global_q] += kve_e[local_p, local_q] + + # Calculate damping matrix + # Currently in the simulations only one material is used + # TODO Add algorithm to calculate cu for every element on its own + # in order to use multiple materials + cu = ( + material_manager.get_alpha_m(0) * mu + + material_manager.get_alpha_k(0) * ku + ) + + # Calculate block matrices + zeros1x1 = np.zeros((number_of_nodes, number_of_nodes)) + zeros2x1 = np.zeros((2*number_of_nodes, number_of_nodes)) + zeros1x2 = np.zeros((number_of_nodes, 2*number_of_nodes)) + + m = sparse.block_array([ + [mu, zeros2x1], + [zeros1x2, zeros1x1], + ]).tolil() + c = sparse.block_array([ + [cu, zeros2x1], + [zeros1x2, zeros1x1], + ]).tolil() + k = sparse.block_array([ + [ku, kuv], + [kuv.T, -1*kv] + ]).tolil() + + return m, c, k + + +def integral_u_nonlinear_quadratic( + node_points: npt.NDArray, + nonlinear_elasticity_tensor: npt.NDArray, + element_order: int +): + """Calculates the nonlinear quadratic integral over the given node points. + + Parameters: + node_points: Points of the current element. + nonlinear_elasticity_tensor: Tensor for a quadratic nonlinear material + model. + element_order: Order of the element. + """ + def inner(s: float, t: float): + dn = gradient_local_shape_functions_2d(s, t, element_order) + jacobian = np.dot(node_points, dn.T) + jacobian_inverted_t = np.linalg.inv(jacobian).T + jacobian_det = np.linalg.det(jacobian) + + b_op = b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + B_TC = np.einsum("im,mjk->ijk", b_op.T, nonlinear_elasticity_tensor) + B_TCB = np.einsum("ink,nj->ijk", B_TC, b_op) + B_TCBB = np.einsum("ijp,pk->ijk", B_TCB, b_op) + + return B_TCBB * r * jacobian_det + + return quadratic_quadrature(inner, element_order) + +# +# -------- Classes -------- +# +class Nonlinearity: + """Class for handling the different nonlinearity types. Adds an additional + layer of abstraction so the solver can be the same for all nonlinear + types. + + Attributes: + mesh_data: Mesh data object. + node_points: List of node points for each element. + nonlinear_type: The type of the nonlinearity. + zeta: Nonlinear data parameter. + custom_data: Nonlinear data parameter. + """ + mesh_data: Union[MeshData, None] + node_points: Union[npt.NDArray, None] + + # Nonlinear data + nonlinear_type: Union[NonlinearType, None] + zeta: Union[float, None] + custom_data: Union[npt.NDArray, None] + + def __init__(self): + self.nonlinear_type = None + self.mesh_data = None + self.node_points = None + + def set_mesh_data(self, mesh_data: MeshData, node_points: npt.NDArray): + """Sets the mesh data. + + Parameters: + mesh_data: MeshData to set. + node_points: List of node points for each element. + """ + self.mesh_data = mesh_data + self.node_points = node_points + + def set_quadratic_rayleigh(self, zeta: float): + """Sets a quadratic rayleigh nonlinearity which only needs one + parameter. + + Parameters: + zeta: Nonlinearity parameter. + """ + self.nonlinear_type = NonlinearType.QuadraticRayleigh + self.zeta = zeta + + def set_quadratic_custom(self, nonlinear_data: npt.NDArray): + """Sets a quadratic nonlinearity with a custom material tensor. + + Parameters: + nonlinear_data: Nonlinear material tensor. + """ + self.nonlinear_type = NonlinearType.QuadraticCustom + + # Reduce to axisymmetric 4x4x4 matrix + # TODO Update this for order 3 or make it n dimensional + voigt_map = [0, 1, 3, 2] + nm_reduced = np.zeros(shape=(4, 4, 4)) + for i_new, i_old in enumerate(voigt_map): + for j_new, j_old in enumerate(voigt_map): + for k_new, k_old in enumerate(voigt_map): + nm_reduced[ + i_new, + j_new, + k_new + ] = nonlinear_data[i_old, j_old, k_old] + + self.custom_data = nm_reduced + + def set_cubic_rayleigh(self, zeta: float): + """Sets a cubic rayleigh nonlinearity which only needs one parameter. + + Parameters: + zeta: Nonlinearity parameter. + """ + self.nonlinear_type = NonlinearType.CubicRayleigh + self.zeta = zeta + + def set_cubic_custom(self, nonlinear_data: npt.NDArray): + """Sets a cubic nonlinearity with a custom nonlinear material tensor + + Parameters: + nonlinear_data: Nonlinear material tensor. + """ + self.nonlinear_type = NonlinearType.CubicCustom + + # Reduce to axisymmetric 4x4x4 matrix + # TODO Update this for order 3 or make it n dimensional + voigt_map = [0, 1, 3, 2] + nm_reduced = np.zeros(shape=(4, 4, 4)) + for i_new, i_old in enumerate(voigt_map): + for j_new, j_old in enumerate(voigt_map): + for k_new, k_old in enumerate(voigt_map): + for l_new, l_old in enumerate(voigt_map): + nm_reduced[ + i_new, + j_new, + k_new, + l_new + ] = nonlinear_data[i_old, j_old, k_old, l_old] + + self.custom_data = nm_reduced + + def evaluate_force_vector( + self, + u: npt.NDArray, + m: sparse.csc_array, + c: sparse.csc_array, + k: sparse.csc_array + ) -> npt.NDArray: + """Evaluates the nonlinear force vector based on the previously set + nonlinearity and the given displacement field u and the FEM matrices. + + + Parameters: + u: Current displacement field. + m: FEM mass matrix. + c: FEM damping matrix. + k: FEM stiffness matrix. + + Returns: + Vector of nonlinear forces. + """ + if self.nonlinear_type is None: + raise ValueError("Cannot evaluate force vector, since no \ + nonlinear type is set") + + match self.nonlinear_type: + case NonlinearType.QuadraticRayleigh: + return self.ln@(u**2) + case NonlinearType.QuadraticCustom: + # Use jit compiler --> wimp + """ + result = np.zeros(3*number_of_nodes) + + for index in range(3*number_of_nodes): + L_k = L[index*3*number_of_nodes:(index+1)*3*number_of_nodes, :] + + result[index] = u @ (L_k @ u) + + return result + """ + raise NotImplementedError() + case NonlinearType.CubicRayleigh: + return self.ln@(u**3) + case NonlinearType.CubicCustom: + raise NotImplementedError() + + def evaluate_jacobian( + self, + u: npt.NDArray, + m: sparse.csc_array, + c: sparse.csc_array, + k: sparse.csc_array + ): + """Evaluates the jacobian of the nonlinear force vector based on the + current displacement u and the FEM matrices. + + Parameters: + u: Current displacement vector. + m: FEM mass matrix. + c: FEM damping matrix. + k: FEM stiffness matrix. + """ + if self.nonlinear_type is None: + raise ValueError("Cannot evaluate jacobian, since no \ + nonlinear type is set") + + match self.nonlinear_type: + case NonlinearType.QuadraticRayleigh: + return 2*self.ln@sparse.diags_array(u) + case NonlinearType.QuadraticCustom: + """ + k_tangent = sparse.lil_matrix(( + 3*number_of_nodes, 3*number_of_nodes + )) + + k_tangent += k + + m = np.arange(3*number_of_nodes) + + # TODO This calculation is very slow + for i in range(3*number_of_nodes): + l_k = ln[3*number_of_nodes*i:3*number_of_nodes*(i+1), :] + l_i = ln[3*number_of_nodes*m+i, :] + + k_tangent += (l_k*u[i]).T + k_tangent += (l_i*u[i]).T + """ + raise NotImplementedError() + case NonlinearType.CubicRayleigh: + return 3*self.ln@sparse.diags_array(u**2) + case NonlinearType.CubicCustom: + raise NotImplementedError() + + def apply_dirichlet_bc(self, dirichlet_nodes: npt.NDArray): + """Applies the dirichlet boundary conditions on the nonlinear matrix.""" + if self.nonlinear_type is None: + raise ValueError("Cannot apply dirichlet boundary \ + conditions, since no nonlinear type is set") + + # TODO Is a different handling for nonlinear types necessary? + match self.nonlinear_type: + case NonlinearType.QuadraticRayleigh | \ + NonlinearType.CubicRayleigh: + self.ln[dirichlet_nodes, :] = 0 + case NonlinearType.QuadraticCustom: + raise NotImplementedError() + case NonlinearType.CubicCustom: + raise NotImplementedError() + + def assemble( + self, + m: sparse.lil_array, + c: sparse.lil_array, + k: sparse.lil_array + ): + """Assembles the nonlinar FEM matrix. It is saved in self.ln. + + Parameters: + m: FEM mass matrix. + c: FEM damping matrix. + k: FEM stiffness matrix. + """ + if self.nonlinear_type is None: + raise ValueError("Cannot assemble matrices, since no \ + nonlinear type is set") + + match self.nonlinear_type: + case NonlinearType.QuadraticRayleigh | \ + NonlinearType.CubicRayleigh: + self.ln = k*self.zeta + case NonlinearType.QuadraticCustom: + self._assemble_quadratic_custom() + case NonlinearType.CubicCustom: + self._assemble_cubic_custom() + + def _assemble_quadratic_custom(self): + """Local function. Assembles the nonlinear FEM matrix for a quadratic + custom nonlinearity type. + """ + nodes = self.mesh_data.nodes + elements = self.mesh_data.elements + element_order = self.mesh_data.element_order + number_of_nodes = len(nodes) + node_points = create_node_points( + nodes, + elements, + element_order + ) + + lu = sparse.lil_array( + (4*number_of_nodes**2, 2*number_of_nodes), + dtype=np.float64 + ) + + for element_index, element in enumerate(elements): + current_node_points = node_points[element_index] + + lu_e = integral_u_nonlinear_quadratic( + current_node_points, + self.custom_data, + element_order + ) * 2 * np.pi + + for local_i, global_i in enumerate(element): + for local_j, global_j in enumerate(element): + for local_k, global_k in enumerate(element): + # lu_e is 6x6x6 but has to be assembled into an + # 4*N**2x2*N sparse matrix + # The k-th plane is therefore append using the i-th + # index + # Since lu_e is 6x6x6, 2x2x2 slices are taken out and + # used for assembling + k_0 = lu_e[ + 2*local_i:2*(local_i+1), + 2*local_j:2*(local_j+1), + 2*local_k + ] + k_1 = lu_e[ + 2*local_i:2*(local_i+1), + 2*local_j:2*(local_j+1), + 2*local_k+1 + ] + + lu[ + ( + 2*global_i + + 3*number_of_nodes*global_k + ):( + 2*(global_i+1) + + 3*number_of_nodes*global_k + ), + 2*global_j:2*(global_j+1), + ] += k_0 + lu[ + ( + 2*global_i + + 3*number_of_nodes*global_k+1 + ):( + 2*(global_i+1) + + 3*number_of_nodes*global_k+1 + ), + 2*global_j:2*(global_j+1), + ] += k_1 + + ln = sparse.lil_array( + (9*number_of_nodes**2, 3*number_of_nodes), + dtype=np.float64 + ) + + for k in range(3*number_of_nodes): + ln[ + k*3*number_of_nodes:k*3*number_of_nodes+2*number_of_nodes, + :2*number_of_nodes + ] = lu[ + k*2*number_of_nodes:(k+1)*2*number_of_nodes, : + ] + + self.ln = ln + + + def _assemble_cubic_custom(self): + pass diff --git a/plutho/simulations/nonlinear/piezo_hb.py b/plutho/simulations/nonlinear/piezo_hb.py new file mode 100644 index 0000000..963a904 --- /dev/null +++ b/plutho/simulations/nonlinear/piezo_hb.py @@ -0,0 +1,503 @@ +"""Module for the simulation of nonlinaer piezoelectric systems using the +harmonic balancing method.""" + +# Python standard libraries +from typing import Union, Tuple, Callable + +# Third party libraries +import numpy as np +import numpy.typing as npt +import scipy +from scipy import sparse +from scipy.sparse import linalg + +# Local libraries +from .base import assemble, NonlinearType +from ...mesh.mesh import MeshData +from plutho.materials import MaterialManager + + +class NLPiezoHB: + """Implementes a nonlinear FEM harmonic balancing simulation. + """ + # Simulation parameters + mesh_data: MeshData + material_manager: MaterialManager + + # Boundary conditions + dirichlet_nodes: npt.NDArray + dirichlet_values: npt.NDArray + + # FEM matrices + k: sparse.lil_array + c: sparse.lil_array + m: sparse.lil_array + lu: sparse.lil_array + + # Resulting fields + u: npt.NDArray + + def __init__( + self, + mesh_data: MeshData, + material_manager: MaterialManager, + harmonic_order: int + ): + self.mesh_data = mesh_data + self.material_manager = material_manager + self.harmonic_order = harmonic_order + self.dirichlet_nodes = np.array([]) + self.dirichlet_valuse = np.array([]) + + def assemble( + self, + nonlinear_order: int, + nonlinear_type: NonlinearType, + **kwargs + ): + """Redirect to general nonlinear assembly function. + + Parameters: + nonlinear_order: Order of the nonlinearity. + nonlinear_type: Type of the nonlinear material. + **kwargs: Parameters for the nonlinear material. + """ + # Get the default FEM matrices + m, c, k = assemble( + self.mesh_data, + self.material_manager, + ) + self.nonlinear_order = nonlinear_order + + self.m = m + self.c = c + self.k = k + + raise NotImplementedError("Nonlinear assembly not implemented for HB") + + def solve_linear(self): + """Solves the linear problem Ku=f (ln is not used). + """ + # Apply boundary conditions + _, _, k = NLPiezoHB.apply_dirichlet_bc( + None, + None, + self.k.copy(), + self.dirichlet_nodes + ) + + f = self.get_load_vector( + self.dirichlet_nodes, + self.dirichlet_values + ) + + # Calculate u using phi as load vector + self.u = linalg.spsolve( + k, + f + ).todense() + + def solve_newton( + self, + angular_frequency: float, + tolerance: float = 1e-10, + max_iter: int = 100, + u_start: Union[npt.NDArray, None] = None, + alpha: float = 1, + load_factor: float = 1 + ): + """Solves the system of nonlinear equations using the Newton-Raphson + method. Saves the result in self.u + + Parameters: + angular_frequency: Angualr frequency at which the simulation is + done. + tolerance: Upper bound for the residuum. Iteration stops when the + residuum is lower than tolerance. + max_iter: Maximum number of iteration when the tolerance is not + met. Then the best solution (lowest residuum) is returned + u_start: Initial guess for the nonlinear problem. If it is None + a linear solution is calculated and used as a guess for + nonlinear u. + alpha: Damping parameter for updating the guess for u. Can be used + to improve convergence. Choose a value between 0 and 1. + load_factor: Multiplied with the load vector. Used to apply less + than than the full excitation. Could improve conversion of the + algorithm. Choose a value between 0 and 1. + """ + number_of_nodes = len(self.mesh_data.nodes) + + # Get FEM matrices + fem_m = self.m.copy() + fem_c = self.c.copy() + fem_k = self.k.copy() + fem_lu = self.lu.copy() + + # Create the harmonic balancing matrices from the FEM matrices + hb_m, hb_c, hb_k = self.get_harmonic_balancing_matrices( + fem_m, + fem_c, + fem_k, + angular_frequency + ) + + # Apply boundary conditions + hb_m, hb_c, hb_k = NLPiezoHB.apply_dirichlet_bc( + hb_m, + hb_c, + hb_k, + self.dirichlet_nodes + ) + + f = self.get_load_vector( + self.dirichlet_nodes, + self.dirichlet_values + ) + + # Set start value + if u_start is None: + # If no start value is given calculate one using a linear + # simulation + current_u = linalg.spsolve( + ( + hb_m + + hb_c + + hb_k + ), + load_factor*f + ) + else: + current_u = u_start + + # Residual of the Newton algorithm + def residual(u, nonlinear_force): + return ( + hb_m@u + + hb_c@u + + hb_k@u + + nonlinear_force + - load_factor*f + ) + + # Nonlinear forces based on starting field guess + nonlinear_force = self.calculate_nonlinear_force( + current_u, + fem_lu, + number_of_nodes + ) + + # best_field tracks the field with the best norm + # best norm tracks the value of the norm itself + # Calculate the first norms based on the starting field guess + best_field = current_u.copy() + best_norm = np.linalg.norm(residual(best_field, nonlinear_force)) + + # Check if the initial value already suffices the tolerance condition + #if best_norm < tolerance: + # print(f"Initial value is already best {best_norm}") + # self.u = best_field + # return + + # Run Newton method + for iteration_count in range(max_iter): + if iteration_count > 0: + # Does not need to be updated at step 0 + nonlinear_force = self.calculate_nonlinear_force( + current_u, + fem_lu, + number_of_nodes + ) + + # Calculate tangential stiffness matrix + k_tangent = self.calculate_tangent_matrix_analytical( + current_u, + hb_m, + hb_c, + hb_k + ) + + # Solve for displacement increment + delta_u = linalg.spsolve( + k_tangent, + residual(current_u, nonlinear_force) + ) + + # Update step + next_u = current_u - alpha * delta_u + + # Check for convergence + norm = scipy.linalg.norm(residual(next_u, nonlinear_force)) + if norm < tolerance: + print( + f"Newton found solution after {iteration_count+1} " + f"steps with residual {norm}" + ) + self.u = next_u + return + elif norm < best_norm: + best_norm = norm + best_field = next_u.copy() + + if iteration_count % 100 == 0 and iteration_count > 0: + print("Iteration:", iteration_count) + + # Update for next iteration + current_u = next_u + + print("Simulation did not converge") + print("Error from best iteration:", best_norm) + + self.u = best_field + + def calculate_nonlinear_force( + self, + u: npt.NDArray, + fem_lu: sparse.lil_array, + number_of_nodes: int + ) -> npt.NDArray: + """Calculates the forces due to the nonlinearity. + + Parameters: + u: Current solution vector. + fem_lu: FEM nonlinearity matrix. + number_of_nodes: Number of nodes of the mesh. + + Returns: + Vector of nonlinear forces. + """ + # TODO Currently only for third order nonlinearity + + u_c = u[::2] + u_s = u[1::2] + + lu_c = fem_lu.dot( + 3/4*u_c**3+3/4*np.multiply(u_c, u_s**2) + ) + lu_s = fem_lu.dot( + 3/4*u_s**3+3/4*np.multiply(u_s, u_c**2) + ) + + nonlinear_force = np.zeros(2*3*number_of_nodes) + nonlinear_force[::2] = lu_c + nonlinear_force[1::2] = lu_s + + """ + nonlinear_force = np.zeros(2*3*number_of_nodes) + for i in range(3*number_of_nodes): + a = u[2*i] + b = u[2*i+1] + nonlinear_force[2*i:2*i+2] = fem_lu[i,i]*np.array([ + 3/4*(a**3+a*b**2), + 3/4*(a**2*b+b**3) + ]) + """ + + # Add dirichlet boundary conditions + nonlinear_force[self.dirichlet_nodes] = 0 + + return nonlinear_force + + def get_harmonic_balancing_matrices( + self, + fem_m: sparse.lil_array, + fem_c: sparse.lil_array, + fem_k: sparse.lil_array, + angular_frequency: float + ) -> Tuple[sparse.lil_array, sparse.lil_array, sparse.lil_array]: + """Calculates the harmonic balancing matrices from the FEM matrices. + The matrices are getting bigger by a factor of (2*harmonic_order)**2 + + Parameters: + fem_m: FEM mass matrix. + fem_c: FEM damping matrix. + fem_k: FEM stiffness matrix. + + Returns: + Tuple of HB mass matrix, HB damping matrix and HB stiffness + matrix. + """ + harmonic_order = self.harmonic_order + + # Create harmonic balance derivative matrices + nabla_m = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + nabla_c = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + nabla_k = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + for i in range(harmonic_order): + k = i + 1 + nabla_m[2*i:2*i+2, 2*i:2*i+2] = np.array([ + [-(angular_frequency**2*k**2), 0], + [0, -(angular_frequency**2*k**2)] + ]) + nabla_c[2*i:2*i+2, 2*i:2*i+2] = np.array([ + [0, angular_frequency*k], + [-(angular_frequency*k), 0] + ]) + nabla_k[2*i:2*i+2, 2*i:2*i+2] = np.array([ + [1, 0], + [0, 1] + ]) + + # Apply HB derivative matrices on FEM matrices + fem_m = sparse.kron(fem_m, nabla_m, format='lil').tolil() + fem_c = sparse.kron(fem_c, nabla_c, format='lil').tolil() + fem_k = sparse.kron(fem_k, nabla_k, format='lil').tolil() + + return fem_m, fem_c, fem_k + + def get_load_vector( + self, + nodes: npt.NDArray, + values: npt.NDArray, + ) -> npt.NDArray: + """Calculates the load vector (right hand side) vector for the + simulation. + + Parameters: + nodes: Nodes at which the dirichlet value shall be set. + values: Dirichlet value which is set at the corresponding node. + + Returns: + Right hand side vector for the simulation. + """ + number_of_nodes = len(self.mesh_data.nodes) + + # Can be initialized to 0 because external load and volume + # charge density is 0. + f = np.zeros( + 3*number_of_nodes*2*self.harmonic_order, + dtype=np.float64 + ) + + for node, value in zip(nodes, values): + f[node] = value + + return f + + def calculate_tangent_matrix_analytical( + self, + u: npt.NDArray, + hb_m: sparse.lil_array, + hb_c: sparse.lil_array, + hb_k: sparse.lil_array + ) -> sparse.csc_array: + """Calculates the tangent matrix using a priori analytical + calculations. + + Parameters: + u: Solution vector containing mechanical and electrical field + values. + hb_m: Harmonic balancing mass matrix. + hb_c: Harmonic balancing damping matrix. + hb_k: Harmonic balancing stiffness matrix. + + Returns: + The tangent stiffness matrix. + """ + # TODO Only valid for harmonic_order = 1 and nonlinearity_order = 3 + fem_lu = self.lu + lu_diag = fem_lu.diagonal() + number_of_nodes = len(self.mesh_data.nodes) + size = number_of_nodes*3*self.harmonic_order*2 + + center_diag = np.zeros(size) + for i in range(len(lu_diag)): + diag_value_first = 3/4*lu_diag[i]*(3*u[2*i]**2+u[2*i+1]**2) + diag_value_second = 3/4*lu_diag[i]*(u[2*i]**2+3*u[2*i+1]**2) + + center_diag[2*i] = diag_value_first + center_diag[2*i+1] = diag_value_second + + left_diag = np.zeros(size-1) + left_diag[::2] = np.multiply(u[1::2], lu_diag) + + right_diag = np.zeros(size-1) + right_diag[::2] = np.multiply(u[:-1:2], lu_diag) + + return hb_m + hb_c + hb_k + sparse.diags( + diagonals=[ + left_diag, + center_diag, + right_diag, + ], + offsets=[-1, 0, 1], + format="csc" + ) + + + # -------- Static functions -------- + + @staticmethod + def calculate_tangent_matrix_numerical( + u: npt.NDArray, + nonlinear_force: npt.NDArray, + residual: Callable + ) -> sparse.csc_array: + """Calculates the tanget stiffness matrix numerically using a first + order finite difference approximation. + Important note: Calculation is very slow. + + Parameters: + u: Current solution vector. + nonlinear_force: Vector of nonlinear forces. + residual: Residual function. + + Returns: + Tangential stiffness matrix in sparse csc format. + """ + h = 1e-14 + k_tangent = sparse.lil_matrix( + (len(u), len(u)), + dtype=np.float64 + ) + for i in range(len(u)): + e_m = np.zeros(len(u)) + e_m[i] = h + k_tangent[:, i] = ( + residual(u+e_m, nonlinear_force) + - residual(u, nonlinear_force) + )/h + + return k_tangent.tocsc() + + + @staticmethod + def apply_dirichlet_bc( + m: sparse.lil_array, + c: sparse.lil_array, + k: sparse.lil_array, + nodes: npt.NDArray + ) -> Tuple[ + sparse.csc_array, + sparse.csc_array, + sparse.csc_array + ]: + """Applies dirichlet boundary conditions to the given matrix based on + the given nodes. Essentially the rows of the m c and k matrices at the + node indices are set to 0. The element of the k matrix at position + (node_index, node_index) is set to 1. + + Parameters: + m: Mass matrix. + c: Damping matrix. + k: Stiffness matrix. + nodes: Nodes at which the dirichlet boundary conditions shall be + applied. + + Returns: + Modified mass, damping and stiffness matrix. + """ + # Set rows of matrices to 0 and diagonal of K to 1 (at node points) + # Matrices for u_r component + for node in nodes: + # Set rows to 0 + if m is not None: + m[node, :] = 0 + if c is not None: + c[node, :] = 0 + if k is not None: + k[node, :] = 0 + + # Set diagonal values to 1 + k[node, node] = 1 + + return m.tocsc(), c.tocsc(), k.tocsc() diff --git a/plutho/simulations/nonlinear/piezo_time.py b/plutho/simulations/nonlinear/piezo_time.py new file mode 100644 index 0000000..26c0e52 --- /dev/null +++ b/plutho/simulations/nonlinear/piezo_time.py @@ -0,0 +1,415 @@ +"""Module for the simulation of nonlinear piezoelectric systems""" + +# Third party libraries +from typing import Tuple, Union +import numpy as np +import numpy.typing as npt +import scipy +from scipy import sparse +from scipy.sparse import linalg + +# Local libraries +from .base import assemble, Nonlinearity +from ...mesh.mesh import MeshData +from ...materials import MaterialData +from ...enums import FieldType +from ..helpers import calculate_charge, create_node_points +from plutho.materials import MaterialManager +from plutho.mesh.mesh import Mesh + + +class NLPiezoTime: + + def __init__( + self, + simulation_name: str, + mesh: Mesh, + nonlinearity: Nonlinearity + ): + self.simulation_name = simulation_name + + # Setup simulation data + self.boundary_conditions = [] + self.dirichlet_nodes = [] + self.dirichlet_values = [] + + # Setup mesh and material data + nodes, elements = mesh.get_mesh_nodes_and_elements() + self.mesh = mesh + self.mesh_data = MeshData(nodes, elements, mesh.element_order) + self.node_points = create_node_points( + nodes, + elements, + mesh.element_order + ) + + self.material_manager = MaterialManager(len(elements)) + self.nonlinearity = nonlinearity + self.nonlinearity.set_mesh_data(self.mesh_data, self.node_points) + + def add_material( + self, + material_name: str, + material_data: MaterialData, + physical_group_name: str + ): + """Adds a material to the simulation + + Parameters: + material_name: Name of the material. + material_data: Material data. + physical_group_name: Name of the physical group for which the + material shall be set. If this is None or "" the material will + be set for the whole model. + """ + if physical_group_name is None or physical_group_name == "": + element_indices = np.arange(len(self.mesh_data.elements)) + else: + element_indices = self.mesh.get_elements_by_physical_groups( + [physical_group_name] + )[physical_group_name] + + self.material_manager.add_material( + material_name, + material_data, + physical_group_name, + element_indices + ) + + def add_dirichlet_bc( + self, + field_type: FieldType, + physical_group_name: str, + values: npt.NDArray + ): + """Adds a dirichlet boundary condition (dirichlet bc) to the + simulation. The given values are for each time step. Therefore each + value is applied to every node from the given physical_group. + + Parameters: + field_type: The type of field for which the bc shall be set. + physical_group_name: Name of the physical group for which the + boundary condition shall be set. + values: List of values per time step. Value of the bc for each time + step. The value for one time step is applied to every element + equally. Length: number_of_time_step. + """ + # Save boundary condition for serialization + self.boundary_conditions.append({ + "field_type": field_type.name, + "physical_group": physical_group_name, + "values": values.tolist() + }) + + # Apply the given values for all nodes from the given physical_group + node_indices = self.mesh.get_nodes_by_physical_groups( + [physical_group_name] + )[physical_group_name] + + number_of_nodes = len(self.mesh_data.nodes) + + for node_index in node_indices: + real_index = 0 + + # Depending on the variable type and the simulation types + # the corresponding field variable may be found at different + # positions of the solution vector + if field_type is FieldType.PHI: + real_index = 2*number_of_nodes+node_index + elif field_type is FieldType.U_R: + real_index = 2*node_index + elif field_type is FieldType.U_Z: + real_index = 2*node_index+1 + else: + raise ValueError(f"Unknown variable type {field_type}") + + self.dirichlet_nodes.append(real_index) + self.dirichlet_values.append(values) + + def clear_dirichlet_bcs(self): + """Resets the dirichlet boundary conditions.""" + self.boundary_conditions = [] + self.dirichlet_nodes = [] + self.dirichlet_values = [] + + def set_electrode(self, electrode_name: str): + self.electrode_elements = self.mesh.get_elements_by_physical_groups( + [electrode_name] + )[electrode_name] + self.electrode_normals = np.tile( + [0, 1], + (len(self.electrode_elements), 1) + ) + + def assemble(self): + """Redirect to general nonlinear assembly function""" + self.material_manager.initialize_materials() + m, c, k = assemble( + self.mesh_data, + self.material_manager + ) + self.m = m + self.c = c + self.k = k + + self.nonlinearity.assemble(m, c, k) + + def simulate( + self, + delta_t: float, + number_of_time_steps: int, + gamma: float, + beta: float, + tolerance: float = 1e-11, + max_iter: int = 300, + u_start: Union[npt.NDArray, None] = None + ): + """Simulates the nonlinear piezo time system. + + Parameters: + delta_t: Difference between time steps. + number_of_time_steps: Total number of time steps. + gamma: Newmark integration parameter. + beta: Newmark integration parameter. + tolerance: Tolerance of the newton iteration. + max_iter: Maximum number of iterations for newton raphson. + u_start: Initial field for u. + """ + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) + number_of_nodes = len(self.mesh_data.nodes) + + m = self.m.copy() + c = self.c.copy() + k = self.k.copy() + + # Time integration constants + a1 = 1/(beta*(delta_t**2)) + a2 = 1/(beta*delta_t) + a3 = (1-2*beta)/(2*beta) + a4 = gamma/(beta*delta_t) + a5 = 1-gamma/beta + a6 = (1-gamma/(2*beta))*delta_t + + # Init arrays + # Displacement u which is calculated + u = np.zeros( + (3*number_of_nodes, number_of_time_steps), + dtype=np.float64 + ) + # Displacement u derived after t (du/dt) + v = np.zeros( + (3*number_of_nodes, number_of_time_steps), + dtype=np.float64 + ) + # v derived after u (d^2u/dt^2) + a = np.zeros( + (3*number_of_nodes, number_of_time_steps), + dtype=np.float64 + ) + # Charge + q = np.zeros(number_of_time_steps, dtype=np.float64) + + # Apply dirichlet bc to matrices + m, c, k = self._apply_dirichlet_bc( + m, + c, + k, + dirichlet_nodes + ) + self.nonlinearity.apply_dirichlet_bc(dirichlet_nodes) + + # Residual for the newton iterations + def residual(next_u, current_u, v, a, f): + return ( + m@(a1*(next_u-current_u)-a2*v-a3*a) + + c@(a4*(next_u-current_u)+a5*v+a6*a) + + k@next_u+self.nonlinearity.evaluate_force_vector( + next_u, + m, + c, + k + )-f + ) + + if u_start is not None: + u[:, 0] = u_start + + print("Starting nonlinear time domain simulation") + for time_index in range(number_of_time_steps-1): + # Calculate load vector + f = self._get_load_vector( + dirichlet_nodes, + dirichlet_values[:, time_index+1] + ) + + # Values of the current time step + current_u = u[:, time_index] + current_v = v[:, time_index] + current_a = a[:, time_index] + + # Current iteration value of u of the next time step + # as a start value this is set to the last converged u of the last + # time step + u_i = current_u.copy() + + # Best values during newton are saved + best_u_i = u_i.copy() + best_norm = scipy.linalg.norm(residual( + u_i, + current_u, + current_v, + current_a, + f + )) + + # Placeholder for result of newton + next_u = np.zeros(3*number_of_nodes) + self.converged = False + + if best_norm > tolerance: + # Start newton iterations + for i in range(max_iter): + # Calculate tangential stiffness matrix + tangent_matrix = self._calculate_tangent_matrix( + u_i, + m, + c, + k + ) + delta_u = linalg.spsolve( + (a1*m+a4*c+tangent_matrix), + residual(u_i, current_u, current_v, current_a, f) + ) + u_i_next = u_i - delta_u + + # Check for convergence + norm = scipy.linalg.norm(residual( + u_i_next, current_u, current_v, current_a, f + )) + + if norm < tolerance: + print( + f"Newton converged at time step {time_index} " + f"after {i+1} iteration(s)" + ) + # print(u_i_next) + next_u = u_i_next + self.converged = True + break + elif norm < best_norm: + best_norm = norm + best_u_i = u_i_next + + if i % 100 == 0 and i > 0: + print("Iteration:", i) + + # Update for next iteration + u_i = u_i_next + if not self.converged: + print( + "Newton did not converge.. Choosing best value: " + f"{best_norm}" + ) + next_u = best_u_i + else: + print("Start value norm already below tolerance") + next_u = best_u_i + + # Calculate next v and a + a[:, time_index+1] = ( + a1*(next_u-current_u) + - a2*current_v + - a3*current_a + ) + v[:, time_index+1] = ( + a4*(next_u-current_u) + + a5*current_v + + a6*current_a + ) + + # Set u array value + u[:, time_index+1] = next_u + + # Calculate charge + if (self.electrode_elements is not None + and self.electrode_elements.shape[0] > 0): + q[time_index] = calculate_charge( + u[:, time_index+1], + self.material_manager, + self.electrode_elements, + self.electrode_normals, + self.mesh_data.nodes, + self.mesh_data.element_order + ) + + self.u = u + self.q = q + + def _get_load_vector( + self, + nodes: npt.NDArray, + values: npt.NDArray + ) -> npt.NDArray: + """Calculates the load vector (right hand side) vector for the + simulation. + + Parameters: + nodes: Nodes at which the dirichlet value shall be set. + values: Dirichlet value which is set at the corresponding node. + + Returns: + Right hand side vector for the simulation. + """ + number_of_nodes = len(self.mesh_data.nodes) + + # Can be initialized to 0 because external load and volume + # charge density is 0. + f = np.zeros(3*number_of_nodes, dtype=np.float64) + + for node, value in zip(nodes, values): + f[node] = value + + return f + + def _apply_dirichlet_bc( + self, + m: sparse.lil_array, + c: sparse.lil_array, + k: sparse.lil_array, + nodes: npt.NDArray + ) -> Tuple[ + sparse.csc_array, + sparse.csc_array, + sparse.csc_array + ]: + # Set rows of matrices to 0 and diagonal of K to 1 (at node points) + m[nodes, :] = 0 + c[nodes, :] = 0 + k[nodes, :] = 0 + k[nodes, nodes] = 1 + + return m.tocsc(), c.tocsc(), k.tocsc() + + def _calculate_tangent_matrix( + self, + u: npt.NDArray, + m: sparse.csc_array, + c: sparse.csc_array, + k: sparse.csc_array + ): + # TODO Duplicate function in piezo_stationary.py + """Calculates the tangent matrix based on an analytically + expression. + + Parameters: + u: Current mechanical displacement. + k: FEM stiffness matrix. + ln: FEM nonlinear stiffness matrix. + """ + linear = k + nonlinear = self.nonlinearity.evaluate_jacobian( + u, m, c, k + ) + + return (linear + nonlinear).tocsc() diff --git a/scripts/example_nonlinear.py b/scripts/example_nonlinear.py index 9a87cfc..f3f7742 100644 --- a/scripts/example_nonlinear.py +++ b/scripts/example_nonlinear.py @@ -54,10 +54,10 @@ def create_chirp( time_values = np.arange(number_of_time_steps)*delta_t return signal.chirp( time_values, - start_frequency, - time_values[-1], - end_frequency, - method="linear" + start_frequency, + time_values[-1], + end_frequency, + method="linear" ) From a7825e5ccda805938df135b1b58b6b37a97f3334 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Wed, 29 Oct 2025 08:48:03 +0100 Subject: [PATCH 3/4] Use FEMSolver for nonlinear sim Move node_points attribute to FEMSolver. Add some docstrings for nonlinear solver. --- plutho/materials.py | 1 + plutho/simulations/__init__.py | 1 + plutho/simulations/nonlinear/__init__.py | 5 +- plutho/simulations/nonlinear/base.py | 9 ++ plutho/simulations/nonlinear/piezo_time.py | 154 +++------------------ plutho/simulations/piezo_freq.py | 8 +- plutho/simulations/piezo_time.py | 9 +- plutho/simulations/solver.py | 6 +- plutho/simulations/thermo_piezo_time.py | 8 +- plutho/simulations/thermo_time.py | 8 +- 10 files changed, 43 insertions(+), 166 deletions(-) diff --git a/plutho/materials.py b/plutho/materials.py index 021a076..aac06a3 100644 --- a/plutho/materials.py +++ b/plutho/materials.py @@ -1,6 +1,7 @@ """Contains different materials which can be used in the simulations.""" # Python standard libraries +import os import json from dataclasses import dataclass, fields from typing import List, Union diff --git a/plutho/simulations/__init__.py b/plutho/simulations/__init__.py index 1990f22..6f2f3fb 100644 --- a/plutho/simulations/__init__.py +++ b/plutho/simulations/__init__.py @@ -3,3 +3,4 @@ from .thermo_piezo_time import * from .thermo_time import * from .helpers import * +from .nonlinear import * diff --git a/plutho/simulations/nonlinear/__init__.py b/plutho/simulations/nonlinear/__init__.py index 92d6ab1..41b601f 100644 --- a/plutho/simulations/nonlinear/__init__.py +++ b/plutho/simulations/nonlinear/__init__.py @@ -1,3 +1,2 @@ -from .base import assemble, Nonlinearity, NonlinearType -from .piezo_time import NLPiezoTime -from .piezo_hb import NLPiezoHB +from .base import * +from .piezo_time import * diff --git a/plutho/simulations/nonlinear/base.py b/plutho/simulations/nonlinear/base.py index 19eab78..81e7044 100644 --- a/plutho/simulations/nonlinear/base.py +++ b/plutho/simulations/nonlinear/base.py @@ -17,6 +17,12 @@ from plutho.materials import MaterialManager from ...enums import NonlinearType + +__all__ = [ + "Nonlinearity" +] + + # # -------- Functions -------- # @@ -434,6 +440,9 @@ def _assemble_quadratic_custom(self): """Local function. Assembles the nonlinear FEM matrix for a quadratic custom nonlinearity type. """ + if self.mesh_data is None: + raise ValueError("Mesh data is not set") + nodes = self.mesh_data.nodes elements = self.mesh_data.elements element_order = self.mesh_data.element_order diff --git a/plutho/simulations/nonlinear/piezo_time.py b/plutho/simulations/nonlinear/piezo_time.py index 26c0e52..068dcbc 100644 --- a/plutho/simulations/nonlinear/piezo_time.py +++ b/plutho/simulations/nonlinear/piezo_time.py @@ -9,16 +9,26 @@ from scipy.sparse import linalg # Local libraries +from ...enums import SolverType from .base import assemble, Nonlinearity -from ...mesh.mesh import MeshData -from ...materials import MaterialData -from ...enums import FieldType -from ..helpers import calculate_charge, create_node_points -from plutho.materials import MaterialManager -from plutho.mesh.mesh import Mesh +from ...mesh.mesh import Mesh +from ..solver import FEMSolver -class NLPiezoTime: +__all__ = [ + "NLPiezoTime" +] + + +class NLPiezoTime(FEMSolver): + """Class for the simulation of nonlinear piezoeletric systems. + + Attributes: + nonlinearity: Nonlinearity object. + """ + # Nonlinear simulation + nonlinearity: Nonlinearity + def __init__( self, @@ -26,124 +36,17 @@ def __init__( mesh: Mesh, nonlinearity: Nonlinearity ): - self.simulation_name = simulation_name - - # Setup simulation data - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - # Setup mesh and material data - nodes, elements = mesh.get_mesh_nodes_and_elements() - self.mesh = mesh - self.mesh_data = MeshData(nodes, elements, mesh.element_order) - self.node_points = create_node_points( - nodes, - elements, - mesh.element_order - ) + super().__init__(simulation_name, mesh) + + self.solver_type = SolverType.PiezoTime - self.material_manager = MaterialManager(len(elements)) self.nonlinearity = nonlinearity self.nonlinearity.set_mesh_data(self.mesh_data, self.node_points) - def add_material( - self, - material_name: str, - material_data: MaterialData, - physical_group_name: str - ): - """Adds a material to the simulation - - Parameters: - material_name: Name of the material. - material_data: Material data. - physical_group_name: Name of the physical group for which the - material shall be set. If this is None or "" the material will - be set for the whole model. - """ - if physical_group_name is None or physical_group_name == "": - element_indices = np.arange(len(self.mesh_data.elements)) - else: - element_indices = self.mesh.get_elements_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - self.material_manager.add_material( - material_name, - material_data, - physical_group_name, - element_indices - ) - - def add_dirichlet_bc( - self, - field_type: FieldType, - physical_group_name: str, - values: npt.NDArray - ): - """Adds a dirichlet boundary condition (dirichlet bc) to the - simulation. The given values are for each time step. Therefore each - value is applied to every node from the given physical_group. - - Parameters: - field_type: The type of field for which the bc shall be set. - physical_group_name: Name of the physical group for which the - boundary condition shall be set. - values: List of values per time step. Value of the bc for each time - step. The value for one time step is applied to every element - equally. Length: number_of_time_step. - """ - # Save boundary condition for serialization - self.boundary_conditions.append({ - "field_type": field_type.name, - "physical_group": physical_group_name, - "values": values.tolist() - }) - - # Apply the given values for all nodes from the given physical_group - node_indices = self.mesh.get_nodes_by_physical_groups( - [physical_group_name] - )[physical_group_name] - - number_of_nodes = len(self.mesh_data.nodes) - - for node_index in node_indices: - real_index = 0 - - # Depending on the variable type and the simulation types - # the corresponding field variable may be found at different - # positions of the solution vector - if field_type is FieldType.PHI: - real_index = 2*number_of_nodes+node_index - elif field_type is FieldType.U_R: - real_index = 2*node_index - elif field_type is FieldType.U_Z: - real_index = 2*node_index+1 - else: - raise ValueError(f"Unknown variable type {field_type}") - - self.dirichlet_nodes.append(real_index) - self.dirichlet_values.append(values) - - def clear_dirichlet_bcs(self): - """Resets the dirichlet boundary conditions.""" - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - - def set_electrode(self, electrode_name: str): - self.electrode_elements = self.mesh.get_elements_by_physical_groups( - [electrode_name] - )[electrode_name] - self.electrode_normals = np.tile( - [0, 1], - (len(self.electrode_elements), 1) - ) - def assemble(self): """Redirect to general nonlinear assembly function""" self.material_manager.initialize_materials() + m, c, k = assemble( self.mesh_data, self.material_manager @@ -207,8 +110,6 @@ def simulate( (3*number_of_nodes, number_of_time_steps), dtype=np.float64 ) - # Charge - q = np.zeros(number_of_time_steps, dtype=np.float64) # Apply dirichlet bc to matrices m, c, k = self._apply_dirichlet_bc( @@ -331,20 +232,7 @@ def residual(next_u, current_u, v, a, f): # Set u array value u[:, time_index+1] = next_u - # Calculate charge - if (self.electrode_elements is not None - and self.electrode_elements.shape[0] > 0): - q[time_index] = calculate_charge( - u[:, time_index+1], - self.material_manager, - self.electrode_elements, - self.electrode_normals, - self.mesh_data.nodes, - self.mesh_data.element_order - ) - self.u = u - self.q = q def _get_load_vector( self, diff --git a/plutho/simulations/piezo_freq.py b/plutho/simulations/piezo_freq.py index d9df1ff..bdf3de9 100644 --- a/plutho/simulations/piezo_freq.py +++ b/plutho/simulations/piezo_freq.py @@ -9,7 +9,7 @@ # Local libraries from .solver import FEMSolver -from .helpers import create_node_points, calculate_volumes, \ +from .helpers import calculate_volumes, \ mat_apply_dbcs from .integrals import integral_m, integral_ku, integral_kuv, integral_kve, \ integral_loss_scs @@ -26,15 +26,11 @@ class PiezoFreq(FEMSolver): """Simulation class for frequency-domain piezoelectric simulations. Attributes: - node_points: List of node points per elements. m: Sparse mass matrix. c: Sparse damping matrix. k: Sparse stiffness matrix. mech_loss: Mechanical loss field. """ - # Internal simulation data - node_points: npt.NDArray - # FEM matrices m: sparse.lil_array c: sparse.lil_array @@ -63,9 +59,7 @@ def assemble(self): raise ValueError("Before assembly some materials must be added.") nodes = self.mesh_data.nodes - elements = self.mesh_data.elements element_order = self.mesh_data.element_order - self.node_points = create_node_points(nodes, elements, element_order) number_of_nodes = len(nodes) mu = sparse.lil_matrix( diff --git a/plutho/simulations/piezo_time.py b/plutho/simulations/piezo_time.py index 735c196..19e6e4c 100644 --- a/plutho/simulations/piezo_time.py +++ b/plutho/simulations/piezo_time.py @@ -3,12 +3,11 @@ # Third party libraries import numpy as np -import numpy.typing as npt import scipy.sparse as sparse import scipy.sparse.linalg as slin # Local libraries -from .helpers import create_node_points, mat_apply_dbcs +from .helpers import mat_apply_dbcs from .integrals import integral_m, integral_ku, integral_kuv, \ integral_kve from ..enums import SolverType @@ -25,14 +24,10 @@ class PiezoTime(FEMSolver): """Class for the simulation of time domain piezoelectric systems. Attributes: - node_points: List of node points per elements. m: Sparse mass matrix. c: Sparse damping matrix. k: Sparse stiffness matrix. """ - # Internal simulation data - node_points: npt.NDArray - # FEM matrices m: sparse.lil_array c: sparse.lil_array @@ -52,9 +47,7 @@ def assemble(self): # Maybe the 2x2 matrix slicing is not very fast self.material_manager.initialize_materials() nodes = self.mesh_data.nodes - elements = self.mesh_data.elements element_order = self.mesh_data.element_order - self.node_points = create_node_points(nodes, elements, element_order) number_of_nodes = len(nodes) mu = sparse.lil_matrix( diff --git a/plutho/simulations/solver.py b/plutho/simulations/solver.py index d180f67..b87eda2 100644 --- a/plutho/simulations/solver.py +++ b/plutho/simulations/solver.py @@ -15,7 +15,7 @@ from ..enums import SolverType, FieldType from ..materials import MaterialManager, MaterialData from ..mesh.mesh import Mesh, MeshData -from .helpers import calculate_charge +from .helpers import calculate_charge, create_node_points __all__ = [ @@ -30,6 +30,7 @@ class FEMSolver(ABC): Attributes: simulation_name: Name of the simulation. mesh: Mesh object, stores all the mesh information. + node_points: Lists of node points per element. material_manager: MaterialManager object, stores all the material informations. dirichlet_nodes: List of nodes for which a dirichlet bc is applied. @@ -45,6 +46,7 @@ class FEMSolver(ABC): """ simulation_name: str mesh: Mesh + node_points: npt.NDArray material_manager: MaterialManager dirichlet_nodes: List dirichlet_values: List @@ -57,8 +59,10 @@ def __init__(self, simulation_name: str, mesh: Mesh): self.solver_type = None nodes, elements = mesh.get_mesh_nodes_and_elements() + element_order = mesh.element_order self.mesh_data = MeshData(nodes, elements, mesh.element_order) self.material_manager = MaterialManager(len(elements)) + self.node_points = create_node_points(nodes, elements, element_order) self.dirichlet_nodes = [] self.dirichlet_values = [] diff --git a/plutho/simulations/thermo_piezo_time.py b/plutho/simulations/thermo_piezo_time.py index 5884c21..cbe95d3 100644 --- a/plutho/simulations/thermo_piezo_time.py +++ b/plutho/simulations/thermo_piezo_time.py @@ -14,7 +14,7 @@ from .solver import FEMSolver from ..enums import SolverType from ..mesh import Mesh -from .helpers import create_node_points, mat_apply_dbcs, calculate_volumes, \ +from .helpers import mat_apply_dbcs, calculate_volumes, \ get_avg_temp_field_per_element from .integrals import integral_m, integral_ku, integral_kuv, \ integral_kve, integral_ktheta, integral_theta_load, integral_loss_scs_time @@ -29,15 +29,11 @@ class ThermoPiezoTime(FEMSolver): """Class for solving time domain themo-piezoelectric systems. Attributes: - node_points: List of node points per elements. m: Sparse mass matrix. c: Sparse damping matrix. k: Sparse stiffness matrix. mech_loss: Mechanical loss field. """ - # Internal simulation data - node_points: npt.NDArray - # FEM matrices m: sparse.lil_array c: sparse.lil_array @@ -60,9 +56,7 @@ def assemble(self): # Maybe the 2x2 matrix slicing is not very fast self.material_manager.initialize_materials() nodes = self.mesh_data.nodes - elements = self.mesh_data.elements element_order = self.mesh_data.element_order - self.node_points = create_node_points(nodes, elements, element_order) number_of_nodes = len(nodes) mu = sparse.lil_matrix( diff --git a/plutho/simulations/thermo_time.py b/plutho/simulations/thermo_time.py index 882aa82..a26a499 100644 --- a/plutho/simulations/thermo_time.py +++ b/plutho/simulations/thermo_time.py @@ -15,7 +15,7 @@ from .integrals import integral_heat_flux, integral_ktheta, \ integral_theta_load, integral_m from .solver import FEMSolver -from .helpers import create_node_points, mat_apply_dbcs, \ +from .helpers import mat_apply_dbcs, \ get_avg_temp_field_per_element from ..enums import SolverType from ..mesh import Mesh @@ -30,7 +30,6 @@ class ThermoTime(FEMSolver): """Class for the simulation of time domain heat conduction problems. Attributes: - node_points: List of node points per elements. c: Sparse damping matrix. k: Sparse stiffness matrix. theta: Resulting thermal field. @@ -41,9 +40,6 @@ class ThermoTime(FEMSolver): convective_outer_temp: Temperature of the material outside of the simulated material. """ - # Internal simulation data - node_points: npt.NDArray - # FEM Matrices c: sparse.lil_array k: sparse.lil_array @@ -73,9 +69,7 @@ def assemble(self): # Maybe the 2x2 matrix slicing is not very fast self.material_manager.initialize_materials() nodes = self.mesh_data.nodes - elements = self.mesh_data.elements element_order = self.mesh_data.element_order - self.node_points = create_node_points(nodes, elements, element_order) number_of_nodes = len(nodes) c = sparse.lil_array( From 797212d46a07db7d5555eceb9bf48403f80c1d1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Wed, 29 Oct 2025 09:39:28 +0100 Subject: [PATCH 4/4] Fix bug in thermo piezo time test --- scripts/example_nonlinear.py | 2 +- tests/test_fields.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/example_nonlinear.py b/scripts/example_nonlinear.py index f3f7742..fa30763 100644 --- a/scripts/example_nonlinear.py +++ b/scripts/example_nonlinear.py @@ -122,7 +122,6 @@ def simulate_nl_time(CWD, mesh, delta_t, number_of_time_steps): "Ground", np.zeros(number_of_time_steps) ) - sim.set_electrode("Electrode") # Run simulation sim.assemble() @@ -134,6 +133,7 @@ def simulate_nl_time(CWD, mesh, delta_t, number_of_time_steps): tolerance=5e-9, max_iter=40 ) + sim.calculate_charge("Electrode") # Save results if not os.path.isdir(CWD): diff --git a/tests/test_fields.py b/tests/test_fields.py index 09ea36c..1085ca8 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -396,7 +396,7 @@ def test_thermo_piezo_time(tmp_path, test=True): ) power = np.zeros(NUMBER_OF_TIME_STEPS) for time_step in range(NUMBER_OF_TIME_STEPS): - power[time_step] = np.dot(sim.mech_loss[:, time_step], volumes) + power[time_step] = np.dot(sim.mech_loss[time_step, :], volumes) total_loss_energy = np.trapezoid( power, @@ -405,7 +405,7 @@ def test_thermo_piezo_time(tmp_path, test=True): ) # Calculate stored thermal energy - theta = sim.u[3*number_of_nodes:] + theta = sim.u[:, 3*number_of_nodes:] stored_thermal_energy = plutho.calculate_stored_thermal_energy( theta[-1, :], nodes,