diff --git a/README.md b/README.md index 5c80ac4..0f82c5e 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,7 @@ the code. Some of those feature can be discussed and are not mandatory: - Check if multiple materials already works or if continuity and boundary conditions are needed - Add test for thermal simulation with convective boundary condition +- Put all nonlinear options into one data class ### Issues diff --git a/plutho/__init__.py b/plutho/__init__.py index fce41d2..861fa9d 100644 --- a/plutho/__init__.py +++ b/plutho/__init__.py @@ -1,6 +1,7 @@ from .simulation import PiezoSimTime, ThermoPiezoSimTime, MeshData, \ SimulationData, MaterialData, SimulationType, ModelType, \ - ThermoSimTime, NonlinearType + 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, \ @@ -9,4 +10,3 @@ from .single_sim import SingleSimulation, FieldType from .mesh import Mesh from .coupled_sim import CoupledThermoPiezoThermoSim, CoupledFreqPiezoTherm -from .nonlinear_sim import PiezoNonlinear diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py deleted file mode 100644 index a688492..0000000 --- a/plutho/nonlinear_sim.py +++ /dev/null @@ -1,528 +0,0 @@ -"""Helper class to organize the usage of the piezo nonlinear simulations""" -# Python standard libraries -import os -from typing import Union, Dict, List -import configparser -import json - -# Third party libraries -import numpy as np -import numpy.typing as npt - -# Local libraries -from plutho import Mesh, MaterialManager, FieldType, MaterialData, \ - MeshData, SimulationData -from .simulation import NonlinearPiezoSimStationary, NonlinearPiezoSimTime, \ - NonlinearType - - -class PiezoNonlinear: - """Wrapper class to handle stationary and time dependent nonlinear - simulations. - - Attributes: - sim_name: Name of the simulation and the simulation folder. - sim_directory: Path to the directory in which the result files - are saved. - mesh: Mesh object containing the gmsh mesh information. - material_manger MaterialManager object. - - Parameters: - sim_name: Name of the simulation and the simulation folder which is - created. - base_directory: Path to a directory in which the simulation directory - is created. - mesh_file: Path to a mesh file which shall be used in this simulation. - """ - - # Simulation settings - sim_name: str - sim_directory: str - mesh: Mesh - material_manager: MaterialManager - solver: Union[None, NonlinearPiezoSimTime, NonlinearPiezoSimStationary] - nonlinear_material_matrix: npt.NDArray - boundary_conditions: List[Dict] - - def __init__( - self, - working_directory: str, - sim_name: str, - mesh: Mesh - ): - self.sim_name = sim_name - self.sim_directory = os.path.join(working_directory, sim_name) - - if not os.path.isdir(self.sim_directory): - os.makedirs(self.sim_directory) - - self.boundary_conditions = [] - self.dirichlet_nodes = [] - self.dirichlet_values = [] - self.nonlinear_material_matrix = np.zeros((4, 4)) - self.solver = None - self.mesh = mesh - nodes, elements = mesh.get_mesh_nodes_and_elements() - self.mesh_data = MeshData(nodes, elements) - self.material_manager = MaterialManager(len(elements)) - self.charge_calculated = False - - 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 setup_stationary_simulation(self): - """Set the simulation type to a nonlinaer stationary simulation. - """ - self.solver = NonlinearPiezoSimStationary( - self.mesh_data, - self.material_manager - ) - - def setup_time_dependent_simulation( - self, - delta_t: float, - number_of_time_steps: int, - gamma: float, - beta: float - ): - """Sets the simulation type to a nonlinear time dependent simulation. - - Parameters: - delta_t: Difference between each time step. - number_of_time_steps: Total number of time steps which are - simulated. - gamma: Time integration parameter (best value = 0.5). - beta: Time integration parameter (best value = 0.25). - """ - self.solver = NonlinearPiezoSimTime( - self.mesh_data, - self.material_manager, - SimulationData( - delta_t, - number_of_time_steps, - gamma, - beta - ) - ) - - def set_nonlinearity_type(self, ntype: NonlinearType, **kwargs): - """Sets the type of nonlinearity for the simulation. - - Parameters: - ntype: NonlinearType object. - kwargs: Can contain different parameters based on the nonlinear - type: - - If ntype=NonlinearType.Rayleigh: - "zeta" (float): Scalar nonlinearity parameter. - - If ntype=NonlinearType.Custom: - "nonlinear_matrix" (npt.NDArray): 6x6 custom nonlinear - material matrix. - """ - # TODO (almost) Duplicate code with nonlinear/piezo_time.py - if ntype is NonlinearType.Rayleigh: - if "zeta" not in kwargs: - raise ValueError( - "Missing 'zeta' parameter for nonlinear type: Rayleigh" - ) - elif ntype is NonlinearType.Custom: - if "nonlinear_matrix" not in kwargs: - raise ValueError( - "Nonlinearity matrix is missing as parameter for curstom " - "nonlinearity" - ) - else: - raise NotImplementedError( - f"Nonlinearity type {NonlinearType.value} is not implemented" - ) - - self.nonlinear_type = ntype - self.nonlinear_params = kwargs - - - def simulate( - self, - **kwargs - ): - """Runs the simulation using the previously set materials, boundary - conditions and simulation type. The simulation results are stored in - the solver object. - - Paramters: - The kwargs parameter can have the following attributes: - - "tolerance" (float): Can be set for either simulation. - Sets the tolerance - for the Newton-Raphson scheme (norm of the reisual must be - smaller in order to stop the algorithm). Typical values can be - in the order of 1e-11. - - "max_iter" (int): Can be set for either simulation. - Sets the maximum number of iterations which are done until - the algorithm stops. If the tolerance is met it is stopped - earlier. - - "electrode_elements" (npt.NDArray): Can only be set for the time - dependent simulation. If the charge shall be calculated it is - necessary to set this parameter. Set it to the indices of the - elements on which the charge shall be calculated. - - "electrode_normals" (npt.NDArray): List of normal vectors - corresponding to each electrode element (same index). - - "u_start" (npt.NDArray): Can only be set for the stationary - simulation. Gives a initial guess for u to achieve a faster - convergence. - - "alpha" (float): Can be set for either simulation. Factor which - is multiplied with the delta_u of each iteration. - Typically a value between 0 and 1. Implements a sort of - damping in each iteration which can result in a better - convergence. - - "load_factor" (float): Can only be set for the stationary - simulation. Is multiplied with the load vector. Therefore the - given load can be adjusted to achieve a faster convergence. - """ - # Check if materials are set - if self.material_manager is None: - raise ValueError( - "Cannot run nonlinear simulation. Please set a" - " material before running the simulation." - ) - - # Check if solvers are set - if self.solver is None: - raise ValueError( - "Cannot run nonlinear simulation. Please set a" - " simulation type before running the simulation" - ) - - self.material_manager.initialize_materials() - self.solver.dirichlet_nodes = np.array(self.dirichlet_nodes) - self.solver.dirichlet_values = np.array(self.dirichlet_values) - - # Assemble - self.solver.assemble( - self.nonlinear_type, - nonlinear_matrix=self.nonlinear_material_matrix - ) - - # Run simulation - if isinstance(self.solver, NonlinearPiezoSimTime): - if kwargs["electrode_elements"] is not None: - self.charge_calculated = True - self.solver.solve_time_implicit( - **kwargs - ) - elif isinstance(self.solver, NonlinearPiezoSimStationary): - self.solver.solve_newton( - **kwargs - ) - else: - raise ValueError( - "Cannot run simulation for simulation type" - f" {type(self.solver)}" - ) - - 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.sim_directory): - os.makedirs(self.sim_directory) - - if file_prefix != "": - file_prefix += "_" - - np.save( - os.path.join( - self.sim_directory, - f"{file_prefix}u.npy" - ), - self.solver.u - ) - - if isinstance(self.solver, NonlinearPiezoSimTime): - if self.charge_calculated: - np.save( - os.path.join( - self.sim_directory, - f"{file_prefix}q.npy" - ), - self.solver.q - ) - - def save_simulation_settings(self, prefix: str = ""): - if prefix != "": - prefix += "_" - - if not os.path.exists(self.sim_directory): - os.makedirs(self.sim_directory) - - settings = configparser.ConfigParser() - - if isinstance(self.solver, NonlinearPiezoSimStationary): - simulation_type = "NonlinearStationary" - elif isinstance(self.solver, NonlinearPiezoSimTime): - simulation_type = "NonlinearTime" - else: - raise ValueError( - f"Cannot save simulation type {type(self.solver)}" - ) - - # General simulation settings - general_settings = { - "name": self.sim_name, - "mesh_file": self.mesh.mesh_file_path, - "simulation_type": simulation_type - } - settings["general"] = general_settings - - # Material setings and data - material_settings = {} - for material in self.material_manager.materials: - material_settings[material.material_name] = material.to_dict() - material_settings["nonlinear_piezo_matrix"] = \ - self.nonlinear_material_matrix.tolist() - - # Simulation data - if isinstance(self.solver, NonlinearPiezoSimTime): - 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.sim_directory, - f"{prefix}{self.sim_name}.cfg" - ), - "w", - encoding="UTF-8" - ) as fd: - settings.write(fd) - - # Save materials to json - with open( - os.path.join( - self.sim_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.sim_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 = "" - ): - 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 = PiezoNonlinear( - simulation_folder, - "", - Mesh(mesh_file) - ) - # Workaround since simulation_folder and simulation_name are - # combined in the constructor, see empty string above for the - # constructor - simulation.sim_name = simulation_name - - if simulation_type == "NonlinearStationary": - simulation.setup_stationary_simulation() - elif simulation_type == "NonlinearTime": - simulation_data = SimulationData( - float(settings["simulation"]["delta_t"]), - int(settings["simulation"]["number_of_time_steps"]), - float(settings["simulation"]["gamma"]), - float(settings["simulation"]["beta"]), - ) - simulation.setup_time_dependent_simulation( - simulation_data.delta_t, - simulation_data.number_of_time_steps, - 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(): - if material_name == "nonlinear_piezo_matrix": - simulation.nonlinear_material_matrix = np.array( - material_settings[material_name] - ) - else: - 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): - if isinstance(self.solver, NonlinearPiezoSimStationary): - u_file = os.path.join( - self.sim_directory, - "u.npy" - ) - self.solver.u = np.load(u_file) - elif isinstance(self.solver, NonlinearPiezoSimTime): - u_file = os.path.join( - self.sim_directory, - "u.npy" - ) - q_file = os.path.join( - self.sim_directory, - "q.npy" - ) - - self.solver.u = np.load(u_file) - if os.path.isfile(q_file): - self.solver.q = np.load(q_file) - else: - raise ValueError( - "Cannot load simulation settings of simulation type", - type(self.solver) - ) diff --git a/plutho/simulation/__init__.py b/plutho/simulation/__init__.py index 096eb32..d404de5 100644 --- a/plutho/simulation/__init__.py +++ b/plutho/simulation/__init__.py @@ -4,5 +4,4 @@ SimulationType from .thermo_time import ThermoSimTime from .piezo_freq import PiezoSimFreq -from .nonlinear import NonlinearPiezoSimTime, NonlinearPiezoSimStationary, \ - NonlinearType +from .nonlinear import NLPiezoTime, NonlinearType, NLPiezoHB, Nonlinearity diff --git a/plutho/simulation/base.py b/plutho/simulation/base.py index 2a350e8..7b4d796 100644 --- a/plutho/simulation/base.py +++ b/plutho/simulation/base.py @@ -15,6 +15,14 @@ # -------- 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 @@ -221,7 +229,7 @@ def gradient_local_shape_functions_2d(s, t, element_order) -> npt.NDArray: ]) case 2: return np.array([ - [ # d_s + [ # d_s -3+4*t+4*s, 4*s-1, 0, @@ -229,7 +237,7 @@ def gradient_local_shape_functions_2d(s, t, element_order) -> npt.NDArray: 4*t, -4*t ], - [ # d_t + [ # d_t -3+4*s+4*t, 0, 4*t-1, @@ -417,7 +425,10 @@ def inner(s, t): ) r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return np.dot(np.dot(b_op.T, elasticity_matrix), b_op)*r*jacobian_det + return np.dot( + np.dot(b_op.T, elasticity_matrix), + b_op + ) * r * jacobian_det return quadratic_quadrature(inner, element_order) @@ -454,7 +465,10 @@ def inner(s, t): global_dn = np.dot(jacobian_inverted_t, dn) r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return np.dot(np.dot(b_op.T, piezo_matrix.T), global_dn)*r*jacobian_det + return np.dot( + np.dot(b_op.T, piezo_matrix.T), + global_dn + ) * r * jacobian_det return quadratic_quadrature(inner, element_order) @@ -491,7 +505,7 @@ def inner(s, t): permittivity_matrix ), global_dn - )*r*jacobian_det + ) * r * jacobian_det return quadratic_quadrature(inner, element_order) @@ -518,7 +532,7 @@ def inner(s, t): n = local_shape_functions_2d(s, t, element_order) r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return np.dot(n.T, theta)*r*jacobian_det + return np.dot(n.T, theta) * r * jacobian_det return quadratic_quadrature(inner, element_order) @@ -543,7 +557,7 @@ def inner(s, t): r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return r*jacobian_det + return r * jacobian_det return quadratic_quadrature(inner, element_order) @@ -815,15 +829,15 @@ def create_node_points( elements: npt.NDArray, element_order: int ) -> npt.NDArray: - """Create the local node data and the corresponding matrices - for every element which are needed in many parts of the simulations. + """Create the local node data for every given element. Parameters: - nodes: Nodes of the mesh - elements: Elements of the mesh + nodes: Nodes of the mesh. + elements: Elements of the mesh. + element_order: Order of the elements. Returns: - List of LocalElementData objects. + List of nodes. """ points_per_element = int(1/2*(element_order+1)*(element_order+2)) diff --git a/plutho/simulation/nonlinear/__init__.py b/plutho/simulation/nonlinear/__init__.py index 962fb19..92d6ab1 100644 --- a/plutho/simulation/nonlinear/__init__.py +++ b/plutho/simulation/nonlinear/__init__.py @@ -1,4 +1,3 @@ -from .base import assemble, NonlinearType -from .piezo_time import NonlinearPiezoSimTime, NonlinearType -from .piezo_stationary import NonlinearPiezoSimStationary - +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 index 1516793..3351c40 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -6,7 +6,6 @@ # Third party libraries import numpy as np -import numpy.typing as npt from scipy import sparse # Local libraries @@ -15,67 +14,29 @@ create_node_points, quadratic_quadrature, gradient_local_shape_functions_2d from plutho.materials import MaterialManager - +# +# -------- Enums -------- +# class NonlinearType(Enum): - Rayleigh = "Rayleigh" - Custom = "Custom" - + QuadraticRayleigh = 0, + QuadraticCustom = 1, + CubicRayleigh = 2, + CubicCustom = 3 +# +# -------- Functions -------- +# def assemble( - mesh_data: MeshData, - material_manager: MaterialManager, - nonlinear_type: NonlinearType, - **kwargs + mesh_data: MeshData, + material_manager: MaterialManager ) -> Tuple[ - sparse.lil_array, sparse.lil_array, sparse.lil_array, sparse.lil_array ]: - """Assembles the FEM matrices based on the set mesh_data and - material_data. Resulting FEM matrices are saved in global variables. - - Parameters: - ntype: NonlinearType object. - kwargs: Can contain different parameters based on the nonlinear - type: - - If ntype=NonlinearType.Rayleigh: - "zeta" (float): Scalar nonlinearity parameter. - - If ntype=NonlinearType.Custom: - "nonlinear_matrix" (npt.NDArray): 6x6 custom nonlinear - material matrix. - - Returns: - Tuple containing the assembled matrices m, c, k, ln - """ - - # Check for nonlinear types - if nonlinear_type is NonlinearType.Rayleigh: - if "zeta" not in kwargs: - raise ValueError( - "Missing 'zeta' parameter for nonlinear type: Rayleigh" - ) - zeta = kwargs["zeta"] - elif nonlinear_type is NonlinearType.Custom: - if "nonlinear_matrix" not in kwargs: - raise ValueError( - "Missing 'nonlinear_matrix' parameter for nonlinaer type:" - " Rayleigh" - ) - nm = kwargs["nonlinear_matrix"] - nonlinear_matrix = np.array([ - [nm[0, 0], nm[0, 2], 0, nm[0, 1]], - [nm[0, 2], nm[2, 2], 0, nm[0, 2]], - [0, 0, nm[3, 3], 0], - [nm[0, 1], nm[0, 2], 0, nm[0, 0]] - ]) - else: - raise NotImplementedError( - f"Nonlinearity type {NonlinearType.value} is not implemented" - ) - # 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 @@ -98,10 +59,6 @@ def assemble( (number_of_nodes, number_of_nodes), dtype=np.float64 ) - lu = sparse.lil_matrix( - (2*number_of_nodes, 2*number_of_nodes), - dtype=np.float64 - ) for element_index, element in enumerate(elements): current_node_points = node_points[element_index] @@ -139,14 +96,6 @@ def assemble( element_order ) * 2 * np.pi ) - if nonlinear_type is NonlinearType.Custom: - lu_e = ( - integral_u_nonlinear( - current_node_points, - nonlinear_matrix, - element_order - ) * 2 * np.pi - ) # Now assemble all element matrices for local_p, global_p in enumerate(element): @@ -171,16 +120,6 @@ def assemble( # be set directly. kv[global_p, global_q] += kve_e[local_p, local_q] - if nonlinear_type is NonlinearType.Custom: - # L_e is similar to Ku_e - lu[ - 2*global_p:2*global_p+2, - 2*global_q:2*global_q+2 - ] += lu_e[ - 2*local_p:2*local_p+2, - 2*local_q:2*local_q+2 - ] - # Calculate damping matrix # Currently in the simulations only one material is used # TODO Add algorithm to calculate cu for every element on its own @@ -190,10 +129,6 @@ def assemble( + material_manager.get_alpha_k(0) * ku ) - if nonlinear_type is NonlinearType.Rayleigh: - # Set Rayleigh type nonlinear matrix - lu = zeta*ku - # Calculate block matrices zeros1x1 = np.zeros((number_of_nodes, number_of_nodes)) zeros2x1 = np.zeros((2*number_of_nodes, number_of_nodes)) @@ -211,32 +146,15 @@ def assemble( [ku, kuv], [kuv.T, -1*kv] ]).tolil() - ln = sparse.block_array([ - [lu, zeros2x1], - [zeros1x2, zeros1x1] - ]).tolil() - return m, c, k, ln + return m, c, k -def integral_u_nonlinear( - node_points: npt.NDArray, - nonlinear_elasticity_matrix: npt.NDArray, - element_order: int +def integral_u_nonlinear_quadratic( + node_points, + nonlinear_elasticity_tensor, + element_order ): - """Calculates the Ku integral - - Parameters: - node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of - one triangle. - jacobian_inverted_t: Jacobian matrix inverted and transposed, needed - for calculation of global derivatives. - elasticity_matrix: Elasticity matrix for the current element - (c matrix). - - Returns: - npt.NDArray: 6x6 Ku 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) @@ -252,12 +170,241 @@ def inner(s, t): ) r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return np.dot( - np.dot( - b_op.T, - nonlinear_elasticity_matrix - ), - b_op - ) * 1/2 * r * jacobian_det + 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 new file mode 100644 index 0000000..6c098dd --- /dev/null +++ b/plutho/simulation/nonlinear/piezo_hb.py @@ -0,0 +1,505 @@ +"""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_stationary.py b/plutho/simulation/nonlinear/piezo_stationary.py deleted file mode 100644 index 1321636..0000000 --- a/plutho/simulation/nonlinear/piezo_stationary.py +++ /dev/null @@ -1,264 +0,0 @@ -"""Module for the simulation of nonlinear piezoeletric systems -for the case of a time-stationary simulation. (d_t=0) -""" - -# Python standard libraries -from typing import Union, Tuple - -# 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 NonlinearPiezoSimStationary: - """Implements a nonlinear time stationary simulation. The nonlinearity - is embed in the hooke law: T_i=C_ij*S_j+1/2*L_ijk*S_j*S_k. - In order to solve this system a Newton-Raphson algorithm is implemented. - For comparison a simulation can also be done using the scipy least squares - function. Additionaly a the corresponding linear system (L=0) can be solved - too. - """ - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - - # Boundary conditions - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray - - # FEM matrices - k: sparse.lil_array - ln: sparse.lil_array - - # Resulting fields - u: npt.NDArray - - def __init__( - self, - mesh_data: MeshData, - material_manager: MaterialManager, - ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.dirichlet_nodes = np.array([]) - self.dirichlet_valuse = np.array([]) - - def assemble(self, nonlinear_type: NonlinearType, **kwargs): - """Redirect to general nonlinear assembly function""" - _, _, k, ln = assemble( - self.mesh_data, - self.material_manager, - nonlinear_type, - **kwargs - ) - self.k = k - self.ln = ln - - def solve_linear(self): - """Solves the linear problem Ku=f (ln is not used). - """ - # Get FEM matrices - k = self.k.copy() - ln = self.ln.copy() # Not used - - # Apply boundary conditions - k, ln = NonlinearPiezoSimStationary.apply_dirichlet_bc( - k, - ln, - 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 - ) - - def solve_newton( - self, - tolerance: float = 1e-11, - max_iter: int = 1000, - 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: - 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. - """ - # Get FEM matrices - k = self.k.copy() - ln = self.ln.copy() - - # Apply boundary conditions - k, ln = NonlinearPiezoSimStationary.apply_dirichlet_bc( - k, - ln, - 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( - k.tocsc(), - f - ) - else: - current_u = u_start - - # Residual of the Newton algorithm - def residual(u): - return k@u+ln@np.square(u)-load_factor*f - - # Add a best found field parameter which can be returned when the - # maximum number of iterations were done - best = current_u.copy() - best_norm = np.linalg.norm(residual(best)) - - # Check if the initial value already suffices the tolerance condition - if best_norm < tolerance: - return best - - for iteration_count in range(max_iter): - # Calculate tangential stiffness matrix - k_tangent = NonlinearPiezoSimStationary. \ - calculate_tangent_matrix_hadamard( - current_u, - k, - ln - ) - - # Solve for displacement increment - delta_u = linalg.spsolve( - k_tangent, - residual(current_u) - ) - - # Update step - next_u = current_u - alpha * delta_u - - # Check for convergence - norm = scipy.linalg.norm(residual(next_u)) - if norm < tolerance: - print( - f"Solve Newton found solution after {iteration_count} " - f"steps with residual {norm}" - ) - return next_u - elif norm < best_norm: - best_norm = norm - best = 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 - - 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 - -# -------- Static functions -------- - - @staticmethod - def calculate_tangent_matrix_hadamard( - u: npt.NDArray, - k: sparse.sparray, - ln: sparse.sparray - ): - # TODO Duplicate function in piezo_time.py - """Calculates the tangent matrix based on an analytically - expression. - - Parameters: - u: Current mechanical displacement. - k: FEM stiffness matrix. - ln: FEM nonlinear stiffness matrix. - """ - k_tangent = k+2*ln@sparse.diags_array(u) - - return k_tangent - - @staticmethod - def apply_dirichlet_bc( - k: sparse.lil_array, - ln: sparse.lil_array, - nodes: npt.NDArray - ) -> Tuple[sparse.csc_array, sparse.csc_array]: - """Sets dirichlet boundary condition for the given matrix. - - Parameters: - k: Stiffness matrix. - l: Nonlinear stiffness matrix. - nodes: List of node indices at which the bc are defined. - - Returns: - Matrix with set boundary conditions. - """ - for node in nodes: - k[node, :] = 0 - ln[node, :] = 0 - - k[node, node] = 1 - - return k.tocsc(), ln.tocsc() diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py index e614b32..a4703fc 100644 --- a/plutho/simulation/nonlinear/piezo_time.py +++ b/plutho/simulation/nonlinear/piezo_time.py @@ -1,7 +1,7 @@ """Module for the simulation of nonlinear piezoelectric systems""" # Third party libraries -from typing import Tuple +from typing import Tuple, Union import numpy as np import numpy.typing as npt import scipy @@ -9,108 +9,174 @@ from scipy.sparse import linalg # Local libraries -from .base import assemble, NonlinearType -from ..base import SimulationData, MeshData +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 NonlinearPiezoSimTime: - """Simulates eletro-mechanical fields for the given mesh, material and - simulation information. The simulation is using a nonlinear stiffness - matrix which is used for a quadratic term in the hooke law. - - Attributes: - mesh_data: Contains the mesh information. - material_manager: Organizes material parameters. Can use materials - with different temperatures. - simulation_data: Contains settings for the simulation. - dirichlet_nodes_u: Nodes at which a dirichlet bc shall be set for the - mechanical field. - dirichlet_values_u: Values which are set at the node from the nodes - list for the mechanical field. - dirichlet_nodes_phi: Nodes at which a dirichlet bc shall be set for - the electrical field. - dirichlet_values_phi: Values which are set at the node from the nodes - list for the electrical field. - m: FEM mass matrix for the mechanical field. - k: FEM stiffness matrix for the mechanical field. - c: FEM damping matrix for the mechanical field. - ln: FEM stiffness matrix for the nonlinear part of the - mechanical field. - u: Mechanical field vector u(element, t). - q: Charge q(t). - - Parameters: - mesh_data: MeshData object containing the mesh. - material_manager: MaterialManager object containing the set materials - for each region. - simulation_data:SimulationData object containing information on the - time domain simulation. - """ - - # Simulation parameters - mesh_data: MeshData - material_manager: MaterialManager - simulation_data: SimulationData - - # Dirichlet boundary condition - dirichlet_nodes: npt.NDArray - dirichlet_values: npt.NDArray - - # Resulting fields - u: npt.NDArray - q: npt.NDArray - - # FEM matrices - m: sparse.lil_array - c: sparse.lil_array - ln: sparse.lil_array - k: sparse.lil_array +class NLPiezoTime: def __init__( self, - mesh_data: MeshData, - material_manager: MaterialManager, - simulation_data: SimulationData + simulation_name: str, + mesh: Mesh, + nonlinearity: Nonlinearity ): - self.mesh_data = mesh_data - self.material_manager = material_manager - self.simulation_data = simulation_data + 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 - self.dirichlet_nodes = np.array([]) - self.dirichlet_values = np.array([]) + 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 assemble(self, nonlinear_type: NonlinearType, **kwargs): + 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""" - m, c, k, ln = assemble( + self.material_manager.initialize_materials() + m, c, k = assemble( self.mesh_data, - self.material_manager, - nonlinear_type, - **kwargs + self.material_manager ) self.m = m self.c = c self.k = k - self.ln = ln - def solve_time_implicit( + 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, - electrode_elements: npt.NDArray = np.array([]), - electrode_normals: npt.NDArray = np.array([]) + u_start: Union[npt.NDArray, None] = None ): - number_of_time_steps = self.simulation_data.number_of_time_steps - delta_t = self.simulation_data.delta_t + 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) - beta = self.simulation_data.beta - gamma = self.simulation_data.gamma m = self.m.copy() c = self.c.copy() k = self.k.copy() - ln = self.ln.copy() # Time integration constants a1 = 1/(beta*(delta_t**2)) @@ -140,30 +206,38 @@ def solve_time_implicit( q = np.zeros(number_of_time_steps, dtype=np.float64) # Apply dirichlet bc to matrices - m, c, k, ln = NonlinearPiezoSimTime.apply_dirichlet_bc( + m, c, k = self._apply_dirichlet_bc( m, c, k, - ln, - self.dirichlet_nodes + 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( - self.dirichlet_nodes, - self.dirichlet_values[:, time_index+1] + f = self._get_load_vector( + dirichlet_nodes, + dirichlet_values[:, time_index+1] ) - # 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+ln@np.square(next_u)-f - ) - # Values of the current time step current_u = u[:, time_index] current_v = v[:, time_index] @@ -186,59 +260,56 @@ def residual(next_u, current_u, v, a, f): # Placeholder for result of newton next_u = np.zeros(3*number_of_nodes) - converged = False + self.converged = False if best_norm > tolerance: # Start newton iterations for i in range(max_iter): # Calculate tangential stiffness matrix - k_tangent = NonlinearPiezoSimTime. \ - calculate_tangent_matrix_hadamard( - u_i, - k, - ln - ) + tangent_matrix = self._calculate_tangent_matrix( + u_i, + m, + c, + k + ) delta_u = linalg.spsolve( - (a1*m+a4*c+k_tangent), + (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) - ) - # print( - # f"Time step: {time_index}, iteration: {i}, " - # f"norm: {norm}" - # ) + 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.copy() - converged = True + next_u = u_i_next + self.converged = True break elif norm < best_norm: best_norm = norm - best_u_i = u_i_next.copy() + 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.copy() - if not converged: + 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.copy() + next_u = best_u_i else: print("Start value norm already below tolerance") - next_u = best_u_i.copy() + next_u = best_u_i # Calculate next v and a a[:, time_index+1] = ( @@ -256,13 +327,13 @@ def residual(next_u, current_u, v, a, f): u[:, time_index+1] = next_u # Calculate charge - if (electrode_elements is not None - and electrode_elements.shape[0] > 0): + 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, - electrode_elements, - electrode_normals, + self.electrode_elements, + self.electrode_normals, self.mesh_data.nodes, self.mesh_data.element_order ) @@ -270,8 +341,7 @@ def residual(next_u, current_u, v, a, f): self.u = u self.q = q - - def get_load_vector( + def _get_load_vector( self, nodes: npt.NDArray, values: npt.NDArray @@ -297,40 +367,31 @@ def get_load_vector( return f - @staticmethod - def apply_dirichlet_bc( + def _apply_dirichlet_bc( + self, m: sparse.lil_array, c: sparse.lil_array, k: sparse.lil_array, - ln: sparse.lil_array, nodes: npt.NDArray ) -> Tuple[ - sparse.csc_array, sparse.csc_array, sparse.csc_array, sparse.csc_array ]: - # TODO Parameters are not really ndarrays # 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 - # Matrices for u_r component - for node in nodes: - # Set rows to 0 - m[node, :] = 0 - c[node, :] = 0 - k[node, :] = 0 - ln[node, :] = 0 + return m.tocsc(), c.tocsc(), k.tocsc() - # Set diagonal values to 1 - k[node, node] = 1 - - return m.tocsc(), c.tocsc(), k.tocsc(), ln.tocsc() - - @staticmethod - def calculate_tangent_matrix_hadamard( + def _calculate_tangent_matrix( + self, u: npt.NDArray, - k: sparse.csc_array, - ln: sparse.csc_array + 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 @@ -341,6 +402,9 @@ def calculate_tangent_matrix_hadamard( k: FEM stiffness matrix. ln: FEM nonlinear stiffness matrix. """ - k_tangent = k+2*ln@sparse.diags_array(u) + linear = k + nonlinear = self.nonlinearity.evaluate_jacobian( + u, m, c, k + ) - return k_tangent + return (linear + nonlinear).tocsc() diff --git a/plutho/simulation/piezo_freq.py b/plutho/simulation/piezo_freq.py index 24f17f3..deeacfc 100644 --- a/plutho/simulation/piezo_freq.py +++ b/plutho/simulation/piezo_freq.py @@ -305,7 +305,7 @@ def solve_frequency( + 1j*angular_frequency*c + k ) - + f = self.get_load_vector( self.dirichlet_nodes, self.dirichlet_values, diff --git a/plutho/single_sim.py b/plutho/single_sim.py index 3bbbdf6..efc0222 100644 --- a/plutho/single_sim.py +++ b/plutho/single_sim.py @@ -1,4 +1,5 @@ """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 @@ -16,20 +17,12 @@ PiezoSimFreq from .materials import MaterialData, MaterialManager from .mesh import Mesh - +from .simulation.base import FieldType class SimulationException(Exception): """Custom exception to simplify errors.""" -class FieldType(Enum): - """Possible field types which are calculated using differnet simulations. - """ - U_R = 0 - U_Z = 1 - PHI = 2 - THETA = 3 - class SingleSimulation: """Base class to handle single simulations. Multiple different simulation diff --git a/pyproject.toml b/pyproject.toml index 0a4ec32..cf951bd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ classifiers = [ "Programming Language ::; Python :: 3", "Topic :: Scientific/Engineering" ] -requires-python = ">=3.10,<3.13" +requires-python = ">=3.10" dependencies = [ "numpy", "matplotlib", diff --git a/scripts/example_nonlinear.py b/scripts/example_nonlinear.py index 4ddeef3..9a87cfc 100644 --- a/scripts/example_nonlinear.py +++ b/scripts/example_nonlinear.py @@ -5,6 +5,8 @@ # Third party libraries import numpy as np +import matplotlib.pyplot as plt +from scipy import signal # Local libraries import plutho @@ -43,6 +45,22 @@ ]) +def create_chirp( + start_frequency, + end_frequency, + number_of_time_steps, + delta_t +): + time_values = np.arange(number_of_time_steps)*delta_t + return signal.chirp( + time_values, + start_frequency, + time_values[-1], + end_frequency, + method="linear" + ) + + def create_sinusoidal_excitation( amplitude, delta_t, @@ -62,26 +80,22 @@ def create_sinusoidal_excitation( return amplitude*np.sin(2*np.pi*frequency*time_steps) -def simulate_nonlinear_stationary(CWD): - sim_name = "nonlinear_stationary_sim" - mesh_file = os.path.join(CWD, "ring_mesh.msh") +def simulate_nl_time(CWD, mesh, delta_t, number_of_time_steps): + # Simulation parameters + GAMMA = 0.5 + BETA = 0.25 + ZETA = 10 + AMPLITUDE = 10 + FREQUENCY = 2.07e6 - # Load/create ring mesh - if not os.path.exists(mesh_file): - plutho.Mesh.generate_rectangular_mesh( - mesh_file, - width=0.00635, - height=0.001, - x_offset=0.0026, - mesh_size=0.0001 - ) - mesh = plutho.Mesh(mesh_file) + nonlinearity = plutho.Nonlinearity() + nonlinearity.set_cubic_rayleigh(ZETA) # Create simulation - sim = plutho.PiezoNonlinear( - sim_name, - CWD, - mesh + sim = plutho.NLPiezoTime( + nl_time_sim_name, + mesh, + nonlinearity ) # Set materials @@ -90,110 +104,64 @@ def simulate_nonlinear_stationary(CWD): material_data=pic181, physical_group_name="" # Means all elements ) - sim.set_nonlinearity_type( - plutho.NonlinearType.Custom, - nonlinear_matrix=c_nonlin_6x6_nonsymmetric - ) # Set boundary conditions - excitation = 1000 + excitation = create_sinusoidal_excitation( + amplitude=AMPLITUDE, + delta_t=delta_t, + number_of_time_steps=number_of_time_steps, + frequency=FREQUENCY + ) sim.add_dirichlet_bc( plutho.FieldType.PHI, "Electrode", - np.array(excitation) + excitation ) sim.add_dirichlet_bc( plutho.FieldType.PHI, "Ground", - np.array(0) + np.zeros(number_of_time_steps) ) + sim.set_electrode("Electrode") - sim.setup_stationary_simulation() - - # For better convergence set u_r and u_z for one node to 0 - sim.dirichlet_nodes.append(0) - sim.dirichlet_nodes.append(1) - sim.dirichlet_values.append(0) - sim.dirichlet_values.append(0) - + # Run simulation + sim.assemble() sim.simulate( - tolerance=1e-10 + delta_t, + number_of_time_steps, + GAMMA, + BETA, + tolerance=5e-9, + max_iter=40 ) + # Save results + if not os.path.isdir(CWD): + os.makedirs(CWD) -def simulate_nonlinaer_time_dep(CWD): - sim_name = "nonlinear_time_dep_sim" - mesh_file = os.path.join(CWD, "ring_mesh.msh") - - # Load/create ring mesh - if not os.path.exists(mesh_file): - plutho.Mesh.generate_rectangular_mesh( - mesh_file, - width=0.00635, - height=0.001, - x_offset=0.0026, - mesh_size=0.0001 - ) - mesh = plutho.Mesh(mesh_file) + u_file = os.path.join(CWD, "u.npy") + q_file = os.path.join(CWD, "q.npy") + np.save(u_file, sim.u) + np.save(q_file, sim.q) - # Create simulation - sim = plutho.PiezoNonlinear( - sim_name, - CWD, - mesh - ) - # Set materials - sim.add_material( - material_name="pic181", - material_data=pic181, - physical_group_name="" # Means all elements - ) - sim.set_nonlinearity_type( - plutho.NonlinearType.Custom, - nonlinear_matrix=c_nonlin_6x6_nonsymmetric - ) - - # Simulation parameters - DELTA_T = 1e-8 - NUMBER_OF_TIME_STEPS = 1000 - - # Set boundary conditions - excitation = create_sinusoidal_excitation( - amplitude=10, - delta_t=DELTA_T, - number_of_time_steps=NUMBER_OF_TIME_STEPS, - frequency=2.07e6 - ) - sim.add_dirichlet_bc( - plutho.FieldType.PHI, - "Electrode", - excitation - ) - sim.add_dirichlet_bc( - plutho.FieldType.PHI, - "Ground", - np.zeros(NUMBER_OF_TIME_STEPS) - ) +def plot_displacement_spectrum( + working_directory, + delta_t, + number_of_time_steps +): + node_index = 129 - sim.setup_time_dependent_simulation( - delta_t=DELTA_T, - number_of_time_steps=NUMBER_OF_TIME_STEPS, - gamma=0.5, - beta=0.25 - ) + u = np.load(os.path.join(working_directory, "u.npy")) + u_r = u[2*node_index, :] + u_z = u[2*node_index+1, :] + U_r_jw = np.fft.fft(u_r) - # Needed in order to get the charge - electrode_elements = mesh.get_elements_by_physical_groups( - ["Electrode"] - )["Electrode"] - electrode_normals = np.tile([0, 1], (len(electrode_elements), 1)) + frequencies = np.fft.fftfreq(number_of_time_steps, delta_t) - # Run simulation - sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals - ) + plt.plot(frequencies, np.abs(U_r_jw)) + plt.grid() + plt.show() if __name__ == "__main__": @@ -205,5 +173,30 @@ def simulate_nonlinaer_time_dep(CWD): if not os.path.isdir(CWD): os.makedirs(CWD) - simulate_nonlinear_stationary(CWD) - simulate_nonlinaer_time_dep(CWD) + + # Load/create ring mesh + mesh_file = os.path.join(CWD, "ring_mesh.msh") + if not os.path.exists(mesh_file): + plutho.Mesh.generate_rectangular_mesh( + mesh_file, + width=0.00635, + height=0.001, + x_offset=0.0026, + mesh_size=0.0001 + ) + mesh = plutho.Mesh(mesh_file, element_order=1) + + DELTA_T = 1e-9 + NUMBER_OF_TIME_STEPS = 20000 + + nl_time_sim_name = "nonlinear_time_dep_sim_20k_1e-9" + nl_time_wd = os.path.join(CWD, nl_time_sim_name) + + ## Simulate + if True: + # simulate_nonlinear_stationary(CWD) + simulate_nl_time(nl_time_wd, mesh, DELTA_T, NUMBER_OF_TIME_STEPS) + + ## Plot + if False: + plot_displacement_spectrum(nl_time_wd, DELTA_T, NUMBER_OF_TIME_STEPS) diff --git a/tests/test_fields.py b/tests/test_fields.py index f5b00f1..bf17e6d 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -10,8 +10,8 @@ # -------- Global variables -------- NUMBER_OF_TIME_STEPS = 1000 -ATOL = 1e-15 -RTOL = 1e-7 +ATOL = 1e-14 +RTOL = 1e-4 pic255 = plutho.MaterialData( **{