From 7fb65c5a6df66eef749bef2d4319683018994d26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 17 Jul 2025 16:45:25 +0200 Subject: [PATCH 01/17] Update solver to use full nonlinearity matrices Add a function for an nonlinear stiffness matrix which uses an analytic quadratic differentiation and a function for a calculating it numerically. Add an example to use quadratic nonlinearities. --- plutho/nonlinear_sim.py | 10 +- plutho/simulation/base.py | 14 +- plutho/simulation/nonlinear/base.py | 304 ++++++++++++++++++++++++++-- scripts/nonlinear_quadratic.py | 208 +++++++++++++++++++ 4 files changed, 508 insertions(+), 28 deletions(-) create mode 100644 scripts/nonlinear_quadratic.py diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index a688492..26dd987 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -41,7 +41,7 @@ class PiezoNonlinear: mesh: Mesh material_manager: MaterialManager solver: Union[None, NonlinearPiezoSimTime, NonlinearPiezoSimStationary] - nonlinear_material_matrix: npt.NDArray + nonlinear_material_matrix: Union[None, npt.NDArray] boundary_conditions: List[Dict] def __init__( @@ -59,11 +59,13 @@ def __init__( self.boundary_conditions = [] self.dirichlet_nodes = [] self.dirichlet_values = [] - self.nonlinear_material_matrix = np.zeros((4, 4)) self.solver = None + self.nonlinear_material_matrix = None + self.mesh = mesh nodes, elements = mesh.get_mesh_nodes_and_elements() - self.mesh_data = MeshData(nodes, elements) + self.mesh_data = MeshData(nodes, elements, mesh.element_order) + self.material_manager = MaterialManager(len(elements)) self.charge_calculated = False @@ -279,7 +281,7 @@ def simulate( # Assemble self.solver.assemble( self.nonlinear_type, - nonlinear_matrix=self.nonlinear_material_matrix + **self.nonlinear_params ) # Run simulation diff --git a/plutho/simulation/base.py b/plutho/simulation/base.py index 2a350e8..4f4b93f 100644 --- a/plutho/simulation/base.py +++ b/plutho/simulation/base.py @@ -221,7 +221,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 +229,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, @@ -815,15 +815,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/base.py b/plutho/simulation/nonlinear/base.py index 1516793..4cf160d 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -2,7 +2,7 @@ # Python standard libraries from enum import Enum -from typing import Tuple +from typing import Tuple, Callable # Third party libraries import numpy as np @@ -10,7 +10,7 @@ from scipy import sparse # Local libraries -from ..base import MeshData, local_to_global_coordinates, b_operator_global, \ +from ..base import MeshData, local_shape_functions_2d, 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 @@ -59,16 +59,19 @@ def assemble( elif nonlinear_type is NonlinearType.Custom: if "nonlinear_matrix" not in kwargs: raise ValueError( - "Missing 'nonlinear_matrix' parameter for nonlinaer type:" + "Missing 'nonlinear_matrix' parameter for nonlinear type:" " Rayleigh" ) + # Reduce to axisymmetric matrix 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]] - ]) + 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] = nm[i_old, j_old, k_old] + + nonlinear_matrix = nm_reduced else: raise NotImplementedError( f"Nonlinearity type {NonlinearType.value} is not implemented" @@ -224,7 +227,7 @@ def integral_u_nonlinear( nonlinear_elasticity_matrix: npt.NDArray, element_order: int ): - """Calculates the Ku integral + """Calculates the nonlinear stiffness integral Parameters: node_points: List of node points [[x1, x2, x3], [y1, y2, y3]] of @@ -235,7 +238,7 @@ def integral_u_nonlinear( (c matrix). Returns: - npt.NDArray: 6x6 Ku matrix for the given element. + npt.NDArray: 6x6 matrix for the given element. """ def inner(s, t): dn = gradient_local_shape_functions_2d(s, t, element_order) @@ -243,6 +246,71 @@ def inner(s, t): jacobian_inverted_t = np.linalg.inv(jacobian).T jacobian_det = np.linalg.det(jacobian) + def t_i(s, t): + b_op = b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + + return np.einsum( + "ijk,jl,kl->il", + nonlinear_elasticity_matrix, + b_op, + b_op + ) + + t_derived = b_opt_global_t_numerical( + node_points, + t_i, + s, + t, + jacobian_inverted_t, + element_order + ) + + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + N = np.diag( + local_shape_functions_2d(s, t, element_order) + ) + + return np.outer( + t_derived, + N + ) * 1/2 * r * jacobian_det + + return quadratic_quadrature(inner, element_order) + + +def integral_u_nonlinear_analytic( + node_points: npt.NDArray, + nonlinear_elasticity_matrix: npt.NDArray, + element_order: int +): + """Calculates the nonlinear stiffness 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 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) + + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + b_op = b_operator_global( node_points, jacobian_inverted_t, @@ -250,14 +318,216 @@ def inner(s, t): t, element_order ) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - return np.dot( - np.dot( - b_op.T, - nonlinear_elasticity_matrix - ), + bb_op = double_b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + + left = np.einsum( + "ijk,jl,kl->il", + nonlinear_elasticity_matrix, + bb_op, b_op + ) + + right = np.einsum( + "ijk,jl,kl->il", + nonlinear_elasticity_matrix, + b_op, + bb_op + ) + + return np.outer( + left+right, + local_shape_functions_2d(s, t, element_order) ) * 1/2 * r * jacobian_det return quadratic_quadrature(inner, element_order) + + +def b_opt_global_t_numerical( + node_points, + func: Callable, + s, + t, + jacobian_inverted_t, + element_order +): + # Idea: + # f contains the function values at the nodes + # first calculate f at the given positions s and t + # then derive f at that point after s and t + # then calculate the derivatives after r and z + + # f should have shape 4 x nodes_per_element + + # TODO Check if this value is suitable + h = 1/20 + + nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) + + # Points at which the function is evaluated + p1 = [s+h,t] + p2 = [s-h,t] + p3 = [s,t+h] + p4 = [s,t-h] + + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + # Local derivatives + f = func(s, t) + ds_f = 1/(2*h)*(func(*p1)-func(*p2)) + dt_f = 1/(2*h)*(func(*p3)-func(*p4)) + + dr_f = np.zeros(shape=ds_f.shape) + dz_f = np.zeros(shape=dt_f.shape) + b_num = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) + + for i in range(ds_f.shape[0]): + for j in range(ds_f.shape[1]): + # TODO Not all derivatives are needed + # Derivatives with respect to local coordinates s and t + local_der = np.array([ds_f[i][j], dt_f[i][j]]) + + # Derivatives with respect to global coordinates r and z + global_der = np.dot(jacobian_inverted_t, local_der) + dr_f[i, j] = global_der[0] + dz_f[i, j] = global_der[1] + + for node_index in range(nodes_per_element): + current_f = f[:, 2*node_index:2*(node_index+1)] + dr = dr_f[:, 2*node_index:2*(node_index+1)] + dz = dz_f[:, 2*node_index:2*(node_index+1)] + + # Get global derivatives + b_num[2*node_index:2*(node_index+1), 2*node_index:2*(node_index+1)] = \ + [ + [ + dr[0][0]+dz[2][1]+current_f[3][0]/r, + dr[0][1]+dz[2][1]+current_f[3][1] + ], + [ + dz[1][0]+dr[2][0], + dz[1][1]+dr[2][1] + ] + ] + + return b_num + + +def second_gradient_local_shape_functions_2d( + s, + t, + element_order +) -> npt.NDArray: + # Since the second gradient has to be done for s and t as well, 3 indices + # are needed [i,j,k]. The first two represent if its derived after s [i=0] + # or t [i=1] and the last index [k] represents the node index of the + # element. + nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) + match element_order: + case 1: + return np.zeros(shape=(2, 2, nodes_per_element)) + case 2: + return np.array([ + [ # (outer) d_s + [4, 4, 0, -8, 0, 0], # (inner) d_s + [4, 0, 0, -4, 4, -4] # (inner) d_t + ], + [ # (outer) d_t + [4, 0, 0, -4, 4, -4], # (inner) d_s + [4, 0, 4, 0, 0, -8] # (inner) d_t + ], + ]) + case 3: + raise NotImplementedError("Not yet implemented for order 3") + + raise ValueError( + "Second gradient of shape functions not implemented for element " + f"order {element_order}" + ) + +def double_b_operator_global( + node_points: npt.NDArray, + jacobian_inverted_t: npt.NDArray, + s: float, + t: float, + element_order: int +): + nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) + + dn = gradient_local_shape_functions_2d(s, t, element_order) + global_dn = np.dot(jacobian_inverted_t, dn) + ddn = second_gradient_local_shape_functions_2d(s, t, element_order) + + # Calculate derivatives with respect to global coordinates + global_ddn = np.zeros(shape=ddn.shape) + + for node_index in range(nodes_per_element): + # Get hessian with shape + # [ds*ds, ds*dt] + # [dt*ds, dt*dt] + # Local hessian since the derivatives are with respect to s and t + local_hessian = np.array([ + [global_ddn[0][0][node_index], global_ddn[0][1][node_index]], + [global_ddn[1][0][node_index], global_ddn[1][1][node_index]] + ]) + + # Global hessian since the derivatives are with respect to r and z + # TODO Verify + global_hessian = np.dot( + np.dot( + jacobian_inverted_t, + local_hessian + ), + jacobian_inverted_t.T + ) + + # TODO Make using slices + global_ddn[0][0][node_index] = global_hessian[0][0] + global_ddn[0][1][node_index] = global_hessian[0][1] + global_ddn[1][0][node_index] = global_hessian[1][0] + global_ddn[1][1][node_index] = global_hessian[1][1] + + r = local_to_global_coordinates(node_points, s, t, element_order)[0] + + bb = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) + for i in range(nodes_per_element): + for j in range(nodes_per_element): + if i == j: + # Diagonal elements have second order derivatives + bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ + [ + global_ddn[0][0][i] + global_ddn[1][1][i]+1/(r**2), + global_ddn[1][0][i] + ], + [ + global_ddn[0][1][i], + global_ddn[0][0][i] + global_ddn[1][1][i] + ] + ] + else: + # Offdiagional elements have 2 first order derivatives + bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ + [ + ( + global_dn[0][i]*global_dn[0][j] + + global_dn[1][i]*global_dn[1][j] + + 1/(r**2) + ), + global_dn[1][i]*global_dn[1][j] + ], + [ + global_dn[0][i]*global_dn[1][j], + ( + global_dn[0][i]*global_dn[0][j] + + global_dn[1][i]*global_dn[1][j] + ) + ] + ] + + return bb diff --git a/scripts/nonlinear_quadratic.py b/scripts/nonlinear_quadratic.py new file mode 100644 index 0000000..4a5d971 --- /dev/null +++ b/scripts/nonlinear_quadratic.py @@ -0,0 +1,208 @@ +"""Implements an example for quadratic nonlinearities in time domain""" + +# Python standard libraries +import os + +# Third party libraries +import numpy as np + +# Local libraries +import plutho + + +c_nonlin_6x6x6 = np.array([ + [[-1.9362e+14, -1.6882e+13, -1.6882e+13, 0.0000e+00, 0.0000e+00, 0.0000e+00], + [-1.6882e+13, -1.6882e+13, 3.1629e+15, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [-1.6882e+13, 3.1629e+15, -1.6882e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, - + 1.5899e+15, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, -4.4184e+13, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13]], + + [[-1.6882e+13, -1.6882e+13, 3.1629e+15, 0.0000e+00, 0.0000e+00, 0.0000e+00], + [-1.6882e+13, -1.9362e+14, -1.6882e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [3.1629e+15, -1.6882e+13, -1.6882e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, - + 4.4184e+13, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, -1.5899e+15, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13]], + + [[-1.6882e+13, 3.1629e+15, -1.6882e+13, 0.0000e+00, 0.0000e+00, 0.0000e+00], + [3.1629e+15, -1.6882e+13, -1.6882e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [-1.6882e+13, -1.6882e+13, -1.9362e+14, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, - + 4.4184e+13, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, -4.4184e+13, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -1.5899e+15]], + + [[0.0000e+00, 0.0000e+00, 0.0000e+00, -1.5899e+15, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, - + 4.4184e+13, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, - + 4.4184e+13, 0.0000e+00, 0.0000e+00], + [-1.5899e+15, -4.4184e+13, -4.4184e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, 0.0000e+00, 7.7285e+14], + [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 7.7285e+14, 0.0000e+00]], + + [[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, -1.5899e+15, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, -4.4184e+13, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, 0.0000e+00, 7.7285e+14], + [-4.4184e+13, -1.5899e+15, -4.4184e+13, + 0.0000e+00, 0.0000e+00, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, 7.7285e+14, 0.0000e+00, 0.0000e+00]], + + [[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, 0.0000e+00, -4.4184e+13], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, 0.0000e+00, -1.5899e+15], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 0.0000e+00, 7.7285e+14, 0.0000e+00], + [0.0000e+00, 0.0000e+00, 0.0000e+00, + 7.7285e+14, 0.0000e+00, 0.0000e+00], + [-4.4184e+13, -4.4184e+13, -1.5899e+15, 0.0000e+00, 0.0000e+00, 0.0000e+00]] +]) + +pic181 = plutho.MaterialData( + **{ + "c11": 141218192791.12772, + "c12": 82292914124.0279, + "c13": 80362329625.75212, + "c33": 137188538228.6311, + "c44": 29848846049.816402, + "e15": 13.216444117013664, + "e31": -4.979419636068149, + "e33": 14.149818966822629, + "eps11": 1.3327347064648263e-08, + "eps33": 5.380490373139249e-09, + "alpha_m": 0.0, + "alpha_k": 1.289813815258054e-10, + "thermal_conductivity": 1.1, + "heat_capacity": 350, + "temperatures": 25, + "density": 7850 + } +) + + +def create_sinusoidal_excitation( + amplitude, + delta_t, + number_of_time_steps, + frequency +): + """Creates a sinusoidal excitation array. + + Parameters: + amplitude: Amplitude of the sin function. + delta_t: Time difference between each time step. + number_of_time_steps: Number of time steps of the excitation. Typically + this is the same as for the simulation + frequency: Frequency of the sin function. + """ + time_steps = np.arange(number_of_time_steps)*delta_t + return amplitude*np.sin(2*np.pi*frequency*time_steps) + + +def simulate(CWD): + sim_name = "nonlinear_time_quadratic" + mesh_file = os.path.join(CWD, "ring_mesh.msh") + element_order = 2 + + # 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.0002, + element_order=element_order + ) + mesh = plutho.Mesh(mesh_file, element_order) + + # 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_6x6x6 + ) + + # Simulation parameters + DELTA_T = 1e-8 + NUMBER_OF_TIME_STEPS = 1000 + + # Set boundary conditions + excitation = create_sinusoidal_excitation( + amplitude=1, + 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) + ) + + sim.setup_time_dependent_simulation( + delta_t=DELTA_T, + number_of_time_steps=NUMBER_OF_TIME_STEPS, + gamma=0.5, + beta=0.25 + ) + + # 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)) + + # Run simulation + sim.simulate( + electrode_elements=electrode_elements, + electrode_normals=electrode_normals + ) + + +if __name__ == "__main__": + CWD = os.path.join( + os.path.dirname(os.path.abspath(__file__)), + "simulations" + ) + + if not os.path.isdir(CWD): + os.makedirs(CWD) + + simulate(CWD) From 724ffe3215caf683c1f2baa2ca43fa644194d5ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 18 Jul 2025 12:54:32 +0200 Subject: [PATCH 02/17] Fix numeric B_T calculation and update quadratic nonlinear simulation --- plutho/nonlinear_sim.py | 11 +++++-- plutho/simulation/nonlinear/base.py | 46 +++++++++++++++++++---------- scripts/nonlinear_quadratic.py | 22 +++++++++----- 3 files changed, 53 insertions(+), 26 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 26dd987..6bf8d64 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -41,7 +41,6 @@ class PiezoNonlinear: mesh: Mesh material_manager: MaterialManager solver: Union[None, NonlinearPiezoSimTime, NonlinearPiezoSimStationary] - nonlinear_material_matrix: Union[None, npt.NDArray] boundary_conditions: List[Dict] def __init__( @@ -362,8 +361,14 @@ def save_simulation_settings(self, prefix: str = ""): 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() + material_settings["nonlinear_type"] = self.nonlinear_type.value + + match self.nonlinear_type: + case NonlinearType.Rayleigh: + material_settings["zeta"] = self.nonlinear_params["zeta"] + case NonlinearType.Custom: + material_settings["nonlinear_matrix"] = \ + self.nonlinear_params["nonlinear_matrix"].tolist() # Simulation data if isinstance(self.solver, NonlinearPiezoSimTime): diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py index 4cf160d..3e15f3d 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -350,13 +350,27 @@ def inner(s, t): def b_opt_global_t_numerical( - node_points, + node_points: npt.NDArray, func: Callable, - s, - t, - jacobian_inverted_t, - element_order + s: float, + t: float, + jacobian_inverted_t: npt.NDArray, + element_order: int ): + """Calculates the transposed B operator on the given func numerically (by + using a central finite difference scheme). + The result is evaluated at the given local coordinates s and t. + The derivatives are done with respect to the global coordinates r and z. + + Parameters: + node_points: Node points of the triangle. + func: Function on which the B_T operator is applied. + s: Local coordinate at which the B_T operator is evaluated. + t: Local coordinate at which the B_T operator is evaluated. + jacobian_inverted_T: Inverted and transposed jacobian of the current + triangle, needed for calculating global derivatives. + element_order: Element order of the triangle. + """ # Idea: # f contains the function values at the nodes # first calculate f at the given positions s and t @@ -366,7 +380,7 @@ def b_opt_global_t_numerical( # f should have shape 4 x nodes_per_element # TODO Check if this value is suitable - h = 1/20 + h = 1/100 nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) @@ -385,13 +399,13 @@ def b_opt_global_t_numerical( dr_f = np.zeros(shape=ds_f.shape) dz_f = np.zeros(shape=dt_f.shape) - b_num = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) + b_num = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) for i in range(ds_f.shape[0]): for j in range(ds_f.shape[1]): # TODO Not all derivatives are needed # Derivatives with respect to local coordinates s and t - local_der = np.array([ds_f[i][j], dt_f[i][j]]) + local_der = np.array([ds_f[i, j], dt_f[i, j]]) # Derivatives with respect to global coordinates r and z global_der = np.dot(jacobian_inverted_t, local_der) @@ -399,20 +413,20 @@ def b_opt_global_t_numerical( dz_f[i, j] = global_der[1] for node_index in range(nodes_per_element): - current_f = f[:, 2*node_index:2*(node_index+1)] - dr = dr_f[:, 2*node_index:2*(node_index+1)] - dz = dz_f[:, 2*node_index:2*(node_index+1)] + current_f = f[:, 2*node_index:2*node_index+2] + dr = dr_f[:, 2*node_index:2*node_index+2] + dz = dz_f[:, 2*node_index:2*node_index+2] # Get global derivatives - b_num[2*node_index:2*(node_index+1), 2*node_index:2*(node_index+1)] = \ + b_num[2*node_index:2*node_index+2, 2*node_index:2*node_index+2] = \ [ [ - dr[0][0]+dz[2][1]+current_f[3][0]/r, - dr[0][1]+dz[2][1]+current_f[3][1] + dr[0, 0]+dz[2, 1]+current_f[3, 0]/r, + dr[0, 1]+dz[2, 1]+current_f[3, 1]/r ], [ - dz[1][0]+dr[2][0], - dz[1][1]+dr[2][1] + dz[1, 0]+dr[2, 0], + dz[1, 1]+dr[2, 1] ] ] diff --git a/scripts/nonlinear_quadratic.py b/scripts/nonlinear_quadratic.py index 4a5d971..c1ac806 100644 --- a/scripts/nonlinear_quadratic.py +++ b/scripts/nonlinear_quadratic.py @@ -120,9 +120,9 @@ def create_sinusoidal_excitation( def simulate(CWD): - sim_name = "nonlinear_time_quadratic" - mesh_file = os.path.join(CWD, "ring_mesh.msh") - element_order = 2 + sim_name = "nonlinear_time_quadratic_1v" + mesh_file = os.path.join(CWD, "ring_mesh_third_order.msh") + element_order = 3 # Load/create ring mesh if not os.path.exists(mesh_file): @@ -138,8 +138,8 @@ def simulate(CWD): # Create simulation sim = plutho.PiezoNonlinear( - sim_name, CWD, + sim_name, mesh ) @@ -151,12 +151,16 @@ def simulate(CWD): ) sim.set_nonlinearity_type( plutho.NonlinearType.Custom, - nonlinear_matrix=c_nonlin_6x6x6 + nonlinear_matrix=c_nonlin_6x6x6/1e3 ) + #sim.set_nonlinearity_type( + # plutho.NonlinearType.Rayleigh, + # zeta=1e3 + #) # Simulation parameters DELTA_T = 1e-8 - NUMBER_OF_TIME_STEPS = 1000 + NUMBER_OF_TIME_STEPS = 10000 # Set boundary conditions excitation = create_sinusoidal_excitation( @@ -192,9 +196,13 @@ def simulate(CWD): # Run simulation sim.simulate( electrode_elements=electrode_elements, - electrode_normals=electrode_normals + electrode_normals=electrode_normals, + tolerance=1e-10 ) + sim.save_simulation_settings() + sim.save_simulation_results() + if __name__ == "__main__": CWD = os.path.join( From 9dbce8a6fe6966dcb79dcbf46a9ce1fb9b8d6124 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 18 Jul 2025 14:39:36 +0200 Subject: [PATCH 03/17] Update nonlinear save/load --- plutho/nonlinear_sim.py | 62 +++++++++++++++++++++++++--------- scripts/nonlinear_quadratic.py | 48 +++++++++++++++++++++----- 2 files changed, 85 insertions(+), 25 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 6bf8d64..6a9096b 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -353,6 +353,7 @@ def save_simulation_settings(self, prefix: str = ""): general_settings = { "name": self.sim_name, "mesh_file": self.mesh.mesh_file_path, + "mesh_order": self.mesh_data.mesh_order, "simulation_type": simulation_type } settings["general"] = general_settings @@ -361,15 +362,18 @@ def save_simulation_settings(self, prefix: str = ""): material_settings = {} for material in self.material_manager.materials: material_settings[material.material_name] = material.to_dict() - material_settings["nonlinear_type"] = self.nonlinear_type.value + nonlinear_materials = {} + nonlinear_materials["nonlinear_type"] = self.nonlinear_type.value match self.nonlinear_type: case NonlinearType.Rayleigh: - material_settings["zeta"] = self.nonlinear_params["zeta"] + nonlinear_materials["zeta"] = self.nonlinear_params["zeta"] case NonlinearType.Custom: - material_settings["nonlinear_matrix"] = \ + nonlinear_materials["nonlinear_matrix"] = \ self.nonlinear_params["nonlinear_matrix"].tolist() + material_settings["nonlinear"] = nonlinear_materials + # Simulation data if isinstance(self.solver, NonlinearPiezoSimTime): settings["simulation"] = self.solver.simulation_data.__dict__ @@ -440,17 +444,19 @@ def load_simulation_settings( # Read general simulation settings mesh_file = settings["general"]["mesh_file"] + mesh_order = int(settings["general"]["mesh_order"]) simulation_type = settings["general"]["simulation_type"] simulation = PiezoNonlinear( simulation_folder, "", - Mesh(mesh_file) + Mesh(mesh_file, mesh_order) ) # 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": @@ -474,21 +480,45 @@ def load_simulation_settings( # Read materials with open(necessary_files[1], "r", encoding="UTF-8") as fd: material_settings = json.load(fd) + + # Load nonlinearities + nonlinear_material = material_settings["nonlinear"] + nonlinear_type = NonlinearType( + nonlinear_material["nonlinear_type"] + ) + + match nonlinear_type: + case NonlinearType.Rayleigh: + nonlinear_params = { + "zeta": float(nonlinear_material["zeta"]) + } + case NonlinearType.Custom: + nonlinear_params = { + "nonlinear_matrix": np.array( + nonlinear_material["nonlinear_matrix"] + ) + } + + # Load other materials 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"] - ) + if material_name == "nonlinear": + continue + + simulation.add_material( + material_name, + MaterialData( + **material_settings[material_name]["material_data"] + ), + material_settings[material_name]["physical_group_name"] + ) + simulation.material_manager.initialize_materials() + simulation.set_nonlinearity_type( + nonlinear_type, + **nonlinear_params + ) + # Read boundary conditions with open(necessary_files[2], "r", encoding="UTF-8") as fd: bc_settings = json.load(fd) diff --git a/scripts/nonlinear_quadratic.py b/scripts/nonlinear_quadratic.py index c1ac806..4d75520 100644 --- a/scripts/nonlinear_quadratic.py +++ b/scripts/nonlinear_quadratic.py @@ -5,6 +5,8 @@ # Third party libraries import numpy as np +import scipy +import matplotlib.pyplot as plt # Local libraries import plutho @@ -116,11 +118,18 @@ def create_sinusoidal_excitation( frequency: Frequency of the sin function. """ time_steps = np.arange(number_of_time_steps)*delta_t - return amplitude*np.sin(2*np.pi*frequency*time_steps) + sin = amplitude*np.sin(2*np.pi*frequency*time_steps) + window = scipy.signal.windows.tukey( + number_of_time_steps + ) + excitation = np.multiply( + sin, + window + ) + return excitation -def simulate(CWD): - sim_name = "nonlinear_time_quadratic_1v" +def simulate(CWD, simulation_name): mesh_file = os.path.join(CWD, "ring_mesh_third_order.msh") element_order = 3 @@ -139,7 +148,7 @@ def simulate(CWD): # Create simulation sim = plutho.PiezoNonlinear( CWD, - sim_name, + simulation_name, mesh ) @@ -151,7 +160,7 @@ def simulate(CWD): ) sim.set_nonlinearity_type( plutho.NonlinearType.Custom, - nonlinear_matrix=c_nonlin_6x6x6/1e3 + nonlinear_matrix=c_nonlin_6x6x6 ) #sim.set_nonlinearity_type( # plutho.NonlinearType.Rayleigh, @@ -159,15 +168,15 @@ def simulate(CWD): #) # Simulation parameters - DELTA_T = 1e-8 + DELTA_T = 4e-8 NUMBER_OF_TIME_STEPS = 10000 # Set boundary conditions excitation = create_sinusoidal_excitation( - amplitude=1, + amplitude=10, delta_t=DELTA_T, number_of_time_steps=NUMBER_OF_TIME_STEPS, - frequency=2.07e6 + frequency=0.12261e6 ) sim.add_dirichlet_bc( plutho.FieldType.PHI, @@ -204,6 +213,25 @@ def simulate(CWD): sim.save_simulation_results() +def plot_displacement_spectrum(simulation_folder): + sim = plutho.PiezoNonlinear.load_simulation_settings(simulation_folder) + sim.load_simulation_results() + + node_index = 129 + + u_r = sim.solver.u[2*node_index, :] + u_z = sim.solver.u[2*node_index+1, :] + U_r_jw = np.fft.fft(u_r) + + number_of_time_steps = sim.solver.simulation_data.number_of_time_steps + delta_t = sim.solver.simulation_data.delta_t + frequencies = np.fft.fftfreq(number_of_time_steps, delta_t) + + plt.plot(frequencies, np.abs(U_r_jw)) + plt.grid() + plt.show() + + if __name__ == "__main__": CWD = os.path.join( os.path.dirname(os.path.abspath(__file__)), @@ -213,4 +241,6 @@ def simulate(CWD): if not os.path.isdir(CWD): os.makedirs(CWD) - simulate(CWD) + sim_name = "nonlinear_time_quadratic_1v" + # simulate(CWD, sim_name) + plot_displacement_spectrum(os.path.join(CWD, sim_name)) From cb7900b143384c7c31dfd4711557289881a3582e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 18 Jul 2025 14:43:48 +0200 Subject: [PATCH 04/17] Remove nonlinear quadratic script --- scripts/nonlinear_quadratic.py | 246 --------------------------------- 1 file changed, 246 deletions(-) delete mode 100644 scripts/nonlinear_quadratic.py diff --git a/scripts/nonlinear_quadratic.py b/scripts/nonlinear_quadratic.py deleted file mode 100644 index 4d75520..0000000 --- a/scripts/nonlinear_quadratic.py +++ /dev/null @@ -1,246 +0,0 @@ -"""Implements an example for quadratic nonlinearities in time domain""" - -# Python standard libraries -import os - -# Third party libraries -import numpy as np -import scipy -import matplotlib.pyplot as plt - -# Local libraries -import plutho - - -c_nonlin_6x6x6 = np.array([ - [[-1.9362e+14, -1.6882e+13, -1.6882e+13, 0.0000e+00, 0.0000e+00, 0.0000e+00], - [-1.6882e+13, -1.6882e+13, 3.1629e+15, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [-1.6882e+13, 3.1629e+15, -1.6882e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - - 1.5899e+15, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, -4.4184e+13, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13]], - - [[-1.6882e+13, -1.6882e+13, 3.1629e+15, 0.0000e+00, 0.0000e+00, 0.0000e+00], - [-1.6882e+13, -1.9362e+14, -1.6882e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [3.1629e+15, -1.6882e+13, -1.6882e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - - 4.4184e+13, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, -1.5899e+15, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13]], - - [[-1.6882e+13, 3.1629e+15, -1.6882e+13, 0.0000e+00, 0.0000e+00, 0.0000e+00], - [3.1629e+15, -1.6882e+13, -1.6882e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [-1.6882e+13, -1.6882e+13, -1.9362e+14, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - - 4.4184e+13, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, -4.4184e+13, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -1.5899e+15]], - - [[0.0000e+00, 0.0000e+00, 0.0000e+00, -1.5899e+15, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - - 4.4184e+13, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - - 4.4184e+13, 0.0000e+00, 0.0000e+00], - [-1.5899e+15, -4.4184e+13, -4.4184e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, 0.0000e+00, 7.7285e+14], - [0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 7.7285e+14, 0.0000e+00]], - - [[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, -1.5899e+15, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, -4.4184e+13, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, 0.0000e+00, 7.7285e+14], - [-4.4184e+13, -1.5899e+15, -4.4184e+13, - 0.0000e+00, 0.0000e+00, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, 7.7285e+14, 0.0000e+00, 0.0000e+00]], - - [[0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, -4.4184e+13], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, 0.0000e+00, -4.4184e+13], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, 0.0000e+00, -1.5899e+15], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 0.0000e+00, 7.7285e+14, 0.0000e+00], - [0.0000e+00, 0.0000e+00, 0.0000e+00, - 7.7285e+14, 0.0000e+00, 0.0000e+00], - [-4.4184e+13, -4.4184e+13, -1.5899e+15, 0.0000e+00, 0.0000e+00, 0.0000e+00]] -]) - -pic181 = plutho.MaterialData( - **{ - "c11": 141218192791.12772, - "c12": 82292914124.0279, - "c13": 80362329625.75212, - "c33": 137188538228.6311, - "c44": 29848846049.816402, - "e15": 13.216444117013664, - "e31": -4.979419636068149, - "e33": 14.149818966822629, - "eps11": 1.3327347064648263e-08, - "eps33": 5.380490373139249e-09, - "alpha_m": 0.0, - "alpha_k": 1.289813815258054e-10, - "thermal_conductivity": 1.1, - "heat_capacity": 350, - "temperatures": 25, - "density": 7850 - } -) - - -def create_sinusoidal_excitation( - amplitude, - delta_t, - number_of_time_steps, - frequency -): - """Creates a sinusoidal excitation array. - - Parameters: - amplitude: Amplitude of the sin function. - delta_t: Time difference between each time step. - number_of_time_steps: Number of time steps of the excitation. Typically - this is the same as for the simulation - frequency: Frequency of the sin function. - """ - time_steps = np.arange(number_of_time_steps)*delta_t - sin = amplitude*np.sin(2*np.pi*frequency*time_steps) - window = scipy.signal.windows.tukey( - number_of_time_steps - ) - excitation = np.multiply( - sin, - window - ) - return excitation - - -def simulate(CWD, simulation_name): - mesh_file = os.path.join(CWD, "ring_mesh_third_order.msh") - element_order = 3 - - # 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.0002, - element_order=element_order - ) - mesh = plutho.Mesh(mesh_file, element_order) - - # Create simulation - sim = plutho.PiezoNonlinear( - CWD, - simulation_name, - 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_6x6x6 - ) - #sim.set_nonlinearity_type( - # plutho.NonlinearType.Rayleigh, - # zeta=1e3 - #) - - # Simulation parameters - DELTA_T = 4e-8 - NUMBER_OF_TIME_STEPS = 10000 - - # Set boundary conditions - excitation = create_sinusoidal_excitation( - amplitude=10, - delta_t=DELTA_T, - number_of_time_steps=NUMBER_OF_TIME_STEPS, - frequency=0.12261e6 - ) - sim.add_dirichlet_bc( - plutho.FieldType.PHI, - "Electrode", - excitation - ) - sim.add_dirichlet_bc( - plutho.FieldType.PHI, - "Ground", - np.zeros(NUMBER_OF_TIME_STEPS) - ) - - sim.setup_time_dependent_simulation( - delta_t=DELTA_T, - number_of_time_steps=NUMBER_OF_TIME_STEPS, - gamma=0.5, - beta=0.25 - ) - - # 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)) - - # Run simulation - sim.simulate( - electrode_elements=electrode_elements, - electrode_normals=electrode_normals, - tolerance=1e-10 - ) - - sim.save_simulation_settings() - sim.save_simulation_results() - - -def plot_displacement_spectrum(simulation_folder): - sim = plutho.PiezoNonlinear.load_simulation_settings(simulation_folder) - sim.load_simulation_results() - - node_index = 129 - - u_r = sim.solver.u[2*node_index, :] - u_z = sim.solver.u[2*node_index+1, :] - U_r_jw = np.fft.fft(u_r) - - number_of_time_steps = sim.solver.simulation_data.number_of_time_steps - delta_t = sim.solver.simulation_data.delta_t - frequencies = np.fft.fftfreq(number_of_time_steps, delta_t) - - plt.plot(frequencies, np.abs(U_r_jw)) - plt.grid() - plt.show() - - -if __name__ == "__main__": - CWD = os.path.join( - os.path.dirname(os.path.abspath(__file__)), - "simulations" - ) - - if not os.path.isdir(CWD): - os.makedirs(CWD) - - sim_name = "nonlinear_time_quadratic_1v" - # simulate(CWD, sim_name) - plot_displacement_spectrum(os.path.join(CWD, sim_name)) From ee39233136a93dd546c25f512b527acbbff40fc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 18 Jul 2025 14:52:19 +0200 Subject: [PATCH 05/17] Fix mesh order in .cfg file --- plutho/nonlinear_sim.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 6a9096b..2dfc013 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -353,7 +353,7 @@ def save_simulation_settings(self, prefix: str = ""): general_settings = { "name": self.sim_name, "mesh_file": self.mesh.mesh_file_path, - "mesh_order": self.mesh_data.mesh_order, + "element_order": self.mesh_data.element_order, "simulation_type": simulation_type } settings["general"] = general_settings @@ -444,12 +444,12 @@ def load_simulation_settings( # Read general simulation settings mesh_file = settings["general"]["mesh_file"] - mesh_order = int(settings["general"]["mesh_order"]) + element_order = int(settings["general"]["element_order"]) simulation_type = settings["general"]["simulation_type"] simulation = PiezoNonlinear( simulation_folder, "", - Mesh(mesh_file, mesh_order) + Mesh(mesh_file, element_order) ) # Workaround since simulation_folder and simulation_name are # combined in the constructor, see empty string above for the From 7cb6fafebcfb1838939f8f5372e8227f02902353 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Wed, 23 Jul 2025 16:39:40 +0200 Subject: [PATCH 06/17] Add nonlinear_order parameter --- plutho/nonlinear_sim.py | 17 ++++++++++--- plutho/simulation/nonlinear/base.py | 25 +++++++++++++++---- .../simulation/nonlinear/piezo_stationary.py | 9 ++++++- plutho/simulation/nonlinear/piezo_time.py | 9 ++++++- 4 files changed, 49 insertions(+), 11 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 2dfc013..1694057 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -188,7 +188,15 @@ def setup_time_dependent_simulation( ) ) - def set_nonlinearity_type(self, ntype: NonlinearType, **kwargs): + def set_nonlinearity_order(self, nl_order: int): + if nl_order < 1 or nl_order > 2: + raise NotImplementedError( + "Nonlinearity order only implemented for 2 or 3" + ) + + self.nonlinear_order = nl_order + + def set_nonlinearity_type(self, nl_type: NonlinearType, **kwargs): """Sets the type of nonlinearity for the simulation. Parameters: @@ -202,12 +210,12 @@ def set_nonlinearity_type(self, ntype: NonlinearType, **kwargs): material matrix. """ # TODO (almost) Duplicate code with nonlinear/piezo_time.py - if ntype is NonlinearType.Rayleigh: + if nl_type is NonlinearType.Rayleigh: if "zeta" not in kwargs: raise ValueError( "Missing 'zeta' parameter for nonlinear type: Rayleigh" ) - elif ntype is NonlinearType.Custom: + elif nl_type is NonlinearType.Custom: if "nonlinear_matrix" not in kwargs: raise ValueError( "Nonlinearity matrix is missing as parameter for curstom " @@ -218,7 +226,7 @@ def set_nonlinearity_type(self, ntype: NonlinearType, **kwargs): f"Nonlinearity type {NonlinearType.value} is not implemented" ) - self.nonlinear_type = ntype + self.nonlinear_type = nl_type self.nonlinear_params = kwargs @@ -279,6 +287,7 @@ def simulate( # Assemble self.solver.assemble( + self.nonlinear_order, self.nonlinear_type, **self.nonlinear_params ) diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py index 3e15f3d..ece0243 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -22,10 +22,11 @@ class NonlinearType(Enum): def assemble( - mesh_data: MeshData, - material_manager: MaterialManager, - nonlinear_type: NonlinearType, - **kwargs + mesh_data: MeshData, + material_manager: MaterialManager, + nonlinear_order: int, + nonlinear_type: NonlinearType, + **kwargs ) -> Tuple[ sparse.lil_array, sparse.lil_array, @@ -62,8 +63,22 @@ def assemble( "Missing 'nonlinear_matrix' parameter for nonlinear type:" " Rayleigh" ) - # Reduce to axisymmetric matrix nm = kwargs["nonlinear_matrix"] + + # Check for matrix shape + if len(nm.shape) != nonlinear_order+2: + raise ValueError( + f"Given nonlinearity matrix shape {nm.shape} does not fit to" + f" the given nonlinearity order {nonlinear_order}" + ) + + if nonlinear_order == 3: + raise NotImplementedError( + "Nonlinearity order 3 not implemented yet" + ) + + # Reduce to axisymmetric 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): diff --git a/plutho/simulation/nonlinear/piezo_stationary.py b/plutho/simulation/nonlinear/piezo_stationary.py index 1321636..e68c4a0 100644 --- a/plutho/simulation/nonlinear/piezo_stationary.py +++ b/plutho/simulation/nonlinear/piezo_stationary.py @@ -51,14 +51,21 @@ def __init__( self.dirichlet_nodes = np.array([]) self.dirichlet_valuse = np.array([]) - def assemble(self, nonlinear_type: NonlinearType, **kwargs): + def assemble( + self, + nonlinear_order: int, + nonlinear_type: NonlinearType, + **kwargs + ): """Redirect to general nonlinear assembly function""" _, _, k, ln = assemble( self.mesh_data, self.material_manager, + nonlinear_order, nonlinear_type, **kwargs ) + self.nonlinear_order = nonlinear_order self.k = k self.ln = ln diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py index e614b32..476cbba 100644 --- a/plutho/simulation/nonlinear/piezo_time.py +++ b/plutho/simulation/nonlinear/piezo_time.py @@ -81,14 +81,21 @@ def __init__( self.dirichlet_nodes = np.array([]) self.dirichlet_values = np.array([]) - def assemble(self, nonlinear_type: NonlinearType, **kwargs): + def assemble( + self, + nonlinear_order: int, + nonlinear_type: NonlinearType, + **kwargs + ): """Redirect to general nonlinear assembly function""" m, c, k, ln = 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 From 569b55a9b15eaf4c24d287a418187a8a6bf5dbcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 24 Jul 2025 08:04:29 +0200 Subject: [PATCH 07/17] Add variable nonlinear initial field --- plutho/nonlinear_sim.py | 16 ++++++++++++---- plutho/simulation/nonlinear/piezo_time.py | 14 +++++++++----- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 2dfc013..b9357e2 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -278,10 +278,18 @@ def simulate( self.solver.dirichlet_values = np.array(self.dirichlet_values) # Assemble - self.solver.assemble( - self.nonlinear_type, - **self.nonlinear_params - ) + if "skip_assemble" in kwargs: + if not kwargs["skip_assemble"]: + self.solver.assemble( + self.nonlinear_type, + **self.nonlinear_params + ) + del kwargs["skip_assemble"] + else: + self.solver.assemble( + self.nonlinear_type, + **self.nonlinear_params + ) # Run simulation if isinstance(self.solver, NonlinearPiezoSimTime): diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py index e614b32..37b023f 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 @@ -99,7 +99,8 @@ def solve_time_implicit( tolerance: float = 1e-11, max_iter: int = 300, electrode_elements: npt.NDArray = np.array([]), - electrode_normals: npt.NDArray = np.array([]) + electrode_normals: npt.NDArray = np.array([]), + u0: Union[npt.NDArray, None] = None ): number_of_time_steps = self.simulation_data.number_of_time_steps delta_t = self.simulation_data.delta_t @@ -148,6 +149,9 @@ def solve_time_implicit( self.dirichlet_nodes ) + if u0 is not None: + u[:, 0] = u0 + print("Starting nonlinear time domain simulation") for time_index in range(number_of_time_steps-1): # Calculate load vector @@ -186,7 +190,7 @@ 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 @@ -219,7 +223,7 @@ def residual(next_u, current_u, v, a, f): ) # print(u_i_next) next_u = u_i_next.copy() - converged = True + self.converged = True break elif norm < best_norm: best_norm = norm @@ -230,7 +234,7 @@ def residual(next_u, current_u, v, a, f): # Update for next iteration u_i = u_i_next.copy() - if not converged: + if not self.converged: print( "Newton did not converge.. Choosing best value: " f"{best_norm}" From 30f8b8cc4e5606505e34d23906732776092effad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 25 Jul 2025 08:35:24 +0200 Subject: [PATCH 08/17] WIP Harmonic Balancing Solver --- plutho/__init__.py | 2 +- plutho/simulation/__init__.py | 2 +- plutho/simulation/nonlinear/__init__.py | 2 +- plutho/simulation/nonlinear/piezo_hb.py | 324 ++++++++++++++++++ .../simulation/nonlinear/piezo_stationary.py | 4 +- 5 files changed, 329 insertions(+), 5 deletions(-) create mode 100644 plutho/simulation/nonlinear/piezo_hb.py diff --git a/plutho/__init__.py b/plutho/__init__.py index fce41d2..5d4a585 100644 --- a/plutho/__init__.py +++ b/plutho/__init__.py @@ -1,6 +1,6 @@ from .simulation import PiezoSimTime, ThermoPiezoSimTime, MeshData, \ SimulationData, MaterialData, SimulationType, ModelType, \ - ThermoSimTime, NonlinearType + ThermoSimTime, NonlinearType, NonlinearPiezoSimHb 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, \ diff --git a/plutho/simulation/__init__.py b/plutho/simulation/__init__.py index 096eb32..d031cc9 100644 --- a/plutho/simulation/__init__.py +++ b/plutho/simulation/__init__.py @@ -5,4 +5,4 @@ from .thermo_time import ThermoSimTime from .piezo_freq import PiezoSimFreq from .nonlinear import NonlinearPiezoSimTime, NonlinearPiezoSimStationary, \ - NonlinearType + NonlinearType, NonlinearPiezoSimHb diff --git a/plutho/simulation/nonlinear/__init__.py b/plutho/simulation/nonlinear/__init__.py index 962fb19..8ad9c91 100644 --- a/plutho/simulation/nonlinear/__init__.py +++ b/plutho/simulation/nonlinear/__init__.py @@ -1,4 +1,4 @@ from .base import assemble, NonlinearType from .piezo_time import NonlinearPiezoSimTime, NonlinearType from .piezo_stationary import NonlinearPiezoSimStationary - +from .piezo_hb import NonlinearPiezoSimHb diff --git a/plutho/simulation/nonlinear/piezo_hb.py b/plutho/simulation/nonlinear/piezo_hb.py new file mode 100644 index 0000000..2c306e3 --- /dev/null +++ b/plutho/simulation/nonlinear/piezo_hb.py @@ -0,0 +1,324 @@ +"""Module for the simulation of nonlinaer piezoelectric systems using the +harmonic balancing method.""" + +# 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 NonlinearPiezoSimHb: + """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 + l: 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_order: int, + nonlinear_type: NonlinearType, + harmonic_order: int, + **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, l = assemble( + self.mesh_data, + self.material_manager, + nonlinear_order, + nonlinear_type, + **kwargs + ) + self.nonlinear_order = nonlinear_order + self.harmonic_order = harmonic_order + + self.m = m + self.c = c + self.k = k + self.l = l + + def solve_linear(self): + """Solves the linear problem Ku=f (ln is not used). + """ + # Get FEM matrices + k = self.k.copy() + ln = self.l.copy() # Not used + + # Apply boundary conditions + k, ln = NonlinearPiezoSimHb.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 + ).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: + 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 + k = self.k.copy() + c = self.c.copy() + m = self.m.copy() + l = self.l.copy() + + # Create the harmonic balancing matrices from the FEM matrices + harmonic_order = self.harmonic_order + nabla_hb_k = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + nabla_hb = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + for _k in range(harmonic_order): + nabla_hb_k[2*_k:2*_k+2, 2*_k:2*_k+2] = np.array([ + [0, _k], + [-_k, 0] + ]) + nabla_hb[2*_k:2*_k+2, 2*_k:2*_k+2] = np.array([ + [1, 0], + [0, 1] + ]) + + print(k.shape) + + m = sparse.kron(m, np.square(nabla_hb_k), format='lil') + c = sparse.kron(c, nabla_hb_k, format='lil') + k = sparse.kron(k, nabla_hb, format='lil') + + print("Number of nodes:", number_of_nodes) + print("M:", m.shape) + print("C", c.shape) + print("K:", k.shape) + + # Apply boundary conditions + m, c, k = NonlinearPiezoSimHb.apply_dirichlet_bc( + m, + c, + k, + self.dirichlet_nodes + ) + + + f = self.get_load_vector( + harmonic_order, + 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, omega): + l_hb = np.zeros(shape=(2*3*number_of_nodes, 2*3*number_of_nodes)) + l_hb[self.dirichlet_nodes] = 0 + for i in range(3*number_of_nodes): + a = u[2*i] + b = u[2*i+1] + l_hb[2*i:2*i+2] = l[i]*np.array([ + 3/4*(a**3+a*b**2), + 3/4*(a**2*b+b**3) + ]) + return omega**2*m@u+omega*c@u+k@u+l_hb-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, angular_frequency)) + + # 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 = NonlinearPiezoSimHb. \ + calculate_tangent_matrix_hadamard( + current_u, + k, + l + ) + + # Solve for displacement increment + delta_u = linalg.spsolve( + k_tangent, + residual(current_u, angular_frequency) + ) + + # Update step + next_u = current_u - alpha * delta_u + + # Check for convergence + norm = scipy.linalg.norm(residual(next_u, angular_frequency)) + 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, + harmonic_order: float, + 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*harmonic_order, 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( + 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, + sparse.csc_array + ]: + # TODO Parameters are not really ndarrays + # 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.tocsc(), c.tocsc(), k.tocsc() diff --git a/plutho/simulation/nonlinear/piezo_stationary.py b/plutho/simulation/nonlinear/piezo_stationary.py index e68c4a0..68180ba 100644 --- a/plutho/simulation/nonlinear/piezo_stationary.py +++ b/plutho/simulation/nonlinear/piezo_stationary.py @@ -92,7 +92,7 @@ def solve_linear(self): self.u = linalg.spsolve( k, f - ) + ).todense() def solve_newton( self, @@ -148,7 +148,7 @@ def solve_newton( # Residual of the Newton algorithm def residual(u): - return k@u+ln@np.square(u)-load_factor*f + return k@u+ln@np.power(u, self.nonlinear_order)-load_factor*f # Add a best found field parameter which can be returned when the # maximum number of iterations were done From e986d42b0fb1807ad81158160c1723045cb4e92f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Fri, 25 Jul 2025 14:16:24 +0200 Subject: [PATCH 09/17] Add harmonic balancing to nonlinear class --- README.md | 1 + plutho/nonlinear_sim.py | 69 ++++- plutho/simulation/nonlinear/piezo_hb.py | 280 ++++++++++++------ .../simulation/nonlinear/piezo_stationary.py | 6 +- 4 files changed, 255 insertions(+), 101 deletions(-) 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/nonlinear_sim.py b/plutho/nonlinear_sim.py index 95d27f9..453fd29 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -12,6 +12,7 @@ # Local libraries from plutho import Mesh, MaterialManager, FieldType, MaterialData, \ MeshData, SimulationData +from plutho.simulation.nonlinear.piezo_hb import NonlinearPiezoSimHb from .simulation import NonlinearPiezoSimStationary, NonlinearPiezoSimTime, \ NonlinearType @@ -40,7 +41,12 @@ class PiezoNonlinear: sim_directory: str mesh: Mesh material_manager: MaterialManager - solver: Union[None, NonlinearPiezoSimTime, NonlinearPiezoSimStationary] + solver: Union[ + None, + NonlinearPiezoSimTime, + NonlinearPiezoSimStationary, + NonlinearPiezoSimHb + ] boundary_conditions: List[Dict] def __init__( @@ -153,8 +159,21 @@ def clear_dirichlet_bcs(self): self.dirichlet_nodes = [] self.dirichlet_values = [] + def setup_harmonic_balancing(self, harmonic_order: int): + """Sets the simulation type to a nonlinear harmonic balancing + simulation. + + Parameters: + harmonic_order: Number of the harmonics used. + """ + self.solver = NonlinearPiezoSimHb( + self.mesh_data, + self.material_manager, + harmonic_order + ) + def setup_stationary_simulation(self): - """Set the simulation type to a nonlinaer stationary simulation. + """Sets the simulation type to a nonlinaer stationary simulation. """ self.solver = NonlinearPiezoSimStationary( self.mesh_data, @@ -189,7 +208,7 @@ def setup_time_dependent_simulation( ) def set_nonlinearity_order(self, nl_order: int): - if nl_order < 1 or nl_order > 2: + if nl_order < 1 or nl_order > 3: raise NotImplementedError( "Nonlinearity order only implemented for 2 or 3" ) @@ -229,7 +248,6 @@ def set_nonlinearity_type(self, nl_type: NonlinearType, **kwargs): self.nonlinear_type = nl_type self.nonlinear_params = kwargs - def simulate( self, **kwargs @@ -266,6 +284,11 @@ def simulate( - "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. + - "angular_frequency" (float): Can only be set for the harmonic + balance simulation. Sets the angular frequency at which the + simulation is done. + - "skip_assembly" (bool) : Set to true if assembly shall be + skipped. Default is false. """ # Check if materials are set if self.material_manager is None: @@ -286,20 +309,18 @@ def simulate( self.solver.dirichlet_values = np.array(self.dirichlet_values) # Assemble - if "skip_assemble" in kwargs: - if not kwargs["skip_assemble"]: + if "skip_assembly" in kwargs: + skip_assemble = kwargs["skip_assembly"] + del kwargs["skip_assembly"] + else: + skip_assemble = False + + if not skip_assemble: self.solver.assemble( self.nonlinear_order, self.nonlinear_type, **self.nonlinear_params - ) - del kwargs["skip_assemble"] - else: - self.solver.assemble( - self.nonlinear_order, - self.nonlinear_type, - **self.nonlinear_params - ) + ) # Run simulation if isinstance(self.solver, NonlinearPiezoSimTime): @@ -312,6 +333,10 @@ def simulate( self.solver.solve_newton( **kwargs ) + elif isinstance(self.solver, NonlinearPiezoSimHb): + self.solver.solve_newton( + **kwargs + ) else: raise ValueError( "Cannot run simulation for simulation type" @@ -362,6 +387,8 @@ def save_simulation_settings(self, prefix: str = ""): simulation_type = "NonlinearStationary" elif isinstance(self.solver, NonlinearPiezoSimTime): simulation_type = "NonlinearTime" + elif isinstance(self.solver, NonlinearPiezoSimHb): + simulation_type = "NonlinearHarmonicBalance" else: raise ValueError( f"Cannot save simulation type {type(self.solver)}" @@ -395,6 +422,10 @@ def save_simulation_settings(self, prefix: str = ""): # Simulation data if isinstance(self.solver, NonlinearPiezoSimTime): settings["simulation"] = self.solver.simulation_data.__dict__ + elif isinstance(self.solver, NonlinearPiezoSimHb): + settings["simulation"] = { + "harmonic_order": self.solver.harmonic_order + } # Boundary conditions boundary_conditions = {} @@ -490,6 +521,9 @@ def load_simulation_settings( simulation_data.gamma, simulation_data.beta ) + elif simulation_type == "NonlinearHarmonicBalance": + harmonic_order = int(settings["simulation"]["harmonic_order"]) + simulation.setup_harmonic_balancing(harmonic_order) else: raise IOError( f"Cannot deserialize {simulation_type} simulation type" @@ -576,6 +610,13 @@ def load_simulation_results(self): self.solver.u = np.load(u_file) if os.path.isfile(q_file): self.solver.q = np.load(q_file) + elif isinstance(self.solver, NonlinearPiezoSimHb): + # TODO Duplicate with stationary + u_file = os.path.join( + self.sim_directory, + "u.npy" + ) + self.solver.u = np.load(u_file) else: raise ValueError( "Cannot load simulation settings of simulation type", diff --git a/plutho/simulation/nonlinear/piezo_hb.py b/plutho/simulation/nonlinear/piezo_hb.py index 2c306e3..c3eca00 100644 --- a/plutho/simulation/nonlinear/piezo_hb.py +++ b/plutho/simulation/nonlinear/piezo_hb.py @@ -2,7 +2,7 @@ harmonic balancing method.""" # Python standard libraries -from typing import Union, Tuple +from typing import Union, Tuple, Callable # Third party libraries import numpy as np @@ -31,7 +31,7 @@ class NonlinearPiezoSimHb: k: sparse.lil_array c: sparse.lil_array m: sparse.lil_array - l: sparse.lil_array + lu: sparse.lil_array # Resulting fields u: npt.NDArray @@ -40,9 +40,11 @@ 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([]) @@ -50,7 +52,6 @@ def assemble( self, nonlinear_order: int, nonlinear_type: NonlinearType, - harmonic_order: int, **kwargs ): """Redirect to general nonlinear assembly function @@ -61,7 +62,7 @@ def assemble( **kwargs: Parameters for the nonlinear material. """ # Get the default FEM matrices - m, c, k, l = assemble( + m, c, k, lu = assemble( self.mesh_data, self.material_manager, nonlinear_order, @@ -69,24 +70,20 @@ def assemble( **kwargs ) self.nonlinear_order = nonlinear_order - self.harmonic_order = harmonic_order self.m = m self.c = c self.k = k - self.l = l + self.lu = lu def solve_linear(self): """Solves the linear problem Ku=f (ln is not used). """ - # Get FEM matrices - k = self.k.copy() - ln = self.l.copy() # Not used - # Apply boundary conditions - k, ln = NonlinearPiezoSimHb.apply_dirichlet_bc( - k, - ln, + _, _, k = NonlinearPiezoSimHb.apply_dirichlet_bc( + None, + None, + self.k.copy(), self.dirichlet_nodes ) @@ -130,47 +127,27 @@ def solve_newton( number_of_nodes = len(self.mesh_data.nodes) # Get FEM matrices - k = self.k.copy() - c = self.c.copy() - m = self.m.copy() - l = self.l.copy() + 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 - harmonic_order = self.harmonic_order - nabla_hb_k = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - nabla_hb = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - for _k in range(harmonic_order): - nabla_hb_k[2*_k:2*_k+2, 2*_k:2*_k+2] = np.array([ - [0, _k], - [-_k, 0] - ]) - nabla_hb[2*_k:2*_k+2, 2*_k:2*_k+2] = np.array([ - [1, 0], - [0, 1] - ]) - - print(k.shape) - - m = sparse.kron(m, np.square(nabla_hb_k), format='lil') - c = sparse.kron(c, nabla_hb_k, format='lil') - k = sparse.kron(k, nabla_hb, format='lil') - - print("Number of nodes:", number_of_nodes) - print("M:", m.shape) - print("C", c.shape) - print("K:", k.shape) + hb_m, hb_c, hb_k = self.get_harmonic_balancing_matrices( + fem_m, + fem_c, + fem_k + ) # Apply boundary conditions - m, c, k = NonlinearPiezoSimHb.apply_dirichlet_bc( - m, - c, - k, + hb_m, hb_c, hb_k = NonlinearPiezoSimHb.apply_dirichlet_bc( + hb_m, + hb_c, + hb_k, self.dirichlet_nodes ) - f = self.get_load_vector( - harmonic_order, self.dirichlet_nodes, self.dirichlet_values ) @@ -180,63 +157,90 @@ def solve_newton( # If no start value is given calculate one using a linear # simulation current_u = linalg.spsolve( - k.tocsc(), - f + ( + angular_frequency**2*hb_m.tocsc() + + angular_frequency**2*hb_c.tocsc() + + hb_k.tocsc() + ), + load_factor*f ) else: current_u = u_start # Residual of the Newton algorithm - def residual(u, omega): - l_hb = np.zeros(shape=(2*3*number_of_nodes, 2*3*number_of_nodes)) - l_hb[self.dirichlet_nodes] = 0 - for i in range(3*number_of_nodes): - a = u[2*i] - b = u[2*i+1] - l_hb[2*i:2*i+2] = l[i]*np.array([ - 3/4*(a**3+a*b**2), - 3/4*(a**2*b+b**3) - ]) - return omega**2*m@u+omega*c@u+k@u+l_hb-load_factor*f + def residual(u, nonlinear_force): + return ( + angular_frequency**2*hb_m@u + + angular_frequency*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 + ) - # 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, angular_frequency)) + print(nonlinear_force) + + # 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: - return best + print(f"Initial value is already best {best_norm}") + self.u = best_field + return + # Run Newton method for iteration_count in range(max_iter): + print(f"Iteration {iteration_count}") + + 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 + ) + + print(nonlinear_force) + # Calculate tangential stiffness matrix k_tangent = NonlinearPiezoSimHb. \ - calculate_tangent_matrix_hadamard( + calculate_tangent_matrix_numerical( current_u, - k, - l + nonlinear_force, + residual ) # Solve for displacement increment delta_u = linalg.spsolve( k_tangent, - residual(current_u, angular_frequency) + residual(current_u, nonlinear_force) ) # Update step next_u = current_u - alpha * delta_u # Check for convergence - norm = scipy.linalg.norm(residual(next_u, angular_frequency)) + norm = scipy.linalg.norm(residual(next_u, nonlinear_force)) if norm < tolerance: print( - f"Solve Newton found solution after {iteration_count} " + f"Newton found solution after {iteration_count} " f"steps with residual {norm}" ) - return next_u + self.u = next_u + return elif norm < best_norm: best_norm = norm - best = next_u.copy() + best_field = next_u.copy() if iteration_count % 100 == 0 and iteration_count > 0: print("Iteration:", iteration_count) @@ -246,11 +250,83 @@ def residual(u, omega): print("Simulation did not converge") print("Error from best iteration:", best_norm) - self.u = best + + 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 + nonlinear_force = np.zeros(2*3*number_of_nodes) + for i in range(3*number_of_nodes): + for j in range(3*number_of_nodes): + a = u[2*j] + b = u[2*j+1] + nonlinear_force[2*i:2*i+2] += fem_lu[i,j]*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, + ) -> 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: + m: FEM mass matrix. + c: FEM damping matrix. + k: FEM stiffness matrix. + + Returns: + Tuple of HB mass matrix, HB damping matrix and HB stiffness + matrix. + """ + harmonic_order = self.harmonic_order + + # Set harmonic balance derivative matrices + nabla_hb = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + nabla_hb_const = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) + for i in range(harmonic_order): + nabla_hb[2*i:2*i+2, 2*i:2*i+2] = np.array([ + [0, i], + [-i, 0] + ]) + nabla_hb_const[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, np.square(nabla_hb), format='lil') + fem_c = sparse.kron(fem_c, nabla_hb, format='lil') + fem_k = sparse.kron(fem_k, nabla_hb_const, format='lil') + + return fem_m, fem_c, fem_k def get_load_vector( self, - harmonic_order: float, nodes: npt.NDArray, values: npt.NDArray, ) -> npt.NDArray: @@ -268,7 +344,10 @@ def get_load_vector( # Can be initialized to 0 because external load and volume # charge density is 0. - f = np.zeros(3*number_of_nodes*2*harmonic_order, dtype=np.float64) + f = np.zeros( + 3*number_of_nodes*2*self.harmonic_order, + dtype=np.float64 + ) for node, value in zip(nodes, values): f[node] = value @@ -277,13 +356,42 @@ def get_load_vector( # -------- Static functions -------- + @staticmethod + def calculate_tangent_matrix_numerical( + u: npt.NDArray, + nonlinear_force: npt.NDArray, + residual: Callable + ): + """Calculates the tanget stiffness matrix numerically using a first + order finite difference approximation. + + Parameters: + u: Current solution vector. + nonlinear_force: Vector of nonlinear forces. + residual: Residual function. + """ + 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 calculate_tangent_matrix_hadamard( u: npt.NDArray, k: sparse.sparray, - ln: sparse.sparray + lu: sparse.sparray ): - # TODO Duplicate function in piezo_time.py + # TODO This matrix is not valid anymore in harmonic balancing approach """Calculates the tangent matrix based on an analytically expression. @@ -292,7 +400,7 @@ def calculate_tangent_matrix_hadamard( k: FEM stiffness matrix. ln: FEM nonlinear stiffness matrix. """ - k_tangent = k+2*ln@sparse.diags_array(u) + k_tangent = k+3*lu@sparse.diags_array(np.square(u)) return k_tangent @@ -308,17 +416,19 @@ def apply_dirichlet_bc( 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) # 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 + 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 index 68180ba..212fda0 100644 --- a/plutho/simulation/nonlinear/piezo_stationary.py +++ b/plutho/simulation/nonlinear/piezo_stationary.py @@ -157,7 +157,8 @@ def residual(u): # Check if the initial value already suffices the tolerance condition if best_norm < tolerance: - return best + self.u = best + return for iteration_count in range(max_iter): # Calculate tangential stiffness matrix @@ -184,7 +185,8 @@ def residual(u): f"Solve Newton found solution after {iteration_count} " f"steps with residual {norm}" ) - return next_u + self.u = next_u + return elif norm < best_norm: best_norm = norm best = next_u.copy() From 656e48fb11eef0630e66912d6b892179b0fdd41c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Mon, 4 Aug 2025 14:36:19 +0200 Subject: [PATCH 10/17] Add analytical tangent stiffness calculation --- plutho/nonlinear_sim.py | 10 +- plutho/simulation/nonlinear/piezo_hb.py | 200 ++++++++++++++++-------- 2 files changed, 140 insertions(+), 70 deletions(-) diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py index 453fd29..621d63f 100644 --- a/plutho/nonlinear_sim.py +++ b/plutho/nonlinear_sim.py @@ -316,11 +316,11 @@ def simulate( skip_assemble = False if not skip_assemble: - self.solver.assemble( - self.nonlinear_order, - self.nonlinear_type, - **self.nonlinear_params - ) + self.solver.assemble( + self.nonlinear_order, + self.nonlinear_type, + **self.nonlinear_params + ) # Run simulation if isinstance(self.solver, NonlinearPiezoSimTime): diff --git a/plutho/simulation/nonlinear/piezo_hb.py b/plutho/simulation/nonlinear/piezo_hb.py index c3eca00..758ce6e 100644 --- a/plutho/simulation/nonlinear/piezo_hb.py +++ b/plutho/simulation/nonlinear/piezo_hb.py @@ -54,7 +54,7 @@ def assemble( nonlinear_type: NonlinearType, **kwargs ): - """Redirect to general nonlinear assembly function + """Redirect to general nonlinear assembly function. Parameters: nonlinear_order: Order of the nonlinearity. @@ -111,6 +111,8 @@ def solve_newton( 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 @@ -136,7 +138,8 @@ def solve_newton( hb_m, hb_c, hb_k = self.get_harmonic_balancing_matrices( fem_m, fem_c, - fem_k + fem_k, + angular_frequency ) # Apply boundary conditions @@ -158,9 +161,9 @@ def solve_newton( # simulation current_u = linalg.spsolve( ( - angular_frequency**2*hb_m.tocsc() - + angular_frequency**2*hb_c.tocsc() - + hb_k.tocsc() + hb_m + + hb_c + + hb_k ), load_factor*f ) @@ -170,8 +173,8 @@ def solve_newton( # Residual of the Newton algorithm def residual(u, nonlinear_force): return ( - angular_frequency**2*hb_m@u - + angular_frequency*hb_c@u + hb_m@u + + hb_c@u + hb_k@u + nonlinear_force - load_factor*f @@ -184,8 +187,6 @@ def residual(u, nonlinear_force): number_of_nodes ) - print(nonlinear_force) - # 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 @@ -193,15 +194,13 @@ def residual(u, nonlinear_force): 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 + #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): - print(f"Iteration {iteration_count}") - if iteration_count > 0: # Does not need to be updated at step 0 nonlinear_force = self.calculate_nonlinear_force( @@ -210,15 +209,13 @@ def residual(u, nonlinear_force): number_of_nodes ) - print(nonlinear_force) - # Calculate tangential stiffness matrix - k_tangent = NonlinearPiezoSimHb. \ - calculate_tangent_matrix_numerical( + k_tangent = self.calculate_tangent_matrix_analytical( current_u, - nonlinear_force, - residual - ) + hb_m, + hb_c, + hb_k + ) # Solve for displacement increment delta_u = linalg.spsolve( @@ -233,7 +230,7 @@ def residual(u, nonlinear_force): norm = scipy.linalg.norm(residual(next_u, nonlinear_force)) if norm < tolerance: print( - f"Newton found solution after {iteration_count} " + f"Newton found solution after {iteration_count+1} " f"steps with residual {norm}" ) self.u = next_u @@ -270,15 +267,31 @@ def calculate_nonlinear_force( 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): - for j in range(3*number_of_nodes): - a = u[2*j] - b = u[2*j+1] - nonlinear_force[2*i:2*i+2] += fem_lu[i,j]*np.array([ - 3/4*(a**3+a*b**2), - 3/4*(a**2*b+b**3) - ]) + 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 @@ -290,14 +303,15 @@ def get_harmonic_balancing_matrices( 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: - m: FEM mass matrix. - c: FEM damping matrix. - k: FEM stiffness matrix. + 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 @@ -305,23 +319,29 @@ def get_harmonic_balancing_matrices( """ harmonic_order = self.harmonic_order - # Set harmonic balance derivative matrices - nabla_hb = np.zeros(shape=(2*harmonic_order, 2*harmonic_order)) - nabla_hb_const = np.zeros(shape=(2*harmonic_order, 2*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): - nabla_hb[2*i:2*i+2, 2*i:2*i+2] = np.array([ - [0, i], - [-i, 0] + 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_hb_const[2*i:2*i+2, 2*i:2*i+2] = np.array([ + 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, np.square(nabla_hb), format='lil') - fem_c = sparse.kron(fem_c, nabla_hb, format='lil') - fem_k = sparse.kron(fem_k, nabla_hb_const, format='lil') + 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 @@ -354,21 +374,76 @@ def get_load_vector( return f -# -------- Static functions -------- + 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( @@ -385,24 +460,6 @@ def calculate_tangent_matrix_numerical( return k_tangent.tocsc() - @staticmethod - def calculate_tangent_matrix_hadamard( - u: npt.NDArray, - k: sparse.sparray, - lu: sparse.sparray - ): - # TODO This matrix is not valid anymore in harmonic balancing approach - """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+3*lu@sparse.diags_array(np.square(u)) - - return k_tangent @staticmethod def apply_dirichlet_bc( @@ -411,13 +468,26 @@ def apply_dirichlet_bc( k: sparse.lil_array, nodes: npt.NDArray ) -> Tuple[ - sparse.csc_array, sparse.csc_array, sparse.csc_array, sparse.csc_array ]: - # Set rows of matrices to 0 and diagonal of K to 1 (at node points) + """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 From 0ab1ab058f83ec0b81b6ffdcfbc93b7eb00ee84e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Tue, 30 Sep 2025 14:41:48 +0200 Subject: [PATCH 11/17] WIP Nonlinear analytic calculation --- plutho/simulation/nonlinear/base.py | 95 ++++++++++++++++++----------- 1 file changed, 61 insertions(+), 34 deletions(-) diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py index ece0243..af3d260 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -66,7 +66,7 @@ def assemble( nm = kwargs["nonlinear_matrix"] # Check for matrix shape - if len(nm.shape) != nonlinear_order+2: + if len(nm.shape) != nonlinear_order+1: raise ValueError( f"Given nonlinearity matrix shape {nm.shape} does not fit to" f" the given nonlinearity order {nonlinear_order}" @@ -159,7 +159,7 @@ def assemble( ) if nonlinear_type is NonlinearType.Custom: lu_e = ( - integral_u_nonlinear( + integral_u_nonlinear_analytic( current_node_points, nonlinear_matrix, element_order @@ -236,8 +236,51 @@ def assemble( return m, c, k, ln +def integral_u_nonlinear_analytic( + node_points, + nonlinear_elasticity_matrix, + 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) + + double_b = double_b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) + single_b = b_operator_global( + node_points, + jacobian_inverted_t, + s, + t, + element_order + ) -def integral_u_nonlinear( + first = np.einsum( + "ijk,jl,kl->il", + nonlinear_elasticity_matrix, + double_b, + single_b + ) + + second = np.einsum( + "ijk,jl,kl->il", + nonlinear_elasticity_matrix, + single_b, + double_b + ) + + return 1/2*(first+second) + + return quadratic_quadrature(inner, element_order) + +def integral_u_nonlinear_numeric( node_points: npt.NDArray, nonlinear_elasticity_matrix: npt.NDArray, element_order: int @@ -342,6 +385,10 @@ def inner(s, t): element_order ) + print(b_op.shape) + print(bb_op.shape) + print(nonlinear_elasticity_matrix.shape) + left = np.einsum( "ijk,jl,kl->il", nonlinear_elasticity_matrix, @@ -480,6 +527,7 @@ def second_gradient_local_shape_functions_2d( f"order {element_order}" ) + def double_b_operator_global( node_points: npt.NDArray, jacobian_inverted_t: npt.NDArray, @@ -495,7 +543,6 @@ def double_b_operator_global( # Calculate derivatives with respect to global coordinates global_ddn = np.zeros(shape=ddn.shape) - for node_index in range(nodes_per_element): # Get hessian with shape # [ds*ds, ds*dt] @@ -524,39 +571,19 @@ def double_b_operator_global( r = local_to_global_coordinates(node_points, s, t, element_order)[0] + # Calculate B^t{B{N_a}} bb = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) for i in range(nodes_per_element): for j in range(nodes_per_element): - if i == j: - # Diagonal elements have second order derivatives - bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ - [ - global_ddn[0][0][i] + global_ddn[1][1][i]+1/(r**2), - global_ddn[1][0][i] - ], - [ - global_ddn[0][1][i], - global_ddn[0][0][i] + global_ddn[1][1][i] - ] - ] - else: - # Offdiagional elements have 2 first order derivatives - bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ - [ - ( - global_dn[0][i]*global_dn[0][j] - + global_dn[1][i]*global_dn[1][j] - + 1/(r**2) - ), - global_dn[1][i]*global_dn[1][j] - ], - [ - global_dn[0][i]*global_dn[1][j], - ( - global_dn[0][i]*global_dn[0][j] - + global_dn[1][i]*global_dn[1][j] - ) - ] + bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ + [ + global_ddn[0][0][i]+global_ddn[1][1][i]+r**2, + global_ddn[1][0][i], + ], + [ + global_ddn[0][1][i], + global_ddn[1][1][i]+global_ddn[0][1][i] ] + ] return bb From 1a0f66dfe0f32af616c8646daf755e5cd1ff38a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 2 Oct 2025 16:55:25 +0200 Subject: [PATCH 12/17] WIP Reworked nonlinear simulation --- plutho/__init__.py | 2 +- plutho/simulation/base.py | 24 +- plutho/simulation/nonlinear/base.py | 417 +-------------- plutho/simulation/nonlinear/piezo_time.py | 606 +++++++++++++++------- plutho/single_sim.py | 11 +- 5 files changed, 459 insertions(+), 601 deletions(-) diff --git a/plutho/__init__.py b/plutho/__init__.py index 5d4a585..c564bed 100644 --- a/plutho/__init__.py +++ b/plutho/__init__.py @@ -9,4 +9,4 @@ from .single_sim import SingleSimulation, FieldType from .mesh import Mesh from .coupled_sim import CoupledThermoPiezoThermoSim, CoupledFreqPiezoTherm -from .nonlinear_sim import PiezoNonlinear +from .nonlinear_sim import NonlinearPiezoSimTime diff --git a/plutho/simulation/base.py b/plutho/simulation/base.py index 4f4b93f..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 @@ -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) diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py index af3d260..e86a4da 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -2,7 +2,7 @@ # Python standard libraries from enum import Enum -from typing import Tuple, Callable +from typing import Tuple, Callable, Union # Third party libraries import numpy as np @@ -10,7 +10,7 @@ from scipy import sparse # Local libraries -from ..base import MeshData, local_shape_functions_2d, local_to_global_coordinates, b_operator_global, \ +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 @@ -23,77 +23,15 @@ class NonlinearType(Enum): def assemble( mesh_data: MeshData, - material_manager: MaterialManager, - nonlinear_order: int, - nonlinear_type: NonlinearType, - **kwargs + 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 nonlinear type:" - " Rayleigh" - ) - nm = kwargs["nonlinear_matrix"] - - # Check for matrix shape - if len(nm.shape) != nonlinear_order+1: - raise ValueError( - f"Given nonlinearity matrix shape {nm.shape} does not fit to" - f" the given nonlinearity order {nonlinear_order}" - ) - - if nonlinear_order == 3: - raise NotImplementedError( - "Nonlinearity order 3 not implemented yet" - ) - - # Reduce to axisymmetric 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] = nm[i_old, j_old, k_old] - - nonlinear_matrix = nm_reduced - 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 @@ -116,10 +54,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] @@ -157,14 +91,6 @@ def assemble( element_order ) * 2 * np.pi ) - if nonlinear_type is NonlinearType.Custom: - lu_e = ( - integral_u_nonlinear_analytic( - current_node_points, - nonlinear_matrix, - element_order - ) * 2 * np.pi - ) # Now assemble all element matrices for local_p, global_p in enumerate(element): @@ -189,16 +115,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 @@ -208,10 +124,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)) @@ -229,146 +141,30 @@ def assemble( [ku, kuv], [kuv.T, -1*kv] ]).tolil() - ln = sparse.block_array([ - [lu, zeros2x1], - [zeros1x2, zeros1x1] - ]).tolil() - return m, c, k, ln + print("m:", m.count_nonzero()) + print("c:", c.count_nonzero()) + print("k:", k.count_nonzero()) -def integral_u_nonlinear_analytic( - node_points, - nonlinear_elasticity_matrix, - 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) + print("mu:", mu.count_nonzero()) + print("ku:", ku.count_nonzero()) + print("kuv:", kuv.count_nonzero()) + print("kv:", kv.count_nonzero()) - double_b = double_b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - single_b = b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - - first = np.einsum( - "ijk,jl,kl->il", - nonlinear_elasticity_matrix, - double_b, - single_b - ) - - second = np.einsum( - "ijk,jl,kl->il", - nonlinear_elasticity_matrix, - single_b, - double_b - ) - - return 1/2*(first+second) - - return quadratic_quadrature(inner, element_order) - -def integral_u_nonlinear_numeric( - node_points: npt.NDArray, - nonlinear_elasticity_matrix: npt.NDArray, - element_order: int -): - """Calculates the nonlinear stiffness integral + return m, c, k - 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 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) - - def t_i(s, t): - b_op = b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - - return np.einsum( - "ijk,jl,kl->il", - nonlinear_elasticity_matrix, - b_op, - b_op - ) - - t_derived = b_opt_global_t_numerical( - node_points, - t_i, - s, - t, - jacobian_inverted_t, - element_order - ) - - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - - N = np.diag( - local_shape_functions_2d(s, t, element_order) - ) - - return np.outer( - t_derived, - N - ) * 1/2 * r * jacobian_det - - return quadratic_quadrature(inner, element_order) - - -def integral_u_nonlinear_analytic( - 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 nonlinear stiffness 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 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) - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - b_op = b_operator_global( node_points, jacobian_inverted_t, @@ -376,123 +172,23 @@ def inner(s, t): t, element_order ) + r = local_to_global_coordinates(node_points, s, t, element_order)[0] - bb_op = double_b_operator_global( - node_points, - jacobian_inverted_t, - s, - t, - element_order - ) - - print(b_op.shape) - print(bb_op.shape) - print(nonlinear_elasticity_matrix.shape) - - left = np.einsum( - "ijk,jl,kl->il", - nonlinear_elasticity_matrix, - bb_op, - b_op - ) - - right = np.einsum( - "ijk,jl,kl->il", - nonlinear_elasticity_matrix, - b_op, - bb_op - ) + 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 np.outer( - left+right, - local_shape_functions_2d(s, t, element_order) - ) * 1/2 * r * jacobian_det + return B_TCBB * r * jacobian_det return quadratic_quadrature(inner, element_order) -def b_opt_global_t_numerical( - node_points: npt.NDArray, - func: Callable, - s: float, - t: float, - jacobian_inverted_t: npt.NDArray, - element_order: int +def integral_u_nonlinear_cubic( + node_points, + nonlinear_elasticity_tensor, + element_order ): - """Calculates the transposed B operator on the given func numerically (by - using a central finite difference scheme). - The result is evaluated at the given local coordinates s and t. - The derivatives are done with respect to the global coordinates r and z. - - Parameters: - node_points: Node points of the triangle. - func: Function on which the B_T operator is applied. - s: Local coordinate at which the B_T operator is evaluated. - t: Local coordinate at which the B_T operator is evaluated. - jacobian_inverted_T: Inverted and transposed jacobian of the current - triangle, needed for calculating global derivatives. - element_order: Element order of the triangle. - """ - # Idea: - # f contains the function values at the nodes - # first calculate f at the given positions s and t - # then derive f at that point after s and t - # then calculate the derivatives after r and z - - # f should have shape 4 x nodes_per_element - - # TODO Check if this value is suitable - h = 1/100 - - nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) - - # Points at which the function is evaluated - p1 = [s+h,t] - p2 = [s-h,t] - p3 = [s,t+h] - p4 = [s,t-h] - - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - - # Local derivatives - f = func(s, t) - ds_f = 1/(2*h)*(func(*p1)-func(*p2)) - dt_f = 1/(2*h)*(func(*p3)-func(*p4)) - - dr_f = np.zeros(shape=ds_f.shape) - dz_f = np.zeros(shape=dt_f.shape) - - b_num = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) - for i in range(ds_f.shape[0]): - for j in range(ds_f.shape[1]): - # TODO Not all derivatives are needed - # Derivatives with respect to local coordinates s and t - local_der = np.array([ds_f[i, j], dt_f[i, j]]) - - # Derivatives with respect to global coordinates r and z - global_der = np.dot(jacobian_inverted_t, local_der) - dr_f[i, j] = global_der[0] - dz_f[i, j] = global_der[1] - - for node_index in range(nodes_per_element): - current_f = f[:, 2*node_index:2*node_index+2] - dr = dr_f[:, 2*node_index:2*node_index+2] - dz = dz_f[:, 2*node_index:2*node_index+2] - - # Get global derivatives - b_num[2*node_index:2*node_index+2, 2*node_index:2*node_index+2] = \ - [ - [ - dr[0, 0]+dz[2, 1]+current_f[3, 0]/r, - dr[0, 1]+dz[2, 1]+current_f[3, 1]/r - ], - [ - dz[1, 0]+dr[2, 0], - dz[1, 1]+dr[2, 1] - ] - ] - - return b_num + raise NotImplementedError("Cubic nonlinearity not implemented yet") def second_gradient_local_shape_functions_2d( @@ -526,64 +222,3 @@ def second_gradient_local_shape_functions_2d( "Second gradient of shape functions not implemented for element " f"order {element_order}" ) - - -def double_b_operator_global( - node_points: npt.NDArray, - jacobian_inverted_t: npt.NDArray, - s: float, - t: float, - element_order: int -): - nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) - - dn = gradient_local_shape_functions_2d(s, t, element_order) - global_dn = np.dot(jacobian_inverted_t, dn) - ddn = second_gradient_local_shape_functions_2d(s, t, element_order) - - # Calculate derivatives with respect to global coordinates - global_ddn = np.zeros(shape=ddn.shape) - for node_index in range(nodes_per_element): - # Get hessian with shape - # [ds*ds, ds*dt] - # [dt*ds, dt*dt] - # Local hessian since the derivatives are with respect to s and t - local_hessian = np.array([ - [global_ddn[0][0][node_index], global_ddn[0][1][node_index]], - [global_ddn[1][0][node_index], global_ddn[1][1][node_index]] - ]) - - # Global hessian since the derivatives are with respect to r and z - # TODO Verify - global_hessian = np.dot( - np.dot( - jacobian_inverted_t, - local_hessian - ), - jacobian_inverted_t.T - ) - - # TODO Make using slices - global_ddn[0][0][node_index] = global_hessian[0][0] - global_ddn[0][1][node_index] = global_hessian[0][1] - global_ddn[1][0][node_index] = global_hessian[1][0] - global_ddn[1][1][node_index] = global_hessian[1][1] - - r = local_to_global_coordinates(node_points, s, t, element_order)[0] - - # Calculate B^t{B{N_a}} - bb = np.zeros(shape=(2*nodes_per_element, 2*nodes_per_element)) - for i in range(nodes_per_element): - for j in range(nodes_per_element): - bb[2*i:2*(i+1), 2*j:2*(j+1)] = [ - [ - global_ddn[0][0][i]+global_ddn[1][1][i]+r**2, - global_ddn[1][0][i], - ], - [ - global_ddn[0][1][i], - global_ddn[1][1][i]+global_ddn[0][1][i] - ] - ] - - return bb diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py index 201e6a1..5e73f6c 100644 --- a/plutho/simulation/nonlinear/piezo_time.py +++ b/plutho/simulation/nonlinear/piezo_time.py @@ -9,111 +9,258 @@ from scipy.sparse import linalg # Local libraries -from .base import assemble, NonlinearType -from ..base import SimulationData, MeshData +from .base import assemble, NonlinearType, integral_u_nonlinear_quadratic +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 def __init__( self, - mesh_data: MeshData, - material_manager: MaterialManager, - simulation_data: SimulationData + simulation_name: str, + mesh: Mesh ): - 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)) + + 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 + ) - self.dirichlet_nodes = np.array([]) - self.dirichlet_values = np.array([]) + 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, - nonlinear_order: int, nonlinear_type: NonlinearType, - **kwargs + nonlinear_data: Union[float, npt.NDArray] ): """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_order, - nonlinear_type, - **kwargs + self.material_manager ) - self.nonlinear_order = nonlinear_order + self.nonlinear_type = nonlinear_type self.m = m self.c = c self.k = k - self.ln = ln + self.ln = self._assemble_nonlinear_quadratic( + nonlinear_type, + nonlinear_data + ) + + def _assemble_nonlinear_quadratic( + self, + nonlinear_type: NonlinearType, + nonlinear_data: Union[float, npt.NDArray] + ) -> sparse.lil_array: + if nonlinear_type is NonlinearType.Rayleigh: + if not isinstance(nonlinear_data, float): + raise ValueError( + "When setting Rayleigh nonlinearity, a float parameter must" + "be given." + ) + + return nonlinear_data * self.k + + # Custom nonlinearity with full nonlinear tensor + if len(nonlinear_data.shape) != 3: + raise ValueError("A 6x6x6 tensor must be given.") + + # 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] + + nodes = self.mesh_data.nodes + elements = self.mesh_data.elements + element_order = self.mesh_data.element_order + number_of_nodes = len(nodes) + node_points = self.node_points + + 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, + nm_reduced, + 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): + global_i_transformed_1 = \ + 3*number_of_nodes*2*global_k+global_i + global_i_transformed_2 = \ + 3*number_of_nodes*2*(global_k+1)+global_i + + lu[ + 2*global_i_transformed_1:2*global_i_transformed_1+1, + 2*global_j:2*global_j+1 + ] += lu_e[ + 2*local_i:2*local_i+1, + 2*local_j:2*local_j+1, + 0 + ] + lu[ + 2*global_i_transformed_2:2*global_i_transformed_2+1, + 2*global_j:2*global_j+1 + ] += lu_e[ + 2*local_i:2*local_i+1, + 2*local_j:2*local_j+1, + 1 + ] + + ln = sparse.lil_array( + (9*number_of_nodes**2, 3*number_of_nodes), + dtype=np.float64 + ) + + ln[:4*number_of_nodes**2, :2*number_of_nodes] = lu + + return ln - def solve_time_implicit( + 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([]), - u0: Union[npt.NDArray, None] = None + 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 + self.dirichlet_nodes = np.array(self.dirichlet_nodes) + self.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() @@ -148,7 +295,7 @@ 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, ln = NonlinearPiezoSimTime._apply_dirichlet_bc( m, c, k, @@ -156,133 +303,169 @@ def solve_time_implicit( self.dirichlet_nodes ) - if u0 is not None: - u[:, 0] = u0 - - 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] - ) - - # 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 + sparse.save_npz("ln.npz", ln) + + """ + def nonlinear_product(u, L): + result = np.zeros(3*number_of_nodes) # TODO Move outside this fctn + + for k in range(3*number_of_nodes): + for i in range(3*number_of_nodes): + for j in range(3*number_of_nodes): + result[k] += u[i]*L[3*number_of_nodes*k+i, j]*u[j] + + return result + """ + def nonlinear_product(u, L): + 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 + + if u_start is not None: + u[:, 0] = u_start + + try: + 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] ) - # 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 - k_tangent = NonlinearPiezoSimTime. \ - calculate_tangent_matrix_hadamard( - u_i, - k, - ln - ) - delta_u = linalg.spsolve( - (a1*m+a4*c+k_tangent), - residual(u_i, current_u, current_v, current_a, f) + # 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+next_u.T@ln@next_u-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+nonlinear_product(next_u, ln)-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}" - # ) - if norm < tolerance: + # 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 + k_tangent = NonlinearPiezoSimTime. \ + _calculate_tangent_matrix( + u_i, + k, + ln, + number_of_nodes + ) + delta_u = linalg.spsolve( + (a1*m+a4*c+k_tangent), + 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}" + # ) + 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() + self.converged = True + break + elif norm < best_norm: + best_norm = norm + best_u_i = u_i_next.copy() + + if i % 100 == 0 and i > 0: + print("Iteration:", i) + + # Update for next iteration + u_i = u_i_next.copy() + if not self.converged: print( - f"Newton converged at time step {time_index} " - f"after {i+1} iteration(s)" + "Newton did not converge.. Choosing best value: " + f"{best_norm}" ) - # print(u_i_next) - next_u = u_i_next.copy() - self.converged = True - break - elif norm < best_norm: - best_norm = norm - best_u_i = u_i_next.copy() - - if i % 100 == 0 and i > 0: - print("Iteration:", i) - - # Update for next iteration - u_i = u_i_next.copy() - if not self.converged: - print( - "Newton did not converge.. Choosing best value: " - f"{best_norm}" - ) + next_u = best_u_i.copy() + else: + print("Start value norm already below tolerance") next_u = best_u_i.copy() - else: - print("Start value norm already below tolerance") - next_u = best_u_i.copy() - - # 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 (electrode_elements is not None - and electrode_elements.shape[0] > 0): - q[time_index] = calculate_charge( - u[:, time_index+1], - self.material_manager, - electrode_elements, - electrode_normals, - self.mesh_data.nodes, - self.mesh_data.element_order + + # 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 + ) + + except ValueError as e: + print("Caugh't Value Error:", e) + np.save("u.npy", u) + np.save("q.npy", q) + return + self.u = u self.q = q - def get_load_vector( + def _get_load_vector( self, nodes: npt.NDArray, values: npt.NDArray @@ -309,7 +492,7 @@ def get_load_vector( return f @staticmethod - def apply_dirichlet_bc( + def _apply_dirichlet_bc( m: sparse.lil_array, c: sparse.lil_array, k: sparse.lil_array, @@ -338,7 +521,7 @@ def apply_dirichlet_bc( return m.tocsc(), c.tocsc(), k.tocsc(), ln.tocsc() @staticmethod - def calculate_tangent_matrix_hadamard( + def _calculate_tangent_matrix_hadamard( u: npt.NDArray, k: sparse.csc_array, ln: sparse.csc_array @@ -355,3 +538,36 @@ def calculate_tangent_matrix_hadamard( k_tangent = k+2*ln@sparse.diags_array(u) return k_tangent + + @staticmethod + def _calculate_tangent_matrix( + u: npt.NDArray, + k: sparse.csc_array, + ln: sparse.csc_array, + number_of_nodes: int + ): + """ + k_tangent = sparse.csc_array((9*number_of_nodes**2, 3*number_of_nodes)) + + for i in range(3*number_of_nodes): + k_tangent[i*3*number_of_nodes:(i+1)*3*number_of_nodes] = \ + k+2*ln[i*3*number_of_nodes] + k_tangent = k+2*ln@u + """ + k_tangent = sparse.lil_matrix(( + 3*number_of_nodes, 3*number_of_nodes + )) + + k_tangent += k + + for index in range(3*number_of_nodes): + l_k = ln[3*number_of_nodes*index:3*number_of_nodes*(index+1), :] + + l_k_diag = sparse.diags(l_k.diagonal(0)) + left = l_k@u + right = u.T@l_k + left[0] = 0 + right[0] = 0 + k_tangent[index, :] = 2*l_k_diag@u+left+right + + return k_tangent.tocsc() 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 From 78d9c52526a0ab7019ca0bc3cdc416329e45de3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Tue, 7 Oct 2025 11:16:54 +0200 Subject: [PATCH 13/17] Update nonlinear time solver --- plutho/simulation/nonlinear/base.py | 9 - plutho/simulation/nonlinear/piezo_time.py | 375 +++++++++++----------- 2 files changed, 186 insertions(+), 198 deletions(-) diff --git a/plutho/simulation/nonlinear/base.py b/plutho/simulation/nonlinear/base.py index e86a4da..ec6194d 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -142,15 +142,6 @@ def assemble( [kuv.T, -1*kv] ]).tolil() - print("m:", m.count_nonzero()) - print("c:", c.count_nonzero()) - print("k:", k.count_nonzero()) - - print("mu:", mu.count_nonzero()) - print("ku:", ku.count_nonzero()) - print("kuv:", kuv.count_nonzero()) - print("kv:", kv.count_nonzero()) - return m, c, k diff --git a/plutho/simulation/nonlinear/piezo_time.py b/plutho/simulation/nonlinear/piezo_time.py index 5e73f6c..7b1c24b 100644 --- a/plutho/simulation/nonlinear/piezo_time.py +++ b/plutho/simulation/nonlinear/piezo_time.py @@ -5,7 +5,7 @@ import numpy as np import numpy.typing as npt import scipy -from scipy import sparse +from scipy import sparse, integrate, optimize from scipy.sparse import linalg # Local libraries @@ -165,8 +165,8 @@ def _assemble_nonlinear_quadratic( if nonlinear_type is NonlinearType.Rayleigh: if not isinstance(nonlinear_data, float): raise ValueError( - "When setting Rayleigh nonlinearity, a float parameter must" - "be given." + "When setting Rayleigh nonlinearity, a float parameter" + "must be given." ) return nonlinear_data * self.k @@ -207,37 +207,60 @@ def _assemble_nonlinear_quadratic( nm_reduced, 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): - global_i_transformed_1 = \ - 3*number_of_nodes*2*global_k+global_i - global_i_transformed_2 = \ - 3*number_of_nodes*2*(global_k+1)+global_i + # 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_transformed_1:2*global_i_transformed_1+1, - 2*global_j:2*global_j+1 - ] += lu_e[ - 2*local_i:2*local_i+1, - 2*local_j:2*local_j+1, - 0 - ] + ( + 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_transformed_2:2*global_i_transformed_2+1, - 2*global_j:2*global_j+1 - ] += lu_e[ - 2*local_i:2*local_i+1, - 2*local_j:2*local_j+1, - 1 - ] + ( + 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 ) - ln[:4*number_of_nodes**2, :2*number_of_nodes] = lu + 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, : + ] return ln @@ -251,8 +274,8 @@ def simulate( max_iter: int = 300, u_start: Union[npt.NDArray, None] = None ): - self.dirichlet_nodes = np.array(self.dirichlet_nodes) - self.dirichlet_values = np.array(self.dirichlet_values) + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) self.simulation_data = SimulationData( delta_t, @@ -295,27 +318,14 @@ def simulate( 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, ln = self._apply_dirichlet_bc( m, c, k, ln, - self.dirichlet_nodes + dirichlet_nodes ) - sparse.save_npz("ln.npz", ln) - - """ - def nonlinear_product(u, L): - result = np.zeros(3*number_of_nodes) # TODO Move outside this fctn - - for k in range(3*number_of_nodes): - for i in range(3*number_of_nodes): - for j in range(3*number_of_nodes): - result[k] += u[i]*L[3*number_of_nodes*k+i, j]*u[j] - - return result - """ def nonlinear_product(u, L): result = np.zeros(3*number_of_nodes) @@ -329,142 +339,138 @@ def nonlinear_product(u, L): if u_start is not None: u[:, 0] = u_start - try: - 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] + 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] + ) + + # 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+next_u.T@ln@next_u-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+nonlinear_product( + next_u, + ln + )-f ) - # 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+next_u.T@ln@next_u-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+nonlinear_product(next_u, ln)-f - ) - - # 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 - k_tangent = NonlinearPiezoSimTime. \ - _calculate_tangent_matrix( - u_i, - k, - ln, - number_of_nodes - ) - delta_u = linalg.spsolve( - (a1*m+a4*c+k_tangent), - residual(u_i, current_u, current_v, current_a, f) + # 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 + k_tangent = NonlinearPiezoSimTime. \ + _calculate_tangent_matrix( + u_i, + k, + ln, + number_of_nodes ) - u_i_next = u_i - delta_u + delta_u = linalg.spsolve( + (a1*m+a4*c+k_tangent), + 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}" - # ) - 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() - self.converged = True - break - elif norm < best_norm: - best_norm = norm - best_u_i = u_i_next.copy() - - if i % 100 == 0 and i > 0: - print("Iteration:", i) - - # Update for next iteration - u_i = u_i_next.copy() - if not self.converged: + # 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}" + # ) + if norm < tolerance: print( - "Newton did not converge.. Choosing best value: " - f"{best_norm}" + f"Newton converged at time step {time_index} " + f"after {i+1} iteration(s)" ) - next_u = best_u_i.copy() - else: - print("Start value norm already below tolerance") + # print(u_i_next) + next_u = u_i_next.copy() + self.converged = True + break + elif norm < best_norm: + best_norm = norm + best_u_i = u_i_next.copy() + + if i % 100 == 0 and i > 0: + print("Iteration:", i) + + # Update for next iteration + u_i = u_i_next.copy() + if not self.converged: + print( + "Newton did not converge.. Choosing best value: " + f"{best_norm}" + ) next_u = best_u_i.copy() - - # 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 + else: + print("Start value norm already below tolerance") + next_u = best_u_i.copy() + + # 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 ) - # 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 - ) - - except ValueError as e: - print("Caugh't Value Error:", e) - np.save("u.npy", u) - np.save("q.npy", q) - return - self.u = u self.q = q - def _get_load_vector( self, nodes: npt.NDArray, @@ -491,8 +497,8 @@ def _get_load_vector( return f - @staticmethod def _apply_dirichlet_bc( + self, m: sparse.lil_array, c: sparse.lil_array, k: sparse.lil_array, @@ -504,19 +510,18 @@ def _apply_dirichlet_bc( sparse.csc_array, sparse.csc_array ]: + number_of_nodes = len(self.mesh_data.nodes) # TODO Parameters are not really ndarrays # 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 - ln[node, :] = 0 + m[nodes, :] = 0 + c[nodes, :] = 0 + k[nodes, :] = 0 + k[nodes, nodes] = 1 - # Set diagonal values to 1 - k[node, node] = 1 + for index in range(3*number_of_nodes): + ln[nodes+index*3*number_of_nodes, :] = 0 return m.tocsc(), c.tocsc(), k.tocsc(), ln.tocsc() @@ -546,28 +551,20 @@ def _calculate_tangent_matrix( ln: sparse.csc_array, number_of_nodes: int ): - """ - k_tangent = sparse.csc_array((9*number_of_nodes**2, 3*number_of_nodes)) - - for i in range(3*number_of_nodes): - k_tangent[i*3*number_of_nodes:(i+1)*3*number_of_nodes] = \ - k+2*ln[i*3*number_of_nodes] - k_tangent = k+2*ln@u - """ k_tangent = sparse.lil_matrix(( 3*number_of_nodes, 3*number_of_nodes )) k_tangent += k - for index in range(3*number_of_nodes): - l_k = ln[3*number_of_nodes*index:3*number_of_nodes*(index+1), :] + 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, :] - l_k_diag = sparse.diags(l_k.diagonal(0)) - left = l_k@u - right = u.T@l_k - left[0] = 0 - right[0] = 0 - k_tangent[index, :] = 2*l_k_diag@u+left+right + k_tangent += (l_k*u[i]).T + k_tangent += (l_i*u[i]).T - return k_tangent.tocsc() + return k_tangent From 59df74665a480382a63d3d20454b683895416140 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Tue, 7 Oct 2025 11:17:08 +0200 Subject: [PATCH 14/17] WIP Add nonlinear piezo time implicit solver --- .../nonlinear/piezo_time_implicit.py | 498 ++++++++++++++++++ 1 file changed, 498 insertions(+) create mode 100644 plutho/simulation/nonlinear/piezo_time_implicit.py diff --git a/plutho/simulation/nonlinear/piezo_time_implicit.py b/plutho/simulation/nonlinear/piezo_time_implicit.py new file mode 100644 index 0000000..e79751c --- /dev/null +++ b/plutho/simulation/nonlinear/piezo_time_implicit.py @@ -0,0 +1,498 @@ +"""Module for the simulation of nonlinear piezoelectric systems using an +implicit formulation""" + +# Third party libraries +from typing import Tuple, Union +import numpy as np +import numpy.typing as npt +import scipy +from scipy import sparse, integrate, optimize +from scipy.sparse import linalg + +# Local libraries +from .base import assemble, NonlinearType, integral_u_nonlinear_quadratic, \ + integral_m, integral_ku, integral_kuv, integral_kve +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 NonlinearPiezoSimTimeImplicit: + def __init__( + self, + simulation_name: str, + mesh: Mesh + ): + 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)) + + 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_linear( + self, + 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 + ) + + return mu, ku, kuv, kv, cu + + def assemble( + self, + nonlinear_type: NonlinearType, + nonlinear_data: Union[float, npt.NDArray] + ): + """Redirect to general nonlinear assembly function""" + self.material_manager.initialize_materials() + mu, ku, kuv, kv, cu = self.assemble_linear( + self.mesh_data, + self.material_manager + ) + self.nonlinear_type = nonlinear_type + self.ln = self._assemble_nonlinear_quadratic( + nonlinear_type, + nonlinear_data + ) + + self.mu = mu + self.ku = ku + self.kuv = kuv + self.kv = kv + self.cu = cu + + def _assemble_nonlinear_quadratic( + self, + nonlinear_type: NonlinearType, + nonlinear_data: Union[float, npt.NDArray] + ) -> sparse.lil_array: + if nonlinear_type is NonlinearType.Rayleigh: + if not isinstance(nonlinear_data, float): + raise ValueError( + "When setting Rayleigh nonlinearity, a float parameter must" + "be given." + ) + + return nonlinear_data * self.k + + # Custom nonlinearity with full nonlinear tensor + if len(nonlinear_data.shape) != 3: + raise ValueError("A 6x6x6 tensor must be given.") + + # 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] + + nodes = self.mesh_data.nodes + elements = self.mesh_data.elements + element_order = self.mesh_data.element_order + number_of_nodes = len(nodes) + node_points = self.node_points + + 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, + nm_reduced, + 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, : + ] + + return ln + + def simulate( + self, + delta_t, + number_of_time_steps + ): + dirichlet_nodes = np.array(self.dirichlet_nodes) + dirichlet_values = np.array(self.dirichlet_values) + + number_of_nodes = len(self.mesh_data.nodes) + + mu = self.mu.copy() + ku = self.ku.copy() + kuv = self.kuv.copy() + kv = self.kv.coppy() + + alpha_k = self.material_manager.get_alpha_k(0) + alpha_m = self.material_manager.get_alpha_m(0) + + c = alpha_m*mu+alpha_k*ku + ln = self.ln.copy() + + # Apply dirichlet bc to matrices + mu, ku, kuv_u, kuv_phi, kv = _apply_dirichlet_bc_implicit( + mu, + ku, + kuv, + kv, + dirichlet_nodes + ) + + # Init arrays + # Displacement u which is calculated + y_0 = np.zeros(6*number_of_nodes) + + time_interval = (0, number_of_time_steps*delta_t) + time_steps = np.arange(number_of_time_steps)*delta_t + + def nonlinear_product(u, L): + 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 + + def find_nearest(array, value): + return (np.abs(array - value)).argmin() + + def rhs(t, y): + + current_u = y[:3*number_of_nodes] + current_v = y[3*number_of_nodes:] + + time_index = find_nearest(time_steps, t) + + f = self._get_load_vector( + dirichlet_nodes, + dirichlet_values[:, time_index] + ) + eq = f-c@current_v-k@current_u-nonlinear_product(current_u, ln) + + a = linalg.spsolve(m, eq) + + return np.concatenate((current_v, a)) + + sol = integrate.solve_ivp( + rhs, + time_interval, + y_0, + t_eval=time_steps, + method="RK45" + ) + + if sol.status != 0: + print("Something went wrong while simulating") + print(sol.message) + + print(sol) + + self.u = sol.y[:3*number_of_nodes, :] + + def _apply_dirichlet_bc_implicit( + self, + mu: sparse.lil_array, + ku: sparse.lil_array, + kuv: sparse.lil_array, + kv: sparse.lil_array, + cu: sparse.lil_array, + nodes: npt.NDArray + ) -> Tuple[ + sparse.csc_array, + sparse.csc_array, + sparse.csc_array, + sparse.csc_array, + sparse.csc_array + ]: + number_of_nodes = len(self.mesh_data.nodes) + nodes_u = [] + nodes_phi = [] + for node in nodes: + if node > 2*number_of_nodes: + nodes_phi.append(node) + else: + nodes_u.append(nodes) + + kuv_u = kuv.copy() + kuv_phi = kuv.copy() + + # For u field matrices + mu[nodes_u, :] = 0 + kuv_u[nodes_u, :] = 0 + ku[nodes, :] = 0 + ku[nodes, nodes] = 1 + + # For phi field matrices + kuv_phi[nodes, :] = 0 + kv[nodes, :] = 0 + kv[nodes, nodes] = 1 + + return mu.tocsc(), ku.tocsc(), kuv_u.tocsc(), kuv_phi.tocsc(), \ + kv.tocsc() From 8f56987eb1b09ba2d57669072eb2de8e2c2e9e1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 23 Oct 2025 14:48:50 +0200 Subject: [PATCH 15/17] Update nonlinaer simulations Add Nonlinearity class to handle different nonlinear types. Update nonlinear example script. Remove stationary nonlinear simulation. --- plutho/__init__.py | 4 +- plutho/nonlinear_sim.py | 624 ------------------ plutho/simulation/__init__.py | 3 +- plutho/simulation/nonlinear/__init__.py | 7 +- plutho/simulation/nonlinear/base.py | 281 ++++++-- plutho/simulation/nonlinear/piezo_hb.py | 7 +- .../simulation/nonlinear/piezo_stationary.py | 273 -------- plutho/simulation/nonlinear/piezo_time.py | 262 ++------ .../nonlinear/piezo_time_implicit.py | 12 +- plutho/simulation/piezo_freq.py | 24 + pyproject.toml | 2 +- scripts/example_nonlinear.py | 199 +++--- 12 files changed, 425 insertions(+), 1273 deletions(-) delete mode 100644 plutho/nonlinear_sim.py delete mode 100644 plutho/simulation/nonlinear/piezo_stationary.py diff --git a/plutho/__init__.py b/plutho/__init__.py index c564bed..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, NonlinearPiezoSimHb + 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 NonlinearPiezoSimTime diff --git a/plutho/nonlinear_sim.py b/plutho/nonlinear_sim.py deleted file mode 100644 index 621d63f..0000000 --- a/plutho/nonlinear_sim.py +++ /dev/null @@ -1,624 +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 plutho.simulation.nonlinear.piezo_hb import NonlinearPiezoSimHb -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, - NonlinearPiezoSimHb - ] - 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.solver = None - self.nonlinear_material_matrix = None - - 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 - - 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_harmonic_balancing(self, harmonic_order: int): - """Sets the simulation type to a nonlinear harmonic balancing - simulation. - - Parameters: - harmonic_order: Number of the harmonics used. - """ - self.solver = NonlinearPiezoSimHb( - self.mesh_data, - self.material_manager, - harmonic_order - ) - - def setup_stationary_simulation(self): - """Sets 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_order(self, nl_order: int): - if nl_order < 1 or nl_order > 3: - raise NotImplementedError( - "Nonlinearity order only implemented for 2 or 3" - ) - - self.nonlinear_order = nl_order - - def set_nonlinearity_type(self, nl_type: 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 nl_type is NonlinearType.Rayleigh: - if "zeta" not in kwargs: - raise ValueError( - "Missing 'zeta' parameter for nonlinear type: Rayleigh" - ) - elif nl_type 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 = nl_type - 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. - - "angular_frequency" (float): Can only be set for the harmonic - balance simulation. Sets the angular frequency at which the - simulation is done. - - "skip_assembly" (bool) : Set to true if assembly shall be - skipped. Default is false. - """ - # 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 - if "skip_assembly" in kwargs: - skip_assemble = kwargs["skip_assembly"] - del kwargs["skip_assembly"] - else: - skip_assemble = False - - if not skip_assemble: - self.solver.assemble( - self.nonlinear_order, - self.nonlinear_type, - **self.nonlinear_params - ) - - # 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 - ) - elif isinstance(self.solver, NonlinearPiezoSimHb): - 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" - elif isinstance(self.solver, NonlinearPiezoSimHb): - simulation_type = "NonlinearHarmonicBalance" - 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, - "element_order": self.mesh_data.element_order, - "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() - nonlinear_materials = {} - nonlinear_materials["nonlinear_type"] = self.nonlinear_type.value - - match self.nonlinear_type: - case NonlinearType.Rayleigh: - nonlinear_materials["zeta"] = self.nonlinear_params["zeta"] - case NonlinearType.Custom: - nonlinear_materials["nonlinear_matrix"] = \ - self.nonlinear_params["nonlinear_matrix"].tolist() - - material_settings["nonlinear"] = nonlinear_materials - - # Simulation data - if isinstance(self.solver, NonlinearPiezoSimTime): - settings["simulation"] = self.solver.simulation_data.__dict__ - elif isinstance(self.solver, NonlinearPiezoSimHb): - settings["simulation"] = { - "harmonic_order": self.solver.harmonic_order - } - - # 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"] - element_order = int(settings["general"]["element_order"]) - simulation_type = settings["general"]["simulation_type"] - simulation = PiezoNonlinear( - simulation_folder, - "", - Mesh(mesh_file, element_order) - ) - # 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 - ) - elif simulation_type == "NonlinearHarmonicBalance": - harmonic_order = int(settings["simulation"]["harmonic_order"]) - simulation.setup_harmonic_balancing(harmonic_order) - 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) - - # Load nonlinearities - nonlinear_material = material_settings["nonlinear"] - nonlinear_type = NonlinearType( - nonlinear_material["nonlinear_type"] - ) - - match nonlinear_type: - case NonlinearType.Rayleigh: - nonlinear_params = { - "zeta": float(nonlinear_material["zeta"]) - } - case NonlinearType.Custom: - nonlinear_params = { - "nonlinear_matrix": np.array( - nonlinear_material["nonlinear_matrix"] - ) - } - - # Load other materials - for material_name in material_settings.keys(): - if material_name == "nonlinear": - continue - - simulation.add_material( - material_name, - MaterialData( - **material_settings[material_name]["material_data"] - ), - material_settings[material_name]["physical_group_name"] - ) - - simulation.material_manager.initialize_materials() - - simulation.set_nonlinearity_type( - nonlinear_type, - **nonlinear_params - ) - - # 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) - elif isinstance(self.solver, NonlinearPiezoSimHb): - # TODO Duplicate with stationary - u_file = os.path.join( - self.sim_directory, - "u.npy" - ) - self.solver.u = np.load(u_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 d031cc9..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, NonlinearPiezoSimHb +from .nonlinear import NLPiezoTime, NonlinearType, NLPiezoHB, Nonlinearity diff --git a/plutho/simulation/nonlinear/__init__.py b/plutho/simulation/nonlinear/__init__.py index 8ad9c91..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 .piezo_hb import NonlinearPiezoSimHb +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 ec6194d..3351c40 100644 --- a/plutho/simulation/nonlinear/base.py +++ b/plutho/simulation/nonlinear/base.py @@ -2,11 +2,10 @@ # Python standard libraries from enum import Enum -from typing import Tuple, Callable, Union +from typing import Tuple # Third party libraries import numpy as np -import numpy.typing as npt from scipy import sparse # Local libraries @@ -15,12 +14,18 @@ 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 @@ -173,43 +178,233 @@ def inner(s, t): return quadratic_quadrature(inner, element_order) +# +# -------- Classes -------- +# +class Nonlinearity: -def integral_u_nonlinear_cubic( - node_points, - nonlinear_elasticity_tensor, - element_order -): - raise NotImplementedError("Cubic nonlinearity not implemented yet") + 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 second_gradient_local_shape_functions_2d( - s, - t, - element_order -) -> npt.NDArray: - # Since the second gradient has to be done for s and t as well, 3 indices - # are needed [i,j,k]. The first two represent if its derived after s [i=0] - # or t [i=1] and the last index [k] represents the node index of the - # element. - nodes_per_element = int(1/2*(element_order+1)*(element_order+2)) - match element_order: - case 1: - return np.zeros(shape=(2, 2, nodes_per_element)) - case 2: - return np.array([ - [ # (outer) d_s - [4, 4, 0, -8, 0, 0], # (inner) d_s - [4, 0, 0, -4, 4, -4] # (inner) d_t - ], - [ # (outer) d_t - [4, 0, 0, -4, 4, -4], # (inner) d_s - [4, 0, 4, 0, 0, -8] # (inner) d_t - ], - ]) - case 3: - raise NotImplementedError("Not yet implemented for order 3") - - raise ValueError( - "Second gradient of shape functions not implemented for element " - f"order {element_order}" - ) + 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 index 758ce6e..6c098dd 100644 --- a/plutho/simulation/nonlinear/piezo_hb.py +++ b/plutho/simulation/nonlinear/piezo_hb.py @@ -16,7 +16,8 @@ from ..base import MeshData from plutho.materials import MaterialManager -class NonlinearPiezoSimHb: + +class NLPiezoHB: """Implementes a nonlinear FEM harmonic balancing simulation. """ # Simulation parameters @@ -80,7 +81,7 @@ def solve_linear(self): """Solves the linear problem Ku=f (ln is not used). """ # Apply boundary conditions - _, _, k = NonlinearPiezoSimHb.apply_dirichlet_bc( + _, _, k = NLPiezoHB.apply_dirichlet_bc( None, None, self.k.copy(), @@ -143,7 +144,7 @@ def solve_newton( ) # Apply boundary conditions - hb_m, hb_c, hb_k = NonlinearPiezoSimHb.apply_dirichlet_bc( + hb_m, hb_c, hb_k = NLPiezoHB.apply_dirichlet_bc( hb_m, hb_c, hb_k, diff --git a/plutho/simulation/nonlinear/piezo_stationary.py b/plutho/simulation/nonlinear/piezo_stationary.py deleted file mode 100644 index 212fda0..0000000 --- a/plutho/simulation/nonlinear/piezo_stationary.py +++ /dev/null @@ -1,273 +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_order: int, - nonlinear_type: NonlinearType, - **kwargs - ): - """Redirect to general nonlinear assembly function""" - _, _, k, ln = assemble( - self.mesh_data, - self.material_manager, - nonlinear_order, - nonlinear_type, - **kwargs - ) - self.nonlinear_order = nonlinear_order - 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 - ).todense() - - 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.power(u, self.nonlinear_order)-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: - self.u = best - return - - 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}" - ) - self.u = next_u - return - 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 7b1c24b..a4703fc 100644 --- a/plutho/simulation/nonlinear/piezo_time.py +++ b/plutho/simulation/nonlinear/piezo_time.py @@ -5,11 +5,11 @@ import numpy as np import numpy.typing as npt import scipy -from scipy import sparse, integrate, optimize +from scipy import sparse from scipy.sparse import linalg # Local libraries -from .base import assemble, NonlinearType, integral_u_nonlinear_quadratic +from .base import assemble, Nonlinearity from ..base import SimulationData, MeshData, MaterialData, FieldType, \ create_node_points from plutho.simulation.piezo_time import calculate_charge @@ -17,12 +17,13 @@ from plutho.mesh.mesh import Mesh -class NonlinearPiezoSimTime: +class NLPiezoTime: def __init__( self, simulation_name: str, - mesh: Mesh + mesh: Mesh, + nonlinearity: Nonlinearity ): self.simulation_name = simulation_name @@ -42,6 +43,8 @@ def __init__( ) self.material_manager = MaterialManager(len(elements)) + self.nonlinearity = nonlinearity + self.nonlinearity.set_mesh_data(self.mesh_data, self.node_points) def add_material( self, @@ -137,132 +140,18 @@ def set_electrode(self, electrode_name: str): (len(self.electrode_elements), 1) ) - def assemble( - self, - nonlinear_type: NonlinearType, - nonlinear_data: Union[float, npt.NDArray] - ): + 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.nonlinear_type = nonlinear_type self.m = m self.c = c self.k = k - self.ln = self._assemble_nonlinear_quadratic( - nonlinear_type, - nonlinear_data - ) - - def _assemble_nonlinear_quadratic( - self, - nonlinear_type: NonlinearType, - nonlinear_data: Union[float, npt.NDArray] - ) -> sparse.lil_array: - if nonlinear_type is NonlinearType.Rayleigh: - if not isinstance(nonlinear_data, float): - raise ValueError( - "When setting Rayleigh nonlinearity, a float parameter" - "must be given." - ) - - return nonlinear_data * self.k - - # Custom nonlinearity with full nonlinear tensor - if len(nonlinear_data.shape) != 3: - raise ValueError("A 6x6x6 tensor must be given.") - - # 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] - - nodes = self.mesh_data.nodes - elements = self.mesh_data.elements - element_order = self.mesh_data.element_order - number_of_nodes = len(nodes) - node_points = self.node_points - - 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, - nm_reduced, - 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, : - ] - return ln + self.nonlinearity.assemble(m, c, k) def simulate( self, @@ -288,7 +177,6 @@ def simulate( 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)) @@ -318,23 +206,26 @@ def simulate( q = np.zeros(number_of_time_steps, dtype=np.float64) # Apply dirichlet bc to matrices - m, c, k, ln = self._apply_dirichlet_bc( + m, c, k = self._apply_dirichlet_bc( m, c, k, - ln, dirichlet_nodes ) - - def nonlinear_product(u, L): - 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 + 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 @@ -347,22 +238,6 @@ def nonlinear_product(u, L): 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+next_u.T@ln@next_u-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+nonlinear_product( - next_u, - ln - )-f - ) - # Values of the current time step current_u = u[:, time_index] current_v = v[:, time_index] @@ -391,55 +266,50 @@ def residual(next_u, current_u, v, a, f): # Start newton iterations for i in range(max_iter): # Calculate tangential stiffness matrix - k_tangent = NonlinearPiezoSimTime. \ - _calculate_tangent_matrix( - u_i, - k, - ln, - number_of_nodes - ) + 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() + 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() + 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] = ( @@ -502,34 +372,26 @@ def _apply_dirichlet_bc( 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 ]: - number_of_nodes = len(self.mesh_data.nodes) - # TODO Parameters are not really ndarrays # Set rows of matrices to 0 and diagonal of K to 1 (at node points) - - # Matrices for u_r component m[nodes, :] = 0 c[nodes, :] = 0 k[nodes, :] = 0 k[nodes, nodes] = 1 - for index in range(3*number_of_nodes): - ln[nodes+index*3*number_of_nodes, :] = 0 - - return m.tocsc(), c.tocsc(), k.tocsc(), ln.tocsc() + return m.tocsc(), c.tocsc(), k.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 @@ -540,31 +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) - - return k_tangent - - @staticmethod - def _calculate_tangent_matrix( - u: npt.NDArray, - k: sparse.csc_array, - ln: sparse.csc_array, - number_of_nodes: int - ): - 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 + linear = k + nonlinear = self.nonlinearity.evaluate_jacobian( + u, m, c, k + ) - return k_tangent + return (linear + nonlinear).tocsc() diff --git a/plutho/simulation/nonlinear/piezo_time_implicit.py b/plutho/simulation/nonlinear/piezo_time_implicit.py index e79751c..c3a12c5 100644 --- a/plutho/simulation/nonlinear/piezo_time_implicit.py +++ b/plutho/simulation/nonlinear/piezo_time_implicit.py @@ -5,20 +5,18 @@ from typing import Tuple, Union import numpy as np import numpy.typing as npt -import scipy -from scipy import sparse, integrate, optimize +from scipy import sparse, integrate from scipy.sparse import linalg # Local libraries -from .base import assemble, NonlinearType, integral_u_nonlinear_quadratic, \ - integral_m, integral_ku, integral_kuv, integral_kve -from ..base import SimulationData, MeshData, MaterialData, FieldType, \ +from .base import integral_u_nonlinear_quadratic, \ + integral_m, integral_ku, integral_kuv, integral_kve, NonlinearType +from ..base import 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 NonlinearPiezoSimTimeImplicit: +class NLPiezoTimeImplicit: def __init__( self, simulation_name: str, diff --git a/plutho/simulation/piezo_freq.py b/plutho/simulation/piezo_freq.py index 24f17f3..e686641 100644 --- a/plutho/simulation/piezo_freq.py +++ b/plutho/simulation/piezo_freq.py @@ -193,6 +193,9 @@ def assemble(self): * 2 * np.pi ) + if 3 in element and 2 in element: + print("Zero element") + # Now assemble all element matrices for local_p, global_p in enumerate(element): for local_q, global_q in enumerate(element): @@ -244,6 +247,15 @@ def assemble(self): self.c = c.tolil() self.k = k.tolil() + content = "" + for i in range(3*number_of_nodes): + for j in range(3*number_of_nodes): + content += f"{self.k[i,j].real} {self.k[i,j].imag} " + content += "\n" + + with open("k_py.txt", "w") as fd: + fd.write(content) + def solve_frequency( self, electrode_elements: npt.NDArray, @@ -306,6 +318,18 @@ def solve_frequency( + k ) + """ + content = "" + for i in range(3*number_of_nodes): + for j in range(3*number_of_nodes): + val = a[i,j] + content += f"{val.real} {val.imag} " + content += "\n" + + with open("a.txt", "w") as fd: + fd.write(content) + """ + f = self.get_load_vector( self.dirichlet_nodes, self.dirichlet_values, 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) From 35d516de67ea3f80ddca17a2d773bcdb67c1477d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 23 Oct 2025 15:21:06 +0200 Subject: [PATCH 16/17] Remove useless code --- plutho/simulation/piezo_freq.py | 24 ------------------------ tests/test_fields.py | 4 ++-- 2 files changed, 2 insertions(+), 26 deletions(-) diff --git a/plutho/simulation/piezo_freq.py b/plutho/simulation/piezo_freq.py index e686641..deeacfc 100644 --- a/plutho/simulation/piezo_freq.py +++ b/plutho/simulation/piezo_freq.py @@ -193,9 +193,6 @@ def assemble(self): * 2 * np.pi ) - if 3 in element and 2 in element: - print("Zero element") - # Now assemble all element matrices for local_p, global_p in enumerate(element): for local_q, global_q in enumerate(element): @@ -247,15 +244,6 @@ def assemble(self): self.c = c.tolil() self.k = k.tolil() - content = "" - for i in range(3*number_of_nodes): - for j in range(3*number_of_nodes): - content += f"{self.k[i,j].real} {self.k[i,j].imag} " - content += "\n" - - with open("k_py.txt", "w") as fd: - fd.write(content) - def solve_frequency( self, electrode_elements: npt.NDArray, @@ -317,18 +305,6 @@ def solve_frequency( + 1j*angular_frequency*c + k ) - - """ - content = "" - for i in range(3*number_of_nodes): - for j in range(3*number_of_nodes): - val = a[i,j] - content += f"{val.real} {val.imag} " - content += "\n" - - with open("a.txt", "w") as fd: - fd.write(content) - """ f = self.get_load_vector( self.dirichlet_nodes, 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( **{ From 0e5f6637436c9ac222bab53923884a340db04018 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6lscher?= Date: Thu, 23 Oct 2025 15:28:17 +0200 Subject: [PATCH 17/17] Remove piezo_time_implicit.py --- .../nonlinear/piezo_time_implicit.py | 496 ------------------ 1 file changed, 496 deletions(-) delete mode 100644 plutho/simulation/nonlinear/piezo_time_implicit.py diff --git a/plutho/simulation/nonlinear/piezo_time_implicit.py b/plutho/simulation/nonlinear/piezo_time_implicit.py deleted file mode 100644 index c3a12c5..0000000 --- a/plutho/simulation/nonlinear/piezo_time_implicit.py +++ /dev/null @@ -1,496 +0,0 @@ -"""Module for the simulation of nonlinear piezoelectric systems using an -implicit formulation""" - -# Third party libraries -from typing import Tuple, Union -import numpy as np -import numpy.typing as npt -from scipy import sparse, integrate -from scipy.sparse import linalg - -# Local libraries -from .base import integral_u_nonlinear_quadratic, \ - integral_m, integral_ku, integral_kuv, integral_kve, NonlinearType -from ..base import MeshData, MaterialData, FieldType, \ - create_node_points -from plutho.materials import MaterialManager -from plutho.mesh.mesh import Mesh - -class NLPiezoTimeImplicit: - def __init__( - self, - simulation_name: str, - mesh: Mesh - ): - 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)) - - 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_linear( - self, - 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 - ) - - return mu, ku, kuv, kv, cu - - def assemble( - self, - nonlinear_type: NonlinearType, - nonlinear_data: Union[float, npt.NDArray] - ): - """Redirect to general nonlinear assembly function""" - self.material_manager.initialize_materials() - mu, ku, kuv, kv, cu = self.assemble_linear( - self.mesh_data, - self.material_manager - ) - self.nonlinear_type = nonlinear_type - self.ln = self._assemble_nonlinear_quadratic( - nonlinear_type, - nonlinear_data - ) - - self.mu = mu - self.ku = ku - self.kuv = kuv - self.kv = kv - self.cu = cu - - def _assemble_nonlinear_quadratic( - self, - nonlinear_type: NonlinearType, - nonlinear_data: Union[float, npt.NDArray] - ) -> sparse.lil_array: - if nonlinear_type is NonlinearType.Rayleigh: - if not isinstance(nonlinear_data, float): - raise ValueError( - "When setting Rayleigh nonlinearity, a float parameter must" - "be given." - ) - - return nonlinear_data * self.k - - # Custom nonlinearity with full nonlinear tensor - if len(nonlinear_data.shape) != 3: - raise ValueError("A 6x6x6 tensor must be given.") - - # 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] - - nodes = self.mesh_data.nodes - elements = self.mesh_data.elements - element_order = self.mesh_data.element_order - number_of_nodes = len(nodes) - node_points = self.node_points - - 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, - nm_reduced, - 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, : - ] - - return ln - - def simulate( - self, - delta_t, - number_of_time_steps - ): - dirichlet_nodes = np.array(self.dirichlet_nodes) - dirichlet_values = np.array(self.dirichlet_values) - - number_of_nodes = len(self.mesh_data.nodes) - - mu = self.mu.copy() - ku = self.ku.copy() - kuv = self.kuv.copy() - kv = self.kv.coppy() - - alpha_k = self.material_manager.get_alpha_k(0) - alpha_m = self.material_manager.get_alpha_m(0) - - c = alpha_m*mu+alpha_k*ku - ln = self.ln.copy() - - # Apply dirichlet bc to matrices - mu, ku, kuv_u, kuv_phi, kv = _apply_dirichlet_bc_implicit( - mu, - ku, - kuv, - kv, - dirichlet_nodes - ) - - # Init arrays - # Displacement u which is calculated - y_0 = np.zeros(6*number_of_nodes) - - time_interval = (0, number_of_time_steps*delta_t) - time_steps = np.arange(number_of_time_steps)*delta_t - - def nonlinear_product(u, L): - 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 - - def find_nearest(array, value): - return (np.abs(array - value)).argmin() - - def rhs(t, y): - - current_u = y[:3*number_of_nodes] - current_v = y[3*number_of_nodes:] - - time_index = find_nearest(time_steps, t) - - f = self._get_load_vector( - dirichlet_nodes, - dirichlet_values[:, time_index] - ) - eq = f-c@current_v-k@current_u-nonlinear_product(current_u, ln) - - a = linalg.spsolve(m, eq) - - return np.concatenate((current_v, a)) - - sol = integrate.solve_ivp( - rhs, - time_interval, - y_0, - t_eval=time_steps, - method="RK45" - ) - - if sol.status != 0: - print("Something went wrong while simulating") - print(sol.message) - - print(sol) - - self.u = sol.y[:3*number_of_nodes, :] - - def _apply_dirichlet_bc_implicit( - self, - mu: sparse.lil_array, - ku: sparse.lil_array, - kuv: sparse.lil_array, - kv: sparse.lil_array, - cu: sparse.lil_array, - nodes: npt.NDArray - ) -> Tuple[ - sparse.csc_array, - sparse.csc_array, - sparse.csc_array, - sparse.csc_array, - sparse.csc_array - ]: - number_of_nodes = len(self.mesh_data.nodes) - nodes_u = [] - nodes_phi = [] - for node in nodes: - if node > 2*number_of_nodes: - nodes_phi.append(node) - else: - nodes_u.append(nodes) - - kuv_u = kuv.copy() - kuv_phi = kuv.copy() - - # For u field matrices - mu[nodes_u, :] = 0 - kuv_u[nodes_u, :] = 0 - ku[nodes, :] = 0 - ku[nodes, nodes] = 1 - - # For phi field matrices - kuv_phi[nodes, :] = 0 - kv[nodes, :] = 0 - kv[nodes, nodes] = 1 - - return mu.tocsc(), ku.tocsc(), kuv_u.tocsc(), kuv_phi.tocsc(), \ - kv.tocsc()