diff --git a/docs/source/index.rst b/docs/source/index.rst index efa0d20..71b969c 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -82,3 +82,9 @@ Consult our :ref:`Roadmap ` for planned feature implementations. .. note:: This project is still in early development and documentation may be incomplete. Version 1.0 is scheduled for release in early 2026. Stay tuned! + + +.. toctree:: + :hidden: + + references diff --git a/docs/source/references.rst b/docs/source/references.rst new file mode 100644 index 0000000..0883887 --- /dev/null +++ b/docs/source/references.rst @@ -0,0 +1,5 @@ +References +========== + +.. bibliography:: ./refs.bib + :style: unsrt \ No newline at end of file diff --git a/docs/source/tutorials/channel_flow.rst b/docs/source/tutorials/channel_flow.rst index 7794725..14e1b75 100644 --- a/docs/source/tutorials/channel_flow.rst +++ b/docs/source/tutorials/channel_flow.rst @@ -171,9 +171,3 @@ Lizzy will save the following fields: "FillFactor", "FreeSurface", "Pressure", " :width: 60% :align: center - -References ----------- - -.. bibliography:: ../refs.bib - :style: unsrt \ No newline at end of file diff --git a/docs/source/user_guide/assigning_boundary_conditions.rst b/docs/source/user_guide/assigning_boundary_conditions.rst index ac9a8b2..6010baf 100644 --- a/docs/source/user_guide/assigning_boundary_conditions.rst +++ b/docs/source/user_guide/assigning_boundary_conditions.rst @@ -4,11 +4,11 @@ Boundary conditions operations In this section we look at various boundary condition operations, including creating boundary conditions and modifying them at runtime. All operations can be performed using the :class:`~lizzy.LizzyModel` user-facing methods. - Creation and assignment of inlets ---------------------------------- -The following operations are to be performed **before** the solver is initialised by calling :meth:`~lizzy.LizzyModel.initialise_solver`. +.. note:: + The following operations are to be performed **before** the solver is initialised by calling :meth:`~lizzy.LizzyModel.initialise_solver`. Creating an inlet ~~~~~~~~~~~~~~~~~~ @@ -29,13 +29,13 @@ Let's say we have imported a mesh with one boundary named "left_edge" and one bo Physical domains: ['domain'] Physical lines: ['left_edge', 'right_edge'] -First, we will create an inlet using the :meth:`~lizzy.LizzyModel.create_inlet` method. This method requires two arguments: the inlet initial pressure (in Pa) and a unique string identifier for the inlet. Let's create an inlet with a pressure of 1.0E05 Pa tagged as "inlet_1": +First, we will create an inlet using the :meth:`~lizzy.LizzyModel.create_inlet` method. This method requires two arguments: a unique name for the inlet and the inlet initial pressure (in Pa). Let's create an inlet named "inlet_1" with a pressure of 1.0E05 Pa: .. code-block:: - model.create_inlet(1e5, "inlet_1") + model.create_inlet("inlet_1", 1e5) -An :class:`~lizzy.gates.Inlet` object is created and stored in the model, but it is not assigned yet. To assign the inlet to a specific boundary, we use the :meth:`~lizzy.LizzyModel.assign_inlet` method, providing the name of the inlet and the name of the mesh boundary where we want to assign it: +An :class:`~lizzy.gates.Inlet` object is created and stored in the model, but it is not assigned yet to any boundary. To do so, we use the :meth:`~lizzy.LizzyModel.assign_inlet` method, providing the name of the inlet and the name of the mesh boundary where we want to assign it: .. code-block:: @@ -47,14 +47,14 @@ An :class:`~lizzy.gates.Inlet` object is created and stored in the model, but it .. code-block:: - new_inlet = model.create_inlet(1e5, "inlet_1") + new_inlet = model.create_inlet("inlet_1", 1e5) model.assign_inlet(new_inlet, "left_edge") Once an inlet is assigned, it is set to "open" state by default. We can check its state at any time using the :attr:`~lizzy.gates.Inlet.is_open` property (read-only). .. code-block:: console - >>> new_inlet = model.create_inlet(1e5, "inlet_1") + >>> new_inlet = model.create_inlet("inlet_1", 1e5) >>> model.assign_inlet(new_inlet, "left_edge") >>> new_inlet.is_open @@ -78,8 +78,8 @@ This expression returns the :class:`~lizzy.gates.Inlet` object with that name. W Runtime operations ------------------- -The following operations can be performed at any time **after** the solver has been initialised by calling :meth:`~lizzy.LizzyModel.initialise_solver`. - +.. note:: + The following operations can be performed at any time **after** the solver has been initialised by calling :meth:`~lizzy.LizzyModel.initialise_solver`. Opening / closing inlets ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -141,7 +141,6 @@ This will produce the same effect as using the :meth:`~lizzy.LizzyModel.open_inl >>> Fatal error: The application has terminated. >>> No inlets are currently open. At least one inlet must be open at all times. - Modifying inlet pressure ~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/source/user_guide/basics.rst b/docs/source/user_guide/basics.rst new file mode 100644 index 0000000..23ba918 --- /dev/null +++ b/docs/source/user_guide/basics.rst @@ -0,0 +1,47 @@ +Basics +====== + +The philosophy of Lizzy is inspired to popular FE libraries like `FEniCS `_: we provide a flexible scripting API that allows the user to define any type of infusion scenario. Everything in Lizzy is accomplished by writing a Python script. + +Lizzy workflow +----------------- + +The basic workflow for using Lizzy in a Python script involves the following key steps: + +1. Import the Lizzy package. +2. Instantiate a :class:`~lizzy.LizzyModel` object. +3. Use the :class:`~lizzy.LizzyModel` API to define and run the simulation. +4. Save results. + +Lizzy is designed so that most of the operations can be performed using the :class:`~lizzy.LizzyModel` class. + +Units and conventions +---------------------- +The solver assumes consistent units and does not enforce any. The user is free to use any system of units, as long as they are consistent throughout the script. We recommend to stick to the SI units: + +- Permeability: m² +- Viscosity: Pa·s +- Length: m +- Time: s +- Pressure: Pa +- Velocity: m/s + + +The LizzyModel class +--------------------- + +The :class:`~lizzy.LizzyModel` class is the main protagonist of any Lizzy script. This class provides APIs to execute all the core fucntionalities of the solver. +When writing a Lizzy script, typically the first step (after importing the library) is to instantiate a :class:`~lizzy.LizzyModel` object: + +.. code-block:: + + model = liz.LizzyModel() + +The Lizzy API is designed so that, in most cases, the :class:`~lizzy.LizzyModel` class is the only one the user needs do access to do any operation: + +.. code-block:: + + operation_output = model.some_method(args) + + +The :ref:`using_lizzy` section of this documentation is dedicated to showing how elementary operations can be done in Lizzy. The :ref:`api_reference_index` provides detailed information about the LizzyModel functionalities anf other namespaces of Lizzy. diff --git a/docs/source/user_guide/components/inlet_operations.rst b/docs/source/user_guide/components/inlet_operations.rst deleted file mode 100644 index fa69138..0000000 --- a/docs/source/user_guide/components/inlet_operations.rst +++ /dev/null @@ -1,24 +0,0 @@ -.. _inlet_operations: - -================== -Inlet operations -================== - - We will limit ourselves to using the :class:`~lizzy.LizzyModel` user-facing methods, as this is the way Lizzy is intended to be used. For more details about the underlying core components, please refer to the :ref:`api_reference_index` documentation. - -Preparing the mesh boundaries ------------------------------ - -Before we can create and assign inlets to the model, we need to ensure that the mesh we are using has named boundaries. Named boundaries are essential to identify where inlets should be assigned. When creating a mesh using GMSH, named boundaries can be created using Physical Groups. For example, in a 2D mesh, we can create a Line Physical Group and give it a name like "left_edge". This name will be used later when assigning the inlet to the mesh. - - -.. image:: ../../../images/physical_boundary_tag.png - :width: 80% - :align: center - -.. important:: - - When using .msh format, ensure the file is exported in Version 4 ASCII, as Lizzy currently only supports this format. - - - diff --git a/docs/source/user_guide/managing_materials.rst b/docs/source/user_guide/managing_materials.rst index 7278d64..85a39ae 100644 --- a/docs/source/user_guide/managing_materials.rst +++ b/docs/source/user_guide/managing_materials.rst @@ -8,25 +8,62 @@ Materials in Lizzy In Lizzy, a material is represented by the :class:`~lizzy.materials.PorousMaterial` class. This class encapsulates the properties of a porous material, including its permeability, porosity, and thickness. Each material is defined by the following properties: +- **Name**: A unique string identifier for the material. - **Permeability (k1, k2, k3)**: The permeability of the material in three principal directions (in m²). - **Porosity**: The porosity of the material (dimensionless, between 0 and 1). -- **Thickness**: The thickness of the material (in meters). A material can represent a single layer of fabric, or a multi-layer laminate. In the latter case, the thickness represents the total thickness of the laminate. -- **Name**: A unique string identifier for the material. +- **Thickness**: The thickness of the material (in m). A material can represent a single layer of fabric, or a multi-layer laminate. In the latter case, the thickness represents the total thickness of the laminate. The current version of Lizzy does not allow to compose a multi-layer laminate automatically by defining its layers. Instead, the user must compute the equivalent permeability, porosity, and thickness of the laminate externally and define a single :class:`~lizzy.materials.PorousMaterial` object representing the entire laminate. This is tipically done using arithmetic average schemes :cite:`calado1996effective` :cite:`bancora2018effective`. We plan to implement an automated multi-layer laminate definition feature in future releases. - Creating materials ------------------- -We can create a porous material using the :meth:`~lizzy.LizzyModel.create_material` method. This method requires the permeability values, porosity, thickness, and a unique name for the material. For example, to create a material with permeability values of 1E-10 m² in all directions, a porosity of 0.5, a thickness of 1.0 m, and named "material_01", we would do: +.. note:: + The following operations are to be performed **before** the solver is initialised by calling :meth:`~lizzy.LizzyModel.initialise_solver`. + +We can create a porous material using the :meth:`~lizzy.LizzyModel.create_material` method. This method requires a unique name for the material, the permeability values (as a tuple), porosity and thickness. For example, to create a material named "material_01" with permeability values of 1E-10 m² in all directions, a porosity of 0.5 and a thickness of 1.0 mm, we would type: + +.. code-block:: + + model.create_material("material_01", (1E-10, 1E-10, 1E-10), 0.5, 0.001) + +A :class:`~lizzy.materials.PorousMaterial` object is created and stored in the model, but it is not assigned yet. + +Creating an orientation Rosette +------------------------------- + +When working with anisotropic materials, it is necessary to define how the principal directions of permeability are oriented in the domain. To do so, we can create a :class:`~lizzy.materials.Rosette` object. A rosette is a local basis :math:`\hat e_1, \hat e_2, \hat e_3` that defines an orientation in space. When associated to a material, the principal directions of permeability will be aligned with the basis defined by the rosette. In Lizzy, we don't define the basis components directly. Instead, for any given rosette, we define a single orientation vector :math:`u_1` that will be projected onto the mesh elements to calculate the local basis :math:`\hat e_1, \hat e_2, \hat e_3` (more on this below). To create a rosette, we can use the :meth:`~lizzy.LizzyModel.create_rosette` method by providing a name for the rosette and the vector :math:`u_1`. For example, we can create a rosette named "rosette_01" with a vector :math:`u_1` oriented along (1, 1, 0) by typing: .. code-block:: - model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "material_01") + model.create_rosette("rosette_01", (1,1,0)) + +Assigning materials +-------------------- + +To assign a material to a labeled domain, we use the :meth:`~lizzy.LizzyModel.assign_material` method, providing the following arguments: + +- the name of the material to assign +- the name of the mesh domain where we want to assign it +- the rosette of orientation + +For example, to assign the material "material_01" to a mesh domain labelled as "domain_01" and orient it using using the rosette "rosette_01", we would type: + +.. code-block:: + + model.assign_material("material_01", "domain_01", "rosette_01") + +When the assignment happens, the remaining components of the local rosette (:math:`\hat e_1, \hat e_2, \hat e_3`) are calculated and registered for as follows: + +- :math:`\hat e_1`: is calculated as the projection of :math:`u_1` along the normals of the elements in the domain +- :math:`\hat e_3`: is calculated as the normal of the elements in the domain +- :math:`\hat e_2`: is calculated as :math:`\hat e_1 \times \hat e_3` -References ----------- +.. note:: + If a rosette is not provided, a default orientation is used with :math:`\hat e_1, \hat e_2, \hat e_3` aligned with global axes :math:`e_1, e_2, e_3`. This can be convenient when working with isotropic materials, since we don't care about their orientation. For example: -.. bibliography:: ../refs.bib - :style: unsrt \ No newline at end of file + .. code-block:: + + k_iso = 1.0E-11 + model.create_material("material_iso", (k_iso, k_iso, k_iso), 0.5, 1.0) + model.assign_material("material_iso", "domain_01") diff --git a/docs/source/user_guide/reading_a_mesh.rst b/docs/source/user_guide/reading_a_mesh.rst index 1d64654..736ab09 100644 --- a/docs/source/user_guide/reading_a_mesh.rst +++ b/docs/source/user_guide/reading_a_mesh.rst @@ -6,14 +6,21 @@ Lizzy can read a mesh file that contains tagged domains and boundaries. .. important:: Currently, the only supported mesh format is `.msh` (Version 4 ASCII). This is a typical mesh format that can be exported using software like GMSH. Continuing development, more formats will gradually be supported. -Mesh preparation ----------------- -The `.msh` file has some requirements to be compatible with Lizzy: +Mesh preparation (GMSH version) +------------------------------- -- The mesh must be composed uniquely of 2D triangle elements (quads are not supported). +The GMSH `.msh` file has some requirements to be compatible with Lizzy: + +- The mesh must be composed uniquely of 2D triangle elements (quads or 3D elements are not supported yet). - Physical groups must be defined for all material domains and for all edges where a boundary condition will be applied. The name of the physical groups will be used to identify the different boundaries and materials in the mesh. +.. figure:: ../../images/physical_boundary_tag.png + :width: 60% + :align: center + + An example of a GMSH mesh, where a physical boundary is being defined. + As an example, we provide a :download:`GMSH script <../../../examples/meshes/Rect.geo>` (.geo) that demonstrated how to create a simple rectangular mesh with physical groups for Lizzy. The file can be opened in GMSH and exported as a `.msh` file. To learn more about how to prepare mesh files and tag physical entities using GMSH, refer to the `GMSH documentation `_ . diff --git a/docs/source/user_guide/saving_results.rst b/docs/source/user_guide/saving_results.rst new file mode 100644 index 0000000..7eccae9 --- /dev/null +++ b/docs/source/user_guide/saving_results.rst @@ -0,0 +1,7 @@ +.. _saving_results: + +Saving results +============== + +In this section we look at how we can save simulation results to files that can be read externally. +All operations can be performed using the :class:`~lizzy.LizzyModel` user-facing methods. For more details about the underlying core components, please refer to the :ref:`api_reference_index` documentation. diff --git a/docs/source/user_guide/solver_operations.rst b/docs/source/user_guide/solver_operations.rst index 8d7cb85..3d91fd3 100644 --- a/docs/source/user_guide/solver_operations.rst +++ b/docs/source/user_guide/solver_operations.rst @@ -1,7 +1,7 @@ Running a simulation ==================== -In this section we look at various solver operations, including picking a solver, running a filling simulation and saving results. +In this section we look at various solver operations, including picking a solver and running a filling simulation. All operations can be performed using the :class:`~lizzy.LizzyModel` user-facing methods. For more details about the underlying core components, please refer to the :ref:`api_reference_index` documentation. The solver is a core component of Lizzy, responsible for solving the governing equations that describe flow through porous media. Lizzy uses a finite element / control volume method (FE/CV) to simulate the filling of the part. According to this method, the simulation is discretised into quasi-static time steps. @@ -12,18 +12,164 @@ At every step, the following operations are carried out: - a time step is calculated - the fill state of each control volume at the flow front is updated based on the velocity field and the time step duration -Currently, Lizzy supports three solvers to compute the pressure solution: +Solver initialization +--------------------- + +Every Lizzy script is broadly divided in two parts: + +- Model definition: the model is created and all elements of the simulation are set. This includes creating mesh, materials, sensors, intlets, etc... +- Solution: the simulation is run and the solution is computed. + +The solver initialisation is what separates these two parts. Before it can run the simulation, the solver must always be initialised by calling the :meth:`~lizzy.LizzyModel.initialise_solver` method: + +.. code-block:: + + model.initialise_solver() + +The solver initialization is a critical point in every Lizzy simulation, because it is the point where a lot of object data types are converted into faster types for performance. As a consequence, several methods do not work anymore after the solver initialisation critical point. + +.. + Maybe add an example pictorial illustrated phases here + + +When we initialise the solver, we can also pick the solver that we want to use. Lizzy currently supports three solvers to compute the pressure solution: - **Direct dense solver:** uses a direct method to solve the linear system. There is generally no benefit to using this solver, and will be soon discontinued. - **Direct sparse solver:** uses a direct method and sparse matrix allocation to solve the linear system. This is the default solver in Lizzy if PETSc is not available. Suitable for small to medium-sized problems. - **Iterative PETSc solver:** uses an iterative method and sparse matrix allocation to solve the linear system. This solver is generally faster and more memory efficient than the direct solvers. This is the default solver in Lizzy if PETSc is available. Users should aim at using this solver if possible. Relies on the PETSc library and the petsc4py package. See the :ref:`installation` page for more details about installing and checking dependencies. +The solver types are available as an Enum :class:`~lizzy.SolverType` in the ``lizzy`` namespace. To pick a solver, we just pass the Enum as an argument to the :meth:`~lizzy.LizzyModel.initialise_solver` method: + +.. code-block:: + + model.initialise_solver(lizzy.SolverType.ITERATIVE_PETSC) + +The :meth:`~lizzy.LizzyModel.initialise_solver` method also accepts other arguments that are specific to certain solvers. Normally, the user does not need to specify any. However, the API reference contains detailed description of all parameters. + Time step vs time interval -------------------------- Throughout this section, we will encounter multiple times the terms "time step" and "time interval". It is important to clarify the difference between these two concepts: -- **Time step:** - The discrete increment of time used by the solver to advance the simulation. The time step is determined by the solver at runtime and **the user has no control over this quantity**. -- **Time interval:** - An amount of time over which the simulation advances. A time interval is tipically composed of multiple time steps. **The user has full control over this quantity**. For example, if a simulation is run until the part is completely filled, then the time interval is the entird fill time. Conversely, we can set our simulation to run for a fixed time interval, e.g., 60 seconds, then pause and do something, and then resume the simulation for another time interval... and so on. \ No newline at end of file +- **Time step:** The discrete increment of time used by the solver to advance the simulation. The time step is determined by the solver at runtime. **The user has no control over this quantity**. +- **Time interval:** An amount of time over which the simulation advances. A time interval is tipically composed of multiple time steps. **The user has full control over this quantity**. For example, if a simulation is run until the part is completely filled, then the time interval is the entird fill time. Conversely, if we set our simulation to run for a fixed time interval, e.g., 60 seconds, then pause and do something, and then resume the simulation, then the time interval is 60 seconds. + +Running a simulation until the part is filled +--------------------------------------------- + +To run the simulation until the part is completely filled, we can call the :meth:`~lizzy.LizzyModel.solve` method: + +.. code-block:: + + model.solve() + +This will start the filling simulation from the latest state (if called for the first time, from the empty part) and keep filling until :attr:`~lizzy.LizzyModel.n_empty_cvs` reaches 0. + +.. tip:: + By default, we get a log of the progress: + + .. code-block:: + + >>> Fill time: 1001.25s, Empty CVs: 2231 + + We can suppress that log by passing: + + .. code-block:: + + model.solve(log="off") + +Running a simulation for a time interval +---------------------------------------- + +We can run the simulation only for an interval of fill time and then pause by using the method :meth:`~lizzy.LizzyModel.solve_time_interval`. This method needs an argument specifying the duration of the time interval. For example: + +.. code-block:: + + model.solve_time_interval(300) + +will advance the filling simulation for 300 seconds of process time. We can call this method multiple times, and even combine it with :meth:`~lizzy.LizzyModel.solve`: + +.. code-block:: + + model.initialise_solver() # the part is empty + model.solve_time_interval(300) # fill for 300 seconds + model.solve_time_interval(600) # advance filling for another 600 seconds + model.solve() # continue filling until part is filled + +Creating dynamic filling scenarios +---------------------------------- + +The :meth:`~lizzy.LizzyModel.solve_time_interval` method is very powerful when combined to gate management methods like :meth:`~lizzy.LizzyModel.open_inlet` or :meth:`~lizzy.LizzyModel.change_inlet_pressure`. We can create dynamic scenarios, for example: + +.. code-block:: + + # assume a model that has 2 inlets "inlet_1" and "inlet_2", both open: + model.initialise_solver() # the part is empty + model.solve_time_interval(600) # fill for 600 seconds with both inlets open + model.close_inlet("inlet_1") # close "inlet_1" + model.solve_time_interval(600) # advance filling with only one inlet open + model.open_inlet("inlet_1") # open "inlet_1" again... + model.change_inlet_pressure("inlet_1", 150000) # ...and set its pressure to 1.5 bar + model.solve() # continue filling until part is completely filled + +The Solution datatype +--------------------- + +Lizzy stores solutions using :class:`~lizzy.datatypes.Solution` objects. A Solution object can be seen as a snapshot that stores the current state of the infusion at the time of its creation. When :meth:`~lizzy.LizzyModel.solve` or :meth:`~lizzy.LizzyModel.solve_time_interval` terminate their execution, a Solution object is created. The Solution is returned by the methods, and also stored automatically in the read-only :attr:`~lizzy.LizzyModel.latest_solution` property of the model. + +.. tip:: + + Two ways to retrieve the current solution: + + - get the latest from the model, after solving: ``solution = model.latest_solution``. + + - catch it from the solving method itself: ``solution = model.solve()`` / ``solution = model.solve_time_interval(100)``. + + Note that, since the same solution object returned is also stored in :attr:`~lizzy.LizzyModel.latest_solution`, usually the manual capture is not necessary (but is given as a possibility for advanced uses). + +A Solution objects contains all the solution fields computed by the solver, for all time states that were marked for write-out. This is better explained by example: + +Let's assume a model is set up, and we specify to save a solution state every 100 seconds of fill time. Then we solve and the simulation completes at 550 seconds of fill time: + +.. code-block:: + + ... + model.assign_simulation_parameters(wo_delta_time=100) + ... + model.solve() + + >>> Fill time: 550.00s, Empty CVs: 0 + +At that point, the Solution object will contain results for all these time steps: + +- the initial time step (t = 0) +- 5 sequent time intervals every 100 seconds (t=100, t=200, ... , t=500) +- the final time step when the part was filled (t=550) + +So in this particular Solution we will have 7 time states. We can get these states by accessing the attributes of the Solution object. For example, let's take a look at the shape of the time, pressure and velocity field: + +.. code-block:: + + solution = model.latest_solution + print(f"Number of time states stored: {solution.n_time_states}") + print(f"Time field shape: {solution.time.shape}") + print(f"Pressure field shape: {solution.p.shape}") + print(f"Velocity field shape: {solution.v.shape}") + + >>> Number of time states stored: 7 + Time field shape: (7,) + Pressure field shape: (7, 735) + Velocity field shape: (7, 1366, 3) + + # Let's print the content of solution.time to see the values: + print(solution.time) + + >>> [ 0. 100. 200. 300. 400. 500. 550.] + +Each field that is stored in the Solution has a leading dimension equal to the number of time states present in the Solution. So for example, if we want the pressure field at t=200s (third time state), we can type: ``solution.p[2]``. +Consult the :class:`~lizzy.datatypes.Solution` API reference to get more information on the fields stored. + +The Solution object is also used by Lizzy to write result files. More information on this in :ref:`saving_results`. + + + diff --git a/docs/source/user_guide/using_lizzy.rst b/docs/source/user_guide/using_lizzy.rst index c64ad6b..fb7f709 100644 --- a/docs/source/user_guide/using_lizzy.rst +++ b/docs/source/user_guide/using_lizzy.rst @@ -1,60 +1,20 @@ +.. _using_lizzy: + Using Lizzy =========== -In this section, we detail how to use Lizzy setting up and running simulations. The philosophy of Lizzy is inspired to popular FE libraries like `FEniCS `_: we provide a flexible scripting API that allows the user to define any type of infusion scenario. Everything in Lizzy is accomplished by writing a Python script. - - -Basic script workflow ----------------------- - -The basic workflow for writing a Python script to use Lizzy involves the following key steps: - -1. Import the Lizzy package. -2. Instantiate a :class:`~lizzy.LizzyModel` object. -3. Use the :class:`~lizzy.LizzyModel` API to define and run the simulation. -4. Save results. - -Lizzy is designed so that most of the operations are performed using the :class:`~lizzy.LizzyModel` class directly. - -Units and conventions ----------------------- -The solver assumes consistent units and does not enforce any. The user is free to use any system of units, as long as they are consistent throughout the script. We recommend to stick to the SI units: - -- Permeability: m² -- Viscosity: Pa·s -- Length: m -- Time: s -- Pressure: Pa -- Velocity: m/s - - -The LizzyModel class ---------------------- - -The :class:`~lizzy.LizzyModel` class is the main protagonist of any Lizzy script. This class provides APIs to execute all the core fucntionalities of the solver. -When writing a Lizzy script, typically the first step (after importing the library) is to instantiate a :class:`~lizzy.LizzyModel` object: - -.. code-block:: - - model = liz.LizzyModel() - -The Lizzy API is designed so that, in most cases, the :class:`~lizzy.LizzyModel` class is the only one the user needs do access to do any operation: - -.. code-block:: - - operation_output = model.some_method(args) - - -In this section, we cover the most common operations that can be performed using a LizzyModel object. +In this section, we cover how to use Lizzy setting up and running simulations. .. toctree:: :maxdepth: 1 + basics reading_a_mesh assigning_boundary_conditions assigning_parameters managing_materials solver_operations + saving_results .. seealso:: @@ -63,11 +23,3 @@ In this section, we cover the most common operations that can be performed using :ref:`LizzyModel API reference ` The API reference of the LizzyModel class. - -See also -^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. toctree:: - :maxdepth: 1 - - components/inlet_operations \ No newline at end of file diff --git a/examples/meshes/quarter_annulus.geo b/examples/meshes/quarter_annulus.geo new file mode 100644 index 0000000..1208ef5 --- /dev/null +++ b/examples/meshes/quarter_annulus.geo @@ -0,0 +1,113 @@ +// ============================================ +// GMSH Geometry File: Quarter Annular Disc +// ============================================ + +// -------------------------------------------- +// Parameters (modify these as needed) +// -------------------------------------------- +r_0 = 1.0; // Inner radius +r_f = 2.0; // Outer radius + +// Mesh parameters +mesh_inner = 0.04; // Mesh size at inner radius (finer) +mesh_outer = 0.12; // Mesh size at outer radius (coarser) + +// -------------------------------------------- +// Geometry Definition +// -------------------------------------------- + +// Center point +Point(1) = {0, 0, 0, mesh_inner}; + +// Inner arc points (quarter circle) +Point(2) = {r_0, 0, 0, mesh_inner}; // Inner radius at 0 degrees +Point(3) = {0, r_0, 0, mesh_inner}; // Inner radius at 90 degrees +Point(4) = {r_0/Sqrt(2), r_0/Sqrt(2), 0, mesh_inner}; // Inner radius at 45 degrees (for arc definition) + +// Outer arc points (quarter circle) +Point(5) = {r_f, 0, 0, mesh_outer}; // Outer radius at 0 degrees +Point(6) = {0, r_f, 0, mesh_outer}; // Outer radius at 90 degrees +Point(7) = {r_f/Sqrt(2), r_f/Sqrt(2), 0, mesh_outer}; // Outer radius at 45 degrees (for arc definition) + +// -------------------------------------------- +// Lines and Arcs +// -------------------------------------------- + +// Inner arc (inlet) - using Circle arc with center point +Circle(1) = {2, 1, 4}; // Inner arc first half (0 to 45 degrees) +Circle(2) = {4, 1, 3}; // Inner arc second half (45 to 90 degrees) + +// Outer arc (outlet) - using Circle arc with center point +Circle(3) = {5, 1, 7}; // Outer arc first half (0 to 45 degrees) +Circle(4) = {7, 1, 6}; // Outer arc second half (45 to 90 degrees) + +// Radial lines connecting inner and outer arcs +Line(5) = {2, 5}; // Radial line at 0 degrees +Line(6) = {3, 6}; // Radial line at 90 degrees + +// -------------------------------------------- +// Curve Loops and Surface +// -------------------------------------------- + +// Define curve loop (counter-clockwise) +// Start at inner radius (0 deg), go along inner arc, then radial, then outer arc back, then radial +Curve Loop(1) = {1, 2, 6, -4, -3, -5}; + +// Create plane surface +Plane Surface(1) = {1}; + +// -------------------------------------------- +// Mesh Control - Finer at inner radius +// -------------------------------------------- + +// Distance field from inner arc for graded mesh +Field[1] = Distance; +Field[1].CurvesList = {1, 2}; +Field[1].Sampling = 100; + +Field[2] = Threshold; +Field[2].InField = 1; +Field[2].SizeMin = mesh_inner; +Field[2].SizeMax = mesh_outer; +Field[2].DistMin = 0; +Field[2].DistMax = r_f - r_0; + +Background Field = 2; + +// -------------------------------------------- +// Mesh Settings +// -------------------------------------------- + +// Set mesh algorithm (6 = Frontal-Delaunay for triangles) +Mesh.Algorithm = 6; + +// Force linear triangular elements (order 1) +Mesh.ElementOrder = 1; + +// Ensure 2D mesh +Mesh.Smoothing = 5; + +// -------------------------------------------- +// Physical Groups (Entity Sets) +// -------------------------------------------- + +// Inlet: Inner arc curves +Physical Curve("inlet", 100) = {1, 2}; + +// Outlet: Outer arc curves +Physical Curve("outlet", 200) = {3, 4}; + +// Radial boundaries (if needed) +Physical Curve("symmetry_x", 300) = {5}; // Line along x-axis +Physical Curve("symmetry_y", 400) = {6}; // Line along y-axis + +// Domain: All surface elements +Physical Surface("domain", 500) = {1}; + +// -------------------------------------------- +// Generate Mesh +// -------------------------------------------- +Mesh 2; + +// Optional: Save mesh (uncomment if needed) +// Save "quarter_annulus.msh"; diff --git a/examples/meshes/quarter_annulus.msh b/examples/meshes/quarter_annulus.msh new file mode 100644 index 0000000..904a125 --- /dev/null +++ b/examples/meshes/quarter_annulus.msh @@ -0,0 +1,2484 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 100 "inlet" +1 200 "outlet" +1 300 "symmetry_x" +1 400 "symmetry_y" +2 500 "domain" +$EndPhysicalNames +$Entities +7 6 1 0 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +4 0.7071067811865475 0.7071067811865475 0 0 +5 2 0 0 0 +6 0 2 0 0 +7 1.414213562373095 1.414213562373095 0 0 +1 0.7071067811865476 0 0 1 0.7071067811865475 0 1 100 2 2 -4 +2 1.110223024625157e-16 0.7071067811865476 0 0.7071067811865476 1 0 1 100 2 4 -3 +3 1.414213562373095 0 0 2 1.414213562373095 0 1 200 2 5 -7 +4 2.220446049250313e-16 1.414213562373095 0 1.414213562373095 2 0 1 200 2 7 -6 +5 1 0 0 2 0 0 1 300 2 2 -5 +6 0 1 0 0 2 0 1 400 2 3 -6 +1 0 0 0 2 2 0 1 500 6 1 2 6 -4 -3 -5 +$EndEntities +$Nodes +13 608 1 608 +0 2 0 1 +1 +1 0 0 +0 3 0 1 +2 +0 1 0 +0 4 0 1 +3 +0.7071067811865475 0.7071067811865475 0 +0 5 0 1 +4 +2 0 0 +0 6 0 1 +5 +0 2 0 +0 7 0 1 +6 +1.414213562373095 1.414213562373095 0 +1 1 0 19 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +0.9992290362375071 0.03925981584091626 0 +0.9969173337188422 0.07845909590936345 0 +0.9930684569202419 0.1175373977508844 0 +0.9876883405308875 0.1564344654458907 0 +0.9807852803043174 0.1950903225133976 0 +0.9723699202480781 0.2334453644790282 0 +0.9624552362447066 0.2714404506059221 0 +0.9510565160271076 0.3090169951999081 0 +0.9381913355961035 0.3461170579621865 0 +0.9238795321039749 0.3826834333484275 0 +0.9081431734219714 0.4186597384118412 0 +0.8910065237881879 0.4539905005249443 0 +0.872496006718287 0.4886212421299784 0 +0.8526401640321339 0.522498565241337 0 +0.8314696120097427 0.5555702334578122 0 +0.8090169941335599 0.5877852526247146 0 +0.785316930717919 0.6190939495163776 0 +0.7604059654846134 0.6494480484653203 0 +0.7343225093741095 0.6788007455995543 0 +1 2 0 19 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +0.6788007455342941 0.7343225094344354 0 +0.6494480482136891 0.7604059656995268 0 +0.6190939490514525 0.7853169310844365 0 +0.5877852519074092 0.8090169946547128 0 +0.5555702325329381 0.8314696126277239 0 +0.5224985640153165 0.8526401647834404 0 +0.488621240660577 0.8724960075411914 0 +0.4539904987928685 0.8910065246707246 0 +0.4186597364693214 0.9081431743174853 0 +0.3826834311746453 0.9238795330043852 0 +0.3461170559968018 0.9381913363211722 0 +0.3090169933486863 0.951056516628606 0 +0.2714404490025435 0.9624552366969061 0 +0.2334453630931222 0.9723699205808047 0 +0.195090321319784 0.980785280541742 0 +0.1564344644781432 0.9876883406841637 0 +0.1175373970741732 0.993068457000336 0 +0.07845909546510149 0.9969173337538064 0 +0.03925981563047087 0.9992290362457756 0 +1 3 0 13 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +1.996853630022355 0.1121408947108161 0 +1.987424419724236 0.2239289527590905 0 +1.97174203688675 0.3350124474906743 0 +1.949855824095288 0.4450418690883855 0 +1.921834643455892 0.5536710243572111 0 +1.887766659995207 0.66055812568656 0 +1.847759064206153 0.7653668667011921 0 +1.801937734993994 0.8677674799188524 0 +1.750446843095188 0.9674377755174186 0 +1.69344839784489 1.064064154004153 0 +1.631121737422011 1.157342593059376 0 +1.563662964658493 1.246979604065524 0 +1.491284329691518 1.332693155987724 0 +1 4 0 13 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +1.332693155760922 1.491284329894201 0 +1.246979603260582 1.563662965300413 0 +1.157342591536983 1.631121738502206 0 +1.064064151765856 1.693448399251306 0 +0.9674377725335979 1.750446844744287 0 +0.8677674761590705 1.801937736804609 0 +0.7653668623514042 1.847759066007895 0 +0.6605581217040573 1.887766661388745 0 +0.553671021061469 1.921834644405379 0 +0.4450418663650489 1.949855824716872 0 +0.3350124453008868 1.97174203725881 0 +0.2239289514715623 1.987424419869306 0 +0.1121408941273829 1.99685363005512 0 +1 5 0 13 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +1.040816703746849 0 0 +1.084965407705819 0 0 +1.132718132745993 0 0 +1.184369059323481 0 0 +1.240236426131326 0 0 +1.300664451639965 0 0 +1.366025409823141 0 0 +1.436722007253737 0 0 +1.513189805963489 0 0 +1.59589993923363 0 0 +1.685361996518571 0 0 +1.782127115571373 0 0 +1.886791505768078 0 0 +1 6 0 13 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +0 1.040816703746849 0 +0 1.084965407705819 0 +0 1.132718132745993 0 +0 1.184369059323481 0 +0 1.240236426131326 0 +0 1.300664451639965 0 +0 1.366025409823141 0 +0 1.436722007253737 0 +0 1.513189805963489 0 +0 1.59589993923363 0 +0 1.685361996518571 0 +0 1.782127115571373 0 +0 1.886791505768078 0 +2 1 0 512 +97 +98 +99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +109 +110 +111 +112 +113 +114 +115 +116 +117 +118 +119 +120 +121 +122 +123 +124 +125 +126 +127 +128 +129 +130 +131 +132 +133 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 +256 +257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 +393 +394 +395 +396 +397 +398 +399 +400 +401 +402 +403 +404 +405 +406 +407 +408 +409 +410 +411 +412 +413 +414 +415 +416 +417 +418 +419 +420 +421 +422 +423 +424 +425 +426 +427 +428 +429 +430 +431 +432 +433 +434 +435 +436 +437 +438 +439 +440 +441 +442 +443 +444 +445 +446 +447 +448 +449 +450 +451 +452 +453 +454 +455 +456 +457 +458 +459 +460 +461 +462 +463 +464 +465 +466 +467 +468 +469 +470 +471 +472 +473 +474 +475 +476 +477 +478 +479 +480 +481 +482 +483 +484 +485 +486 +487 +488 +489 +490 +491 +492 +493 +494 +495 +496 +497 +498 +499 +500 +501 +502 +503 +504 +505 +506 +507 +508 +509 +510 +511 +512 +513 +514 +515 +516 +517 +518 +519 +520 +521 +522 +523 +524 +525 +526 +527 +528 +529 +530 +531 +532 +533 +534 +535 +536 +537 +538 +539 +540 +541 +542 +543 +544 +545 +546 +547 +548 +549 +550 +551 +552 +553 +554 +555 +556 +557 +558 +559 +560 +561 +562 +563 +564 +565 +566 +567 +568 +569 +570 +571 +572 +573 +574 +575 +576 +577 +578 +579 +580 +581 +582 +583 +584 +585 +586 +587 +588 +589 +590 +591 +592 +593 +594 +595 +596 +597 +598 +599 +600 +601 +602 +603 +604 +605 +606 +607 +608 +0.9750717696158071 0.3366952773707361 0 +0.3343480130020396 0.9797950420390762 0 +0.8253157435959376 0.6240147444992401 0 +0.5250560719899192 0.8910258896619225 0 +0.7183358447234928 0.7426884800511779 0 +0.9136301666816947 0.4873476828015773 0 +0.06415310390502063 1.396622865816687 0 +1.400628297535887 0.05746546118469935 0 +0.6223362376403836 0.827318868485367 0 +0.7887937614224111 1.742904717702139 0 +1.732546226235368 0.7876680382982323 0 +1.451817742682154 1.220232453479168 0 +1.009957883348557 0.2189657304382422 0 +0.21858759088985 1.008996519421558 0 +1.145085800634323 1.531130047634773 0 +0.04451116845186544 1.217954376226799 0 +1.214881463630588 0.04863943843733343 0 +0.4891073262199337 1.834003045589661 0 +1.842127345593992 0.4691718947222702 0 +0.07530834904694154 1.642343612011399 0 +1.650625378417631 0.06948560464422536 0 +0.9467527098204658 0.4143198080031569 0 +0.7724689437707021 0.6840066382731227 0 +0.8695812190346195 0.5566034944911754 0 +0.4530351177152666 0.9325698710869142 0 +1.02159602576616 0.1391670809135997 0 +0.1440607114278919 1.03021339976804 0 +1.311830102613035 1.379989607455256 0 +1.570610398115707 1.064886524020957 0 +0.9719479908546648 1.64462162998192 0 +1.104735901487112 0.04231040753333555 0 +0.03498785817379871 1.107273192242248 0 +0.2662165237713827 1.869919188267699 0 +1.885168244869961 0.2532909836553252 0 +0.06081393878480611 1.031458767893607 0 +1.029893443704415 0.06053850703538683 0 +1.836429869765211 0.07960373655911218 0 +0.08722153604095265 1.834459310669726 0 +0.07282018610668095 1.554012758955574 0 +0.1464829347815063 1.596811852014825 0 +0.1476222913807396 1.503590924711242 0 +0.2188785137371423 1.555786232776368 0 +0.2256588108450855 1.478371112014669 0 +0.2879032380761391 1.517435044179453 0 +0.2894757424847788 1.442005929185004 0 +0.3543796674980318 1.478774490839553 0 +0.3545221792551996 1.404677161686404 0 +0.4204645117047413 1.440336982539778 0 +0.4207967388523445 1.369167785249609 0 +0.4834348052389895 1.400106108435671 0 +0.4857928750663499 1.333132927417778 0 +0.5468628384219774 1.364338739269415 0 +0.5505951089926039 1.297122883030203 0 +0.6110502163720875 1.331780915635299 0 +0.6151080966698622 1.26019526499878 0 +0.6768257105943442 1.295419193253333 0 +0.677697216763074 1.222068148982743 0 +0.7408469903473172 1.258539842563131 0 +0.7411583031617314 1.187414111308634 0 +0.8029768987882655 1.220582817764308 0 +0.805586164232246 1.15323633014166 0 +0.8664267287702871 1.185232466387691 0 +0.8679634193851469 1.116029782853439 0 +0.931632664645019 1.151139947004919 0 +0.9314026277813531 1.077891044560816 0 +0.3576195196718622 1.335338033817056 0 +0.9982858108817579 1.115711607838983 0 +0.9980431178774944 1.040640119588015 0 +0.9336826028405106 1.004283989552674 0 +0.9966423084836289 0.9660838533660465 0 +1.062437330686895 1.00243889601783 0 +1.05880835239419 0.9280160398667427 0 +1.12525259561008 0.9648603434855412 0 +1.120140979838883 0.8911559933788884 0 +1.059481524985423 0.8551008151622084 0 +1.119689509541492 0.8189116780639535 0 +1.183435430311209 0.853850358877002 0 +1.178738465999546 0.7799286578310367 0 +1.241436894723046 0.8131027515944242 0 +1.237507008665568 0.7405154716714947 0 +1.173680926754597 0.7059621675570236 0 +1.233336069181417 0.6673630939602715 0 +0.4899390239835744 1.266078550848807 0 +0.7465116102264342 1.121510934591656 0 +1.172414869398377 0.6370638926141308 0 +1.226807971339517 0.5935237010764862 0 +1.291947775552961 0.6258965028139676 0 +1.284511731584737 0.5519382093863987 0 +1.354537676197181 0.5839157975337883 0 +1.346329511830674 0.5091090092215204 0 +1.278106056789575 0.4796061161995143 0 +1.336970846430419 0.4361368752054932 0 +1.271560264560173 0.4108076179560232 0 +0.6181193779210366 1.193641370408596 0 +1.407385036239083 0.4640956635501527 0 +0.9354365604533006 0.9332292053546394 0 +1.111084930987003 0.6726146752104083 0 +1.167911584976023 0.57357817604063 0 +1.39725838069294 0.3903266058626269 0 +0.4882124059840066 1.4689995608343 0 +0.6731409832355133 1.370211048868217 0 +0.8602328306411192 1.247082155424371 0 +1.468242739490023 0.4176657825045612 0 +1.457487294140843 0.3438069956621714 0 +1.301322633957986 0.7761721816224194 0 +1.131842698991889 1.037828097308547 0 +1.253230456261213 0.8899695971483341 0 +1.480132599971486 0.4921624962586962 0 +1.388186029651802 0.3155811368551571 0 +0.99762038117408 1.194088806504861 0 +1.527864968739368 0.3718329949499179 0 +0.8767223990874405 0.9675206648193148 0 +0.8783464959906684 0.903846232875535 0 +0.8252062434032752 0.9323600417436395 0 +0.3527092905576642 1.553483826465075 0 +0.8263298677735764 0.8732365888048259 0 +0.8760246102769393 0.8444085785530806 0 +0.776219126792741 0.9025364705898271 0 +1.44747130364136 0.2688525748014542 0 +1.380543915506079 0.2401146105802878 0 +0.2299711810721651 1.410813304711242 0 +0.825842523411688 0.809588677124766 0 +0.8776078823889014 0.7873577622389236 0 +0.2213795999459891 1.643234867030123 0 +0.7738318889389126 0.9585503463759789 0 +0.7268964697742901 0.9291744455322712 0 +0.7291670247059241 0.8752672250910287 0 +0.7235033801150023 0.9843621643105742 0 +0.6778100516145584 0.9540266669929789 0 +0.6732064051947253 1.007586697719035 0 +0.6288514188931755 0.9780056699303731 0 +0.6245726540808529 1.032310770064526 0 +0.5806353780413281 1.004062713550636 0 +0.5785892816179792 1.058926138831759 0 +0.5328365721929792 1.031011895349724 0 +0.5275524684333691 1.083360356904136 0 +0.484412140597256 1.055892495072849 0 +0.4783798391204953 1.105516942147971 0 +0.4367399310758698 1.079099735734734 0 +0.4316961409994102 1.128660680622334 0 +0.3881438489571998 1.102571890372614 0 +0.3837359776723542 1.154874912555332 0 +0.3396333138893559 1.13031246256062 0 +0.3342966720508148 1.179422848306405 0 +0.2894256809768301 1.153205163011988 0 +0.2869271550763744 1.205261999336247 0 +0.2392375884111763 1.179498960820843 0 +0.2379877542970392 1.236970352699562 0 +0.1890391092341191 1.207880670338746 0 +0.1900813995771821 1.153276182781194 0 +0.1865288106057882 1.266548274273805 0 +0.2961696769729216 1.10060499123254 0 +1.216804666875589 0.4487445800331787 0 +1.210368100612826 0.3862661378352557 0 +1.160864872341217 0.420852133432793 0 +1.154380972074237 0.3623917480782608 0 +1.201135285484459 0.3259575519432497 0 +1.145502026104667 0.3043504528792477 0 +1.190559668222042 0.2676655933338872 0 +1.132460165843933 0.2432487601226246 0 +1.181871393509594 0.2116990191211017 0 +1.108642075261353 0.3973519611457933 0 +1.134229901323663 0.188693102305209 0 +1.113185173759765 0.4517952093707419 0 +1.091173210465958 0.2889357326578693 0 +0.9202801839363083 0.8144639147459749 0 +0.919211149542974 0.7597492786536122 0 +0.8781599931354183 0.7350043545207978 0 +0.9185325772922786 0.706420253905698 0 +0.961672173249015 0.730469905125141 0 +0.9606079694530463 0.6765307696747024 0 +1.065393417281894 0.4309359013330981 0 +1.069254205565182 0.4818646054036834 0 +0.6333360304746689 0.9246959824316021 0 +0.4430690800028255 1.027992468189344 0 +1.175370213577283 0.1621719278274584 0 +1.225024444544895 0.1781960479882186 0 +1.24756282224411 0.2846005118328726 0 +1.126839535188133 0.1441535596338993 0 +0.9159814669721807 0.6534161603737441 0 +0.9573883995689119 0.6241423265548959 0 +1.43571591306987 0.1932052106143322 0 +1.504946770373345 0.2221761405077029 0 +0.3772984227980178 1.207711621377413 0 +1.002670105723065 0.6455954273271476 0 +0.9969477387944838 0.5934317328376186 0 +1.492571532793412 0.1432735055847609 0 +1.058405104894839 0.3783357510179134 0 +0.1429752150006053 1.17760106206556 0 +0.1454924084822835 1.125892610568476 0 +0.189334163466356 1.101199399806917 0 +1.068674618825445 1.156148880184051 0 +1.066920751819472 1.239063186226783 0 +1.142538073650294 1.196640000828064 0 +1.141078680399858 1.28303735271115 0 +1.061823973537896 1.325470635919135 0 +1.219428664833714 1.237436713193385 0 +1.218536119996336 1.150113787078126 0 +0.5359156294985142 0.9789203667993324 0 +1.367738621695559 0.1637018943920193 0 +1.190766348255617 0.9993466921915543 0 +0.952687920789989 0.5738323745672218 0 +0.9908013487351671 0.5426667040985989 0 +1.035270163074184 0.5635982357997935 0 +1.023852777431684 0.4629347918832364 0 +1.361489522230069 0.6607025762739778 0 +1.426727284334255 0.6169251930457498 0 +1.433388749902001 0.6971605284261995 0 +1.502059971220608 0.6511568918862312 0 +1.508249294221028 0.7347278320905142 0 +1.580372333228862 0.6868384266441186 0 +1.572042433418347 0.6025348367777777 0 +1.654839795819387 0.6387038194120858 0 +1.644848441733309 0.5494523636937676 0 +1.436985917507242 0.7804843732873517 0 +1.513612119578366 0.8219155830282501 0 +1.436590248618384 0.8695490007819574 0 +0.7676378121042881 1.012331721855892 0 +1.061619006505992 0.7862444749419926 0 +0.5180320099936435 1.135656610479865 0 +1.564836697922241 0.174395301059501 0 +1.579568077598444 0.2529797583155292 0 +1.640499452139587 0.205699979084631 0 +1.657893010662176 0.2877273315989123 0 +1.722163091571033 0.2347440967620715 0 +1.740632654792654 0.3244679887700906 0 +1.668180325123994 0.3816528494375324 0 +0.6322852927892915 1.090512547114313 0 +0.2443480673306609 1.301443273881876 0 +1.296933086068858 1.188308518274949 0 +1.296570966229539 1.096597046823363 0 +1.087818370075722 0.1716274222060415 0 +1.73309349408522 0.5853398425520919 0 +1.518964326699472 0.9110393508047836 0 +1.439937777117054 0.9626488974095239 0 +1.59520786992432 0.857591889430415 0 +1.139190497707465 1.371263265794764 0 +1.054468018322577 1.411509255966372 0 +0.9794481106680856 1.368884254069838 0 +0.975606905134528 1.456770355096074 0 +0.8969342287915857 1.411999400121515 0 +0.8927894832377127 1.502071263321296 0 +0.8155352257766576 1.457356046166858 0 +0.8101433950837771 1.545664391945406 0 +0.7347226709866914 1.499272620440403 0 +0.7270724619539154 1.586264994053826 0 +0.6517005936769631 1.539898961651285 0 +0.6442074261429749 1.6274887768468 0 +0.5688497429587234 1.584206959791459 0 +0.5630939429512792 1.672126330768489 0 +0.4888305360204228 1.631066785414296 0 +0.0995494611672061 1.148153092395699 0 +0.1026880752637668 1.097697129012315 0 +1.085538550386784 0.1291036702893406 0 +0.6352451430183194 1.710219834593238 0 +1.215168071033625 0.125073621700699 0 +1.268421787234955 0.1439133029844726 0 +1.275774580311387 0.196157728101181 0 +0.8175963572326777 1.371287203560819 0 +1.345864859301421 0.9134260403591684 0 +0.1791248947917659 1.326619893143129 0 +0.130567941644559 1.293243564833191 0 +1.114208886928072 0.5032687636312463 0 +0.4865601150630425 1.722927543246075 0 +0.407732374993505 1.672281090427933 0 +1.2612790947271 0.08668823572267063 0 +1.178297257204718 0.0870845313620411 0 +0.397362482859799 1.762313141983895 0 +0.3100455886085291 1.699304369433979 0 +1.006969867066336 0.8183031559025364 0 +1.729277331079046 0.4878152208450949 0 +1.100816529730125 0.3426129370381317 0 +1.046795789579697 0.3243119764497056 0 +1.045236563081774 0.2760881829914157 0 +1.07637177464051 0.2432774330821306 0 +1.319094614669997 0.1156231508521268 0 +0.9916396895435595 0.887924129341686 0 +0.1483261364185847 1.077428797349172 0 +0.1932030321140971 1.050125659013123 0 +0.2266228742749443 1.077607715776918 0 +0.1229818508235762 1.358286265105221 0 +0.06337522799904137 1.319898101462418 0 +0.1245336609952814 1.427704582334475 0 +1.016989874888987 0.2975464163027901 0 +0.6810666052354886 0.9003350420340186 0 +0.6838619487998828 0.8476885258595913 0 +0.7307429378400747 0.8252424806668462 0 +0.696508776716198 0.8058535108584484 0 +0.6362688563007322 0.8705671174636433 0 +0.5937800747788486 0.8949452581768409 0 +0.5830467119644682 0.9452014403281023 0 +0.4895223335571213 1.004811776955072 0 +0.4231625024958504 1.302033815738966 0 +1.117357989626791 0.7474889604177181 0 +1.059546652478764 0.7212975555065591 0 +0.2422176304451849 1.120297644633965 0 +0.2845480080916609 1.362345363312832 0 +0.5552924125556575 1.229926317538007 0 +1.121975573012293 0.6096390182545403 0 +1.019436900048708 0.4136908099518517 0 +0.6733802419618825 1.148645498771958 0 +1.212938848222505 0.5210634820151683 0 +1.379945404315192 1.140462166507198 0 +1.695877417162402 0.1481830378977935 0 +0.8084882422864169 1.088664352104843 0 +0.8640793760388396 1.040387706213497 0 +0.7765583461084123 0.8467560013111399 0 +0.77104859612647 0.7946645655895468 0 +0.8050459907532166 0.7580679579819321 0 +0.7604041612031949 0.7512571414593318 0 +1.364064413909809 1.027613477207641 0 +1.457334645598812 1.065211964363801 0 +0.4493209643204615 0.9784679022795709 0 +0.4058838760536905 0.9975698900433096 0 +0.3947033103961011 1.050911008348653 0 +0.3421522902230988 1.067301539366369 0 +0.2919071108394329 1.054478551413668 0 +1.264732830896838 0.9873347477827898 0 +0.06821080020980082 1.472488796113883 0 +1.558600669359585 0.0844870929781029 0 +1.318230491556022 0.2856962272032679 0 +1.328536291825525 0.3646824469923361 0 +0.9437515157737607 0.5226766937879459 0 +0.3621548552077039 1.013513983006491 0 +0.4915303985224727 0.9565122929675558 0 +0.3774388195253128 1.865204622799509 0 +0.8756515235668996 0.6835525614898466 0 +0.8341610845916058 0.7142601398345738 0 +0.08148752456144381 1.736610870767709 0 +0.1753300681159163 1.790003955042747 0 +1.470423044598341 0.06808489621737779 0 +0.04936805772370569 1.158518538502977 0 +0.9149097941732809 0.6049599524788509 0 +0.871560156037305 0.6285947153323894 0 +0.1381050963365854 1.233131890746329 0 +1.028117579313072 0.09974497837459814 0 +0.9308766285275418 0.4503401672391823 0 +0.8940270795345389 0.5226783750605783 0 +0.5581339036209131 0.8706060914536886 0 +0.797645859957188 0.6540319953629787 0 +0.8489708810096144 0.592238192380148 0 +1.017563525876297 0.1789127675338507 0 +0.1013314713111217 1.02883569727624 0 +0.2587249332108494 0.9988030661515932 0 +0.2958529208976051 0.9888724832054925 0 +0.3746339258259126 0.9664185829581078 0 +0.572796452339944 1.488731935877131 0 +0.6628494145945181 1.452667705341444 0 +0.9633610694371049 0.3747849291976811 0 +0.4905366779988183 0.9125686318171905 0 +0.5921029838791355 0.8514116512804908 0 +0.5732340273333154 1.76577486529218 0 +1.156123794829242 0.04276959448906141 0 +0.9898377823099197 1.279066500745902 0 +0.8058700231376632 1.638384825528617 0 +0.08514165577463596 1.255436630601037 0 +1.422798763419875 0.1225205594823216 0 +1.740064106360412 0.691421420085675 0 +1.851743814737654 0.3572793215291755 0 +1.390008428226446 1.308789174308583 0 +1.216997378616021 1.444303354199705 0 +1.061407020763173 1.591079962077938 0 +1.520805830647002 1.147940677774434 0 +0.4955990443501759 1.194662749794343 0 +1.126274911554028 0.09321233201943803 0 +1.06926682174648 0.08330713735902377 0 +1.04149225279646 0.6146899442165161 0 +1.36459270471978 0.8202487006752738 0 +1.558289775852931 0.5198678972794938 0 +0.9307055116483222 0.8730445133520938 0 +0.5308516680446539 0.931735208867255 0 +0.9766863436925248 0.4464114540432426 0 +0.985662723844144 0.4946726993713682 0 +0.414036175259883 0.9511748256149621 0 +1.796976710279657 0.1737893001445605 0 +1.900043719971946 0.1600012015220335 0 +1.610037155730009 0.9606968731190765 0 +0.1618664684056866 1.90207968594875 0 +1.164649098517722 0.4800532657792035 0 +0.0940236889071869 1.19871516126472 0 +0.6948145281460931 1.088049796461561 0 +0.0375152010834033 1.063467354784323 0 +1.059358585612403 0.03551595032289604 0 +1.518113986860071 0.2982073214697316 0 +1.690709184057843 0.8739942042156701 0 +1.259011738061792 0.3495462188714486 0 +1.003382949376849 0.2622313590618652 0 +1.233505142925866 0.230068484919671 0 +1.066190662936479 1.077499598196786 0 +0.2864280833035561 1.596287773142907 0 +0.26077327063091 1.775267198504735 0 +1.189593106862306 0.9310862888106849 0 +0.7135408437861834 1.034849094039974 0 +1.00845696220035 0.3651287021091478 0 +0.4213653230938347 1.512909345401651 0 +1.585106439593642 0.7712896297129277 0 +1.660790333923994 0.7243818889490712 0 +0.742284412886948 1.33364180189666 0 +0.7410215107855009 1.4140727377124 0 +1.326576970664633 0.05730531128097838 0 +0.8235831157980494 0.9910168131759044 0 +0.3068110079566557 1.305200353061049 0 +0.9305210994349519 1.223625836221116 0 +1.221580017384016 1.332483011219602 0 +0.5725668633654364 1.112689638304212 0 +0.4315355514481705 1.179230231341497 0 +1.417904517432741 0.5393510850076291 0 +0.3558943247303962 1.268632101388988 0 +1.297793959926762 0.7013231641567605 0 +1.300954780134004 1.278867773963349 0 +0.1510651941648524 1.68536923479244 0 +0.8889964672156963 1.595650616456043 0 +0.8833598662904392 1.695444311714339 0 +0.909536629075023 0.5615142933431218 0 +1.053642885059979 0.2033729573506799 0 +1.029697577374025 0.5130284912717156 0 +1.1406115119394 1.114071494971224 0 +0.8969742676711139 1.313489510055457 0 +1.052199818882176 0.6658627411441217 0 +0.1813281516866603 1.012173044830439 0 +0.7478337546689857 0.7130498664895968 0 +1.622114319753331 0.4634157814229805 0 +0.8042651585409768 1.292480146693285 0 +0.9743975033153862 1.548052378362292 0 +1.827511523455146 0.5696956636529338 0 +1.771581292535207 0.4053731239010422 0 +0.7275410767318867 0.7814259719450345 0 +0.2497120731552854 1.03972941202752 0 +0.08049605038675402 1.061801520795098 0 +1.008220274768456 0.7564842811437653 0 +1.492703877700883 0.5702088620597213 0 +0.8355138633726098 0.663619744971877 0 +1.591345506667063 0.3295954645816138 0 +1.365928914707307 0.739348587406997 0 +1.73374455604497 0.06525951339709518 0 +0.7196909016649504 1.672139071850307 0 +1.215886897235959 1.066198606231771 0 +1.007354067214663 0.6992836776387482 0 +1.055662137902702 1.498132415782658 0 +0.2237187442932394 1.718635924960805 0 +0.6037045099571375 1.403978288770452 0 +0.6585491110769115 0.8048770286476694 0 +1.518400789228701 0.9938338610594798 0 +0.9495216573498428 0.4815541644853261 0 +0.03853360238790567 1.270450438885647 0 +1.074329866470616 0.5336156350304411 0 +1.268811266321652 0.03795889194559587 0 +0.4198180448075109 1.589667235748983 0 +0.4933120842007404 1.545930303844636 0 +0.538967854916064 1.425262529837267 0 +0.4291837015933722 1.239618534118342 0 +1.801099304389725 0.2685969430764412 0 +0.7927481149476536 0.7177947801934678 0 +1.052055872782193 0.1595935021879605 0 +0.683407628473333 1.795376523907553 0 +1.539725815698376 0.4436984762465818 0 +1.118719099629847 0.5550881608826638 0 +0.688514329539976 0.7692067094322909 0 +1.366162580904363 1.226373622084272 0 +0.9616424111947303 0.7866056310052817 0 +0.2618214878637801 1.078110328138807 0 +0.3182846993877029 1.020727142553471 0 +0.5623938234166928 1.167857086147311 0 +0.07857240423336669 1.918895610859213 0 +1.918895610955244 0.07857240429098397 0 +1.654017357740762 0.8030812769293624 0 +0.1763116981257575 1.385795244698437 0 +1.301207879384284 0.8425615943459119 0 +1.280237948723169 0.242092009140765 0 +0.9807026007574648 0.4074185226041115 0 +0.6146499526170059 1.140543930408408 0 +1.592364229032327 0.3980078810354126 0 +1.36953565993469 0.1034014392617561 0 +1.026471305598188 0.02694406440552334 0 +0.029393172717359 1.0292396403285 0 +1.622658286588352 0.1361703433408553 0 +0.2916451199586982 1.258242992904494 0 +1.325921133421859 0.2136506196154959 0 +0.07250422736588273 1.120446755321898 0 +0.1154655729773193 1.059276918931037 0 +0.2235711281631183 1.35538657963349 0 +0.3553466764541531 1.622204859043775 0 +0.4717000212756925 1.14985288958122 0 +1.123007579246225 1.4502524955585 0 +0.6652919645127456 1.052346185945562 0 +0.6500804875767984 0.8383486891140259 0 +0.5589012769463485 0.9061157451487396 0 +1.311590826718427 0.1663472298877524 0 +0.882262245143949 0.5883348696729584 0 +0.1802616987225899 1.440975599423022 0 +1.216381108894212 0.08685107339619863 0 +0.3292123389228603 1.223854312662709 0 +0.7567909726442468 1.064417997947223 0 +0.283293006156337 1.020401256289731 0 +1.09815101634351 0.2073302693932885 0 +0.8052270826204004 0.6869508848917331 0 +1.452125234238628 1.138534572481589 0 +1.15511681326342 0.5257294899909277 0 +0.9834760608468145 0.2995665970132234 0 +1.082893686830119 0.6346436515858649 0 +1.037718411101534 0.2407871325848667 0 +0.8099979111915888 1.036403360881582 0 +1.051082379174267 0.1227793741884137 0 +1.079113440302253 0.5852124409616366 0 +0.9604167480223177 0.8372405416045872 0 +1.006708562165286 0.3258497285534657 0 +0.5905386605412853 1.850125063226075 0 +1.811355723892352 0.6470680646013771 0 +0.06516204290010691 1.090073256873067 0 +0.8461288340300803 0.7609317698353109 0 +1.163122520035724 0.1237574234693646 0 +0.3291656429494038 1.797712398450035 0 +$EndNodes +$Elements +7 1214 1 1214 +1 1 1 20 +1 1 7 +2 7 8 +3 8 9 +4 9 10 +5 10 11 +6 11 12 +7 12 13 +8 13 14 +9 14 15 +10 15 16 +11 16 17 +12 17 18 +13 18 19 +14 19 20 +15 20 21 +16 21 22 +17 22 23 +18 23 24 +19 24 25 +20 25 3 +1 2 1 20 +21 3 26 +22 26 27 +23 27 28 +24 28 29 +25 29 30 +26 30 31 +27 31 32 +28 32 33 +29 33 34 +30 34 35 +31 35 36 +32 36 37 +33 37 38 +34 38 39 +35 39 40 +36 40 41 +37 41 42 +38 42 43 +39 43 44 +40 44 2 +1 3 1 14 +41 4 45 +42 45 46 +43 46 47 +44 47 48 +45 48 49 +46 49 50 +47 50 51 +48 51 52 +49 52 53 +50 53 54 +51 54 55 +52 55 56 +53 56 57 +54 57 6 +1 4 1 14 +55 6 58 +56 58 59 +57 59 60 +58 60 61 +59 61 62 +60 62 63 +61 63 64 +62 64 65 +63 65 66 +64 66 67 +65 67 68 +66 68 69 +67 69 70 +68 70 5 +1 5 1 14 +69 1 71 +70 71 72 +71 72 73 +72 73 74 +73 74 75 +74 75 76 +75 76 77 +76 77 78 +77 78 79 +78 79 80 +79 80 81 +80 81 82 +81 82 83 +82 83 4 +1 6 1 14 +83 2 84 +84 84 85 +85 85 86 +86 86 87 +87 87 88 +88 88 89 +89 89 90 +90 90 91 +91 91 92 +92 92 93 +93 93 94 +94 94 95 +95 95 96 +96 96 5 +2 1 2 1118 +97 375 524 110 +98 414 533 327 +99 203 414 356 +100 356 564 203 +101 53 473 54 +102 365 486 220 +103 455 548 322 +104 421 467 446 +105 407 414 327 +106 53 481 473 +107 90 378 89 +108 317 416 283 +109 420 442 98 +110 101 523 406 +111 322 522 455 +112 103 378 90 +113 460 559 316 +114 457 500 333 +115 416 427 283 +116 229 387 295 +117 426 474 129 +118 198 519 514 +119 412 413 248 +120 399 407 327 +121 396 566 468 +122 261 370 369 +123 408 459 125 +124 274 482 417 +125 367 521 115 +126 129 487 426 +127 239 412 248 +128 374 375 123 +129 471 531 400 +130 410 442 420 +131 329 521 367 +132 504 547 280 +133 203 488 414 +134 446 467 100 +135 514 519 355 +136 248 557 392 +137 414 488 297 +138 406 523 404 +139 125 539 408 +140 428 476 348 +141 127 461 449 +142 376 524 375 +143 198 514 499 +144 394 559 460 +145 445 566 490 +146 399 408 407 +147 392 557 376 +148 258 368 284 +149 513 533 202 +150 499 514 450 +151 133 531 471 +152 208 402 165 +153 248 392 241 +154 382 582 385 +155 396 468 301 +156 218 405 404 +157 368 369 284 +158 417 482 418 +159 51 454 107 +160 377 378 103 +161 389 547 504 +162 409 421 121 +163 397 477 324 +164 379 415 137 +165 261 369 368 +166 124 500 457 +167 297 533 414 +168 102 434 419 +169 311 530 464 +170 464 530 201 +171 324 567 397 +172 360 364 114 +173 166 373 192 +174 351 551 532 +175 51 604 454 +176 317 572 416 +177 220 536 365 +178 196 545 443 +179 538 554 27 +180 387 467 295 +181 384 538 382 +182 218 404 403 +183 381 385 270 +184 172 390 315 +185 365 578 486 +186 105 447 385 +187 434 510 419 +188 165 402 161 +189 443 545 345 +190 182 398 194 +191 193 391 390 +192 385 447 386 +193 502 547 460 +194 532 551 106 +195 227 387 229 +196 367 522 323 +197 382 385 381 +198 496 543 362 +199 112 476 428 +200 410 470 442 +201 526 534 266 +202 538 582 382 +203 432 599 462 +204 449 461 363 +205 409 470 410 +206 375 516 123 +207 54 473 125 +208 188 418 189 +209 432 462 132 +210 223 403 383 +211 137 415 135 +212 444 495 197 +213 390 391 315 +214 341 495 444 +215 256 371 261 +216 225 381 270 +217 114 448 360 +218 214 403 223 +219 391 534 526 +220 252 368 258 +221 162 393 143 +222 427 453 283 +223 179 389 147 +224 295 388 231 +225 460 547 179 +226 424 549 405 +227 393 498 325 +228 445 490 97 +229 190 394 151 +230 372 496 362 +231 448 551 351 +232 490 566 396 +233 202 533 297 +234 113 449 363 +235 411 412 237 +236 373 466 192 +237 412 558 413 +238 193 515 391 +239 370 380 369 +240 180 397 155 +241 473 481 332 +242 237 412 239 +243 135 415 92 +244 295 467 421 +245 110 516 375 +246 405 606 424 +247 388 421 409 +248 256 591 371 +249 295 421 388 +250 420 558 412 +251 184 398 182 +252 174 390 172 +253 241 392 243 +254 384 554 538 +255 168 373 166 +256 405 406 404 +257 333 580 457 +258 403 404 383 +259 195 418 188 +260 537 546 443 +261 325 577 393 +262 287 375 374 +263 127 462 461 +264 235 411 237 +265 103 379 377 +266 444 537 443 +267 448 603 551 +268 268 396 301 +269 181 395 193 +270 208 497 402 +271 130 472 471 +272 276 430 423 +273 323 518 367 +274 204 527 465 +275 369 490 284 +276 385 386 270 +277 323 522 322 +278 287 374 286 +279 92 415 91 +280 465 527 308 +281 93 135 92 +282 378 541 89 +283 94 425 116 +284 212 403 214 +285 134 474 426 +286 137 586 379 +287 130 548 455 +288 284 396 268 +289 95 425 94 +290 143 393 141 +291 389 504 162 +292 209 212 210 +293 180 477 397 +294 252 258 251 +295 250 251 249 +296 162 498 393 +297 223 382 381 +298 411 420 412 +299 222 381 225 +300 364 422 114 +301 471 548 130 +302 250 252 251 +303 192 209 208 +304 423 424 264 +305 256 261 254 +306 388 409 271 +307 161 402 159 +308 261 371 370 +309 383 384 382 +310 189 250 249 +311 477 589 489 +312 147 389 145 +313 258 284 268 +314 275 328 259 +315 189 249 187 +316 194 395 181 +317 275 350 328 +318 294 533 513 +319 298 429 277 +320 181 193 177 +321 214 223 222 +322 276 423 265 +323 213 218 212 +324 462 599 350 +325 225 270 227 +326 63 509 106 +327 104 453 427 +328 315 366 171 +329 212 214 210 +330 127 479 462 +331 218 403 212 +332 151 394 149 +333 229 295 231 +334 145 162 143 +335 209 210 208 +336 246 287 286 +337 246 286 285 +338 149 179 147 +339 141 217 139 +340 286 374 349 +341 286 349 348 +342 223 381 222 +343 68 422 129 +344 166 192 165 +345 231 388 233 +346 153 190 151 +347 265 423 264 +348 270 387 227 +349 192 208 165 +350 233 271 235 +351 141 393 217 +352 277 429 276 +353 350 550 328 +354 358 452 378 +355 258 268 260 +356 145 389 162 +357 418 482 189 +358 337 355 339 +359 157 180 155 +360 67 422 68 +361 159 402 401 +362 268 301 269 +363 155 397 153 +364 149 394 179 +365 239 248 241 +366 406 517 101 +367 172 315 171 +368 267 276 265 +369 233 388 271 +370 254 368 252 +371 178 181 177 +372 187 398 184 +373 360 448 346 +374 182 194 181 +375 216 417 205 +376 164 165 161 +377 350 599 550 +378 409 410 271 +379 79 427 416 +380 243 246 245 +381 385 582 105 +382 177 390 174 +383 417 574 565 +384 74 449 113 +385 153 397 190 +386 263 264 219 +387 426 507 425 +388 223 383 382 +389 362 543 113 +390 170 171 168 +391 386 447 435 +392 121 470 409 +393 253 254 252 +394 190 559 394 +395 282 298 277 +396 188 189 187 +397 429 430 276 +398 341 444 343 +399 205 418 195 +400 116 135 93 +401 255 256 254 +402 159 401 157 +403 287 392 376 +404 257 259 256 +405 209 213 212 +406 115 522 367 +407 206 450 289 +408 272 275 259 +409 157 401 180 +410 179 547 389 +411 278 296 216 +412 171 373 168 +413 391 526 315 +414 114 422 67 +415 299 419 298 +416 421 446 121 +417 215 216 205 +418 106 509 451 +419 113 587 362 +420 136 137 135 +421 176 177 174 +422 186 187 184 +423 243 392 246 +424 286 348 285 +425 430 528 423 +426 346 448 351 +427 213 219 218 +428 356 414 407 +429 138 139 137 +430 478 525 131 +431 246 285 245 +432 408 593 459 +433 363 607 352 +434 511 591 328 +435 140 141 139 +436 261 368 254 +437 142 143 141 +438 164 166 165 +439 370 483 380 +440 214 222 221 +441 337 514 355 +442 258 260 251 +443 56 108 57 +444 51 107 52 +445 63 106 64 +446 59 111 60 +447 218 606 405 +448 200 205 195 +449 398 594 194 +450 144 145 143 +451 416 572 117 +452 267 277 276 +453 46 130 47 +454 68 129 69 +455 263 265 264 +456 66 114 67 +457 48 115 49 +458 222 225 224 +459 249 398 187 +460 313 356 331 +461 363 587 113 +462 209 466 213 +463 123 576 374 +464 464 564 356 +465 6 124 58 +466 61 126 62 +467 96 134 95 +468 82 133 83 +469 430 437 99 +470 146 147 145 +471 225 227 226 +472 99 528 430 +473 54 125 55 +474 424 606 264 +475 170 172 171 +476 282 299 298 +477 419 510 298 +478 289 450 292 +479 268 269 260 +480 94 116 93 +481 80 117 81 +482 343 443 345 +483 180 589 477 +484 148 149 147 +485 91 415 103 +486 213 466 262 +487 134 425 95 +488 287 376 375 +489 80 416 117 +490 348 575 428 +491 79 416 80 +492 178 182 181 +493 227 229 228 +494 193 390 177 +495 91 103 90 +496 77 104 78 +497 88 112 87 +498 74 113 75 +499 86 128 85 +500 72 127 73 +501 216 574 417 +502 249 475 398 +503 328 550 511 +504 150 151 149 +505 247 431 358 +506 271 411 235 +507 311 464 313 +508 335 337 336 +509 462 479 132 +510 417 418 205 +511 250 253 252 +512 337 339 338 +513 176 178 177 +514 400 531 117 +515 229 231 230 +516 152 153 151 +517 343 444 443 +518 259 591 256 +519 394 460 179 +520 339 341 340 +521 517 549 119 +522 112 428 87 +523 292 335 334 +524 275 461 350 +525 410 420 411 +526 231 233 232 +527 292 450 335 +528 154 155 153 +529 366 373 171 +530 87 428 86 +531 86 428 128 +532 461 462 350 +533 405 549 406 +534 395 596 193 +535 108 456 57 +536 186 188 187 +537 59 457 111 +538 111 458 60 +539 56 459 108 +540 215 278 216 +541 341 343 342 +542 233 235 234 +543 494 495 355 +544 130 455 47 +545 450 514 335 +546 48 455 115 +547 200 215 205 +548 253 255 254 +549 255 257 256 +550 257 272 259 +551 116 136 135 +552 6 456 124 +553 386 387 270 +554 156 157 155 +555 124 457 58 +556 61 458 126 +557 355 495 339 +558 452 541 378 +559 235 237 236 +560 246 392 287 +561 352 587 363 +562 45 472 46 +563 260 475 251 +564 125 459 55 +565 343 345 344 +566 158 159 157 +567 356 407 331 +568 237 239 238 +569 136 138 137 +570 345 347 346 +571 46 472 130 +572 299 469 419 +573 160 161 159 +574 138 140 139 +575 214 221 210 +576 347 361 360 +577 239 241 240 +578 140 142 141 +579 163 164 161 +580 77 496 104 +581 167 168 166 +582 423 528 424 +583 361 365 364 +584 128 478 85 +585 72 479 127 +586 241 243 242 +587 142 144 143 +588 431 452 358 +589 372 584 296 +590 243 245 244 +591 358 378 377 +592 359 475 260 +593 452 476 112 +594 219 606 218 +595 222 224 221 +596 313 464 356 +597 144 146 145 +598 98 558 420 +599 134 426 425 +600 262 263 219 +601 326 327 294 +602 33 446 32 +603 117 572 400 +604 266 267 265 +605 169 170 168 +606 22 437 21 +607 192 466 209 +608 225 226 224 +609 146 148 147 +610 16 445 15 +611 36 442 35 +612 39 440 38 +613 247 358 357 +614 31 435 30 +615 18 433 17 +616 339 495 341 +617 78 427 79 +618 50 604 51 +619 278 453 296 +620 404 523 383 +621 293 294 290 +622 73 449 74 +623 148 150 149 +624 335 336 334 +625 82 531 133 +626 227 228 226 +627 250 482 253 +628 431 476 452 +629 413 590 524 +630 337 338 336 +631 173 174 172 +632 369 602 490 +633 281 282 277 +634 339 340 338 +635 150 152 151 +636 410 411 271 +637 320 323 322 +638 478 605 525 +639 229 230 228 +640 341 342 340 +641 301 512 269 +642 152 154 153 +643 426 536 507 +644 183 184 182 +645 231 232 230 +646 469 512 301 +647 284 490 396 +648 175 176 174 +649 47 455 48 +650 57 456 6 +651 58 457 59 +652 233 234 232 +653 60 458 61 +654 154 156 155 +655 55 459 56 +656 343 344 342 +657 129 474 69 +658 296 569 372 +659 311 313 312 +660 304 311 306 +661 199 200 195 +662 309 310 308 +663 251 475 249 +664 307 308 305 +665 185 186 184 +666 283 453 278 +667 235 236 234 +668 156 158 157 +669 310 465 308 +670 329 367 310 +671 103 415 379 +672 345 346 344 +673 313 331 330 +674 191 195 188 +675 285 431 245 +676 279 283 278 +677 164 167 166 +678 347 360 346 +679 158 160 159 +680 273 352 272 +681 124 506 500 +682 237 238 236 +683 253 482 274 +684 518 552 465 +685 365 536 487 +686 353 362 352 +687 239 240 238 +688 213 262 219 +689 361 364 360 +690 85 478 84 +691 71 479 72 +692 160 163 161 +693 21 437 120 +694 76 496 77 +695 32 446 100 +696 321 400 319 +697 8 132 7 +698 44 131 43 +699 241 242 240 +700 17 433 118 +701 468 469 301 +702 170 173 172 +703 294 513 290 +704 504 588 573 +705 293 326 294 +706 221 497 210 +707 440 441 38 +708 435 447 30 +709 15 445 97 +710 167 169 168 +711 326 399 327 +712 442 470 35 +713 19 102 18 +714 26 101 3 +715 32 100 31 +716 23 99 22 +717 243 244 242 +718 43 439 42 +719 9 432 8 +720 5 560 96 +721 83 561 4 +722 425 507 116 +723 107 481 52 +724 245 247 244 +725 263 266 265 +726 69 474 70 +727 245 431 247 +728 267 281 277 +729 189 482 250 +730 10 122 9 +731 42 123 41 +732 358 377 357 +733 289 292 291 +734 40 110 39 +735 12 109 11 +736 24 436 23 +737 20 434 19 +738 37 98 36 +739 15 97 14 +740 282 300 299 +741 406 549 517 +742 11 438 10 +743 25 119 24 +744 21 120 20 +745 292 334 333 +746 367 518 310 +747 38 441 37 +748 17 118 16 +749 35 470 34 +750 269 359 260 +751 34 121 33 +752 475 594 398 +753 30 447 29 +754 441 558 98 +755 373 601 466 +756 29 105 28 +757 178 183 182 +758 117 531 81 +759 206 499 450 +760 52 481 53 +761 348 476 285 +762 407 408 331 +763 106 551 64 +764 173 175 174 +765 215 279 278 +766 518 568 552 +767 175 201 176 +768 112 541 452 +769 99 437 22 +770 563 577 357 +771 100 435 31 +772 102 433 18 +773 131 439 43 +774 8 432 132 +775 307 309 308 +776 164 485 167 +777 104 427 78 +778 136 220 138 +779 199 207 200 +780 311 312 306 +781 309 329 310 +782 320 529 323 +783 186 191 188 +784 110 440 39 +785 13 483 12 +786 98 442 36 +787 313 330 312 +788 121 446 33 +789 127 449 73 +790 118 445 16 +791 504 573 498 +792 140 211 142 +793 257 273 272 +794 191 199 195 +795 253 274 255 +796 281 463 282 +797 42 439 123 +798 122 432 9 +799 23 436 99 +800 210 497 208 +801 19 434 102 +802 298 510 429 +803 285 476 431 +804 62 509 63 +805 183 185 184 +806 419 540 102 +807 114 603 448 +808 66 603 114 +809 144 196 146 +810 10 438 122 +811 314 497 221 +812 321 471 400 +813 109 438 11 +814 279 317 283 +815 119 436 24 +816 120 434 20 +817 309 454 329 +818 342 451 340 +819 417 565 274 +820 170 488 173 +821 282 463 300 +822 37 441 98 +823 150 197 152 +824 291 293 290 +825 330 332 312 +826 206 288 163 +827 108 593 399 +828 289 290 288 +829 191 204 199 +830 163 485 164 +831 292 333 291 +832 224 314 221 +833 320 322 321 +834 273 353 352 +835 346 351 344 +836 88 541 112 +837 353 372 362 +838 113 543 75 +839 399 555 108 +840 279 318 317 +841 156 198 158 +842 456 506 124 +843 34 470 121 +844 304 306 305 +845 306 307 305 +846 413 557 248 +847 318 320 319 +848 12 483 109 +849 304 305 303 +850 169 488 170 +851 302 304 303 +852 160 206 163 +853 247 357 325 +854 304 530 311 +855 330 473 332 +856 230 324 228 +857 433 468 118 +858 49 521 50 +859 29 447 105 +860 308 527 305 +861 173 203 175 +862 167 202 169 +863 202 297 169 +864 234 316 232 +865 302 303 185 +866 524 557 413 +867 345 545 347 +868 183 302 185 +869 28 538 27 +870 320 321 319 +871 240 280 238 +872 247 325 244 +873 289 291 290 +874 201 505 176 +875 485 513 202 +876 352 607 272 +877 76 543 496 +878 273 354 353 +879 422 608 129 +880 206 289 288 +881 318 319 317 +882 455 522 115 +883 275 607 461 +884 200 480 215 +885 207 480 200 +886 335 514 337 +887 126 509 62 +888 1 570 71 +889 84 571 2 +890 167 485 202 +891 471 472 133 +892 7 570 1 +893 2 571 44 +894 176 505 178 +895 27 554 26 +896 315 526 366 +897 101 517 3 +898 443 546 196 +899 123 516 41 +900 194 553 395 +901 40 516 110 +902 186 503 191 +903 226 489 224 +904 458 520 126 +905 197 537 444 +906 288 513 485 +907 25 517 119 +908 185 503 186 +909 191 503 204 +910 386 583 387 +911 451 509 508 +912 496 569 104 +913 70 560 5 +914 4 561 45 +915 428 575 128 +916 138 486 140 +917 173 488 203 +918 116 507 136 +919 97 602 595 +920 293 506 326 +921 300 512 299 +922 152 494 154 +923 142 491 144 +924 255 484 257 +925 207 529 480 +926 288 485 163 +927 383 523 384 +928 139 586 137 +929 461 607 363 +930 41 516 40 +931 215 480 279 +932 340 508 338 +933 3 517 25 +934 120 510 434 +935 355 519 494 +936 238 502 236 +937 178 505 183 +938 397 567 190 +939 232 501 230 +940 266 534 267 +941 380 602 369 +942 115 521 49 +943 297 488 169 +944 220 486 138 +945 109 511 438 +946 400 572 319 +947 158 499 160 +948 307 493 309 +949 371 591 511 +950 201 530 505 +951 457 580 111 +952 136 507 220 +953 309 493 454 +954 211 491 142 +955 528 592 424 +956 140 486 211 +957 306 492 307 +958 312 492 306 +959 129 608 487 +960 332 492 312 +961 274 484 255 +962 162 504 498 +963 257 484 273 +964 303 503 185 +965 279 480 318 +966 144 491 196 +967 468 540 469 +968 299 512 469 +969 291 500 293 +970 366 601 373 +971 224 489 314 +972 197 494 152 +973 451 508 340 +974 433 540 468 +975 197 495 494 +976 374 576 349 +977 111 535 458 +978 370 597 483 +979 454 493 107 +980 361 578 365 +981 81 531 82 +982 333 500 291 +983 160 499 206 +984 183 505 302 +985 198 499 158 +986 500 506 293 +987 280 502 238 +988 230 501 324 +989 273 484 354 +990 371 597 370 +991 387 583 467 +992 508 509 126 +993 316 501 232 +994 492 493 307 +995 362 587 352 +996 480 529 318 +997 290 513 288 +998 316 579 460 +999 465 552 204 +1000 408 539 331 +1001 154 519 156 +1002 372 569 496 +1003 503 527 204 +1004 338 520 336 +1005 156 519 198 +1006 131 525 439 +1007 97 595 14 +1008 393 577 217 +1009 565 574 354 +1010 201 564 464 +1011 310 518 465 +1012 327 533 294 +1013 263 556 266 +1014 303 527 503 +1015 281 515 463 +1016 342 532 451 +1017 13 595 483 +1018 451 532 106 +1019 347 544 361 +1020 351 532 344 +1021 494 519 154 +1022 505 530 302 +1023 436 528 99 +1024 217 577 563 +1025 508 520 338 +1026 75 543 76 +1027 439 576 123 +1028 264 606 219 +1029 89 541 88 +1030 330 539 473 +1031 429 585 430 +1032 110 524 440 +1033 150 537 197 +1034 489 581 477 +1035 196 546 146 +1036 453 569 296 +1037 483 597 109 +1038 454 604 329 +1039 377 563 357 +1040 402 598 401 +1041 477 581 324 +1042 460 579 502 +1043 467 583 100 +1044 401 589 180 +1045 318 529 320 +1046 305 527 303 +1047 296 574 216 +1048 321 548 471 +1049 204 552 199 +1050 302 530 304 +1051 399 593 408 +1052 511 550 438 +1053 458 535 520 +1054 108 555 456 +1055 344 532 342 +1056 217 586 139 +1057 469 540 419 +1058 524 590 440 +1059 498 573 325 +1060 175 564 201 +1061 325 573 244 +1062 266 556 526 +1063 336 535 334 +1064 267 534 281 +1065 595 602 380 +1066 506 555 326 +1067 520 535 336 +1068 45 561 472 +1069 365 608 364 +1070 148 537 150 +1071 473 539 125 +1072 459 593 108 +1073 324 581 228 +1074 430 585 437 +1075 234 579 316 +1076 349 575 348 +1077 329 604 521 +1078 328 591 259 +1079 105 582 538 +1080 589 598 314 +1081 64 551 65 +1082 281 534 515 +1083 226 581 489 +1084 487 608 365 +1085 487 536 426 +1086 349 605 575 +1087 364 608 422 +1088 240 588 280 +1089 331 539 330 +1090 379 563 377 +1091 280 588 504 +1092 456 555 506 +1093 529 568 323 +1094 120 585 510 +1095 502 579 236 +1096 354 584 353 +1097 483 595 380 +1098 525 605 349 +1099 468 566 118 +1100 105 538 28 +1101 269 542 359 +1102 128 605 478 +1103 126 520 508 +1104 146 546 148 +1105 134 560 474 +1106 102 540 433 +1107 111 580 535 +1108 507 536 220 +1109 438 550 122 +1110 436 592 528 +1111 322 548 321 +1112 148 546 537 +1113 262 556 263 +1114 100 583 435 +1115 272 607 275 +1116 481 562 332 +1117 104 569 453 +1118 96 560 134 +1119 133 561 83 +1120 474 560 70 +1121 26 554 101 +1122 199 552 207 +1123 551 603 65 +1124 510 585 429 +1125 535 580 334 +1126 326 555 399 +1127 65 603 66 +1128 101 554 523 +1129 497 598 402 +1130 478 571 84 +1131 71 570 479 +1132 300 542 512 +1133 472 561 133 +1134 491 545 196 +1135 211 544 491 +1136 194 594 553 +1137 300 600 542 +1138 521 604 50 +1139 515 534 391 +1140 490 602 97 +1141 131 571 478 +1142 479 570 132 +1143 550 599 122 +1144 466 601 262 +1145 542 600 553 +1146 280 547 502 +1147 132 570 7 +1148 44 571 131 +1149 526 556 366 +1150 107 562 481 +1151 523 554 384 +1152 401 598 589 +1153 319 572 317 +1154 332 562 492 +1155 14 595 13 +1156 574 584 354 +1157 211 578 544 +1158 489 589 314 +1159 203 564 175 +1160 544 578 361 +1161 511 597 371 +1162 244 573 242 +1163 118 566 445 +1164 349 576 525 +1165 236 579 234 +1166 334 580 333 +1167 492 562 493 +1168 228 581 226 +1169 359 594 475 +1170 437 585 120 +1171 316 559 501 +1172 512 542 269 +1173 357 577 325 +1174 553 594 359 +1175 435 583 386 +1176 440 590 441 +1177 552 568 207 +1178 193 596 515 +1179 274 565 484 +1180 556 601 366 +1181 544 545 491 +1182 119 592 436 +1183 323 568 518 +1184 353 584 372 +1185 122 599 432 +1186 242 588 240 +1187 573 588 242 +1188 463 600 300 +1189 542 553 359 +1190 484 565 354 +1191 501 567 324 +1192 493 562 107 +1193 596 600 463 +1194 486 578 211 +1195 207 568 529 +1196 262 601 556 +1197 525 576 439 +1198 376 557 524 +1199 314 598 497 +1200 109 597 511 +1201 563 586 217 +1202 395 600 596 +1203 347 545 544 +1204 559 567 501 +1205 515 596 463 +1206 575 605 128 +1207 190 567 559 +1208 549 592 119 +1209 441 590 558 +1210 296 584 574 +1211 424 592 549 +1212 379 586 563 +1213 553 600 395 +1214 558 590 413 +$EndElements diff --git a/examples/scripts/benchmark_solver_comparison.py b/examples/scripts/benchmark_solver_comparison.py deleted file mode 100644 index 850204c..0000000 --- a/examples/scripts/benchmark_solver_comparison.py +++ /dev/null @@ -1,97 +0,0 @@ -import lizzy as liz -import time - -# Setup model -def setup_model(): - model = liz.LizzyModel() - model.read_mesh_file("../meshes/Rect1M_R2.msh") - model.assign_simulation_parameters(mu=0.1, wo_delta_time=100, fill_tolerance=0.01) - model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "domain_material") - model.assign_material("domain_material", 'domain') - model.create_sensor(0.2, 0.25, 0) - model.create_inlet(100000, "inlet_left") - model.assign_inlet("inlet_left", "left_edge") - return model - -# Test a solver configuration -def test_solver_config(solver_type, use_masked, name, **kwargs): - print(f"\nTesting {name}...") - model = setup_model() - - try: - start = time.time() - model.initialise_solver(solver_type=solver_type, use_masked_solver=use_masked, **kwargs) - solution = model.solve() - elapsed = time.time() - start - print(f" ✓ Success: {elapsed:.2f}s") - print(f" - Fill time: {solution.time[-1]:.5f}s") - return True, elapsed, solution - except Exception as e: - print(f" ✗ Failed: {e}") - import traceback - traceback.print_exc() - return False, None, None - -# Define solver configurations to test -solver_configs = [ - (liz.SolverType.ITERATIVE_PETSC, True, "Masked Solver (PETSc CG+GAMG)", {"solver_tol": 1e-4, "solver_max_iter": 1000, "ksp_type": "cg", "pc_type": "gamg"}), - (liz.SolverType.ITERATIVE_PETSC, True, "Masked Solver (PETSc GMRES+ILU)", {"solver_tol": 1e-4, "solver_max_iter": 1000, "ksp_type": "gmres", "pc_type": "ilu"}), - (liz.SolverType.DIRECT_SPARSE, True, "Masked Solver (DIRECT_SPARSE)", {}), - (liz.SolverType.DIRECT_DENSE, True, "Masked Solver (DIRECT_DENSE)", {}), - (liz.SolverType.DIRECT_SPARSE, False, "Traditional Solver (DIRECT_SPARSE)", {}), -] - -# Run tests -print("=" * 70) -print("Solver Performance Comparison: Traditional vs Masked (Optimized)") -print("=" * 70) - -results = [] -baseline_time = None - -for solver_type, use_masked, name, kwargs in solver_configs: - success, time_elapsed, solution = test_solver_config(solver_type, use_masked, name, **kwargs) - if success: - # Store baseline time from traditional solver - if not use_masked and baseline_time is None: - baseline_time = time_elapsed - - results.append((name, time_elapsed, use_masked)) - -# Summary -print("\n" + "=" * 70) -print("Performance Summary") -print("=" * 70) - -if results: - # Sort by time - results.sort(key=lambda x: x[1]) - - print(f"\n{'Rank':<6} {'Solver Configuration':<45} {'Time (s)':<12} {'Speedup':<10}") - print("-" * 70) - - for i, (name, t, is_masked) in enumerate(results, 1): - if baseline_time and baseline_time > 0: - speedup = baseline_time / t - speedup_str = f"{speedup:.2f}x" - else: - speedup_str = "N/A" - - marker = "🚀" if is_masked else "📊" - print(f"{i:<6} {marker} {name:<43} {t:<12.2f} {speedup_str:<10}") - - # Highlight the winner - if len(results) > 1: - winner_name, winner_time, winner_masked = results[0] - print("\n" + "=" * 70) - print(f"🏆 Fastest: {winner_name}") - if baseline_time and winner_masked: - improvement = ((baseline_time - winner_time) / baseline_time) * 100 - print(f" Performance improvement over traditional: {improvement:.1f}%") - print("=" * 70) -else: - print("No solvers succeeded") - -print("\nLegend:") -print(" 📊 = Traditional solver (full matrix)") -print(" 🚀 = Masked/optimized solver (submatrix extraction)") diff --git a/examples/scripts/channel_infusion.py b/examples/scripts/channel_infusion.py index e79a155..5625d70 100644 --- a/examples/scripts/channel_infusion.py +++ b/examples/scripts/channel_infusion.py @@ -1,16 +1,19 @@ -import lizzy as liz +import lizzy -model = liz.LizzyModel() +model = lizzy.LizzyModel() model.read_mesh_file("../meshes/Rect1M_R1.msh") -model.assign_simulation_parameters(mu=0.1, wo_delta_time=100, fill_tolerance=0.01) +model.assign_simulation_parameters(wo_delta_time=100) -model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "domain_material") +model.create_resin("resin", 0.1) +model.assign_resin("resin") + +model.create_material("domain_material", (1E-10, 1E-10, 1E-10), 0.5, 1.0) model.assign_material("domain_material", 'domain') -model.create_inlet(100000, "inlet_left") +model.create_inlet("inlet_left", 100000) model.assign_inlet("inlet_left", "left_edge") model.initialise_solver() -solution = model.solve() +model.solve() -model.save_results(solution, "Rect1M_R1") \ No newline at end of file +model.save_results() \ No newline at end of file diff --git a/examples/scripts/channel_infusion_each_solver_test.py b/examples/scripts/channel_infusion_each_solver_test.py deleted file mode 100644 index b34264f..0000000 --- a/examples/scripts/channel_infusion_each_solver_test.py +++ /dev/null @@ -1,62 +0,0 @@ -import lizzy as liz -import time - -# Setup model -def setup_model(): - model = liz.LizzyModel() - model.read_mesh_file("../meshes/Rect1M_R2.msh") - model.assign_simulation_parameters(mu=0.1, wo_delta_time=100, fill_tolerance=0.01) - model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "domain_material") - model.assign_material("domain_material", 'domain') - model.create_sensor(0.2, 0.25, 0) - model.create_inlet(100000, "inlet_left") - model.assign_inlet("inlet_left", "left_edge") - return model - -# Test a solver -def test_solver(solver_type, name, **kwargs): - print(f"\nTesting {name}...") - model = setup_model() - - try: - start = time.time() - model.initialise_solver(solver_type=solver_type, **kwargs) - model.solve() - elapsed = time.time() - start - print(f" ✓ Success: {elapsed:.2f}s") - return True, elapsed - except Exception as e: - print(f" ✗ Failed: {e}") - return False, None - -# Define solvers to test -solvers = [ - (liz.SolverType.DIRECT_DENSE, "Direct Dense", {}), - (liz.SolverType.DIRECT_SPARSE, "Direct Sparse", {}), - (liz.SolverType.ITERATIVE_PETSC, "PETSc CG+GAMG", - {"solver_tol": 1e-8, "solver_max_iter": 1000, "ksp_type": "cg", "pc_type": "gamg"}), - (liz.SolverType.ITERATIVE_PETSC, "PETSc GMRES+ILU", - {"solver_tol": 1e-8, "solver_max_iter": 1000, "ksp_type": "gmres", "pc_type": "ilu"}), -] - -# Run tests -print("=" * 50) -print("Solver Comparison Test") -print("=" * 50) - -results = [] -for solver_type, name, kwargs in solvers: - success, time_elapsed = test_solver(solver_type, name, **kwargs) - if success: - results.append((name, time_elapsed)) - -# Summary -print("\n" + "=" * 50) -print("Summary") -print("=" * 50) -if results: - results.sort(key=lambda x: x[1]) - for i, (name, t) in enumerate(results, 1): - print(f"{i}. {name}: {t:.2f}s") -else: - print("No solvers succeeded") \ No newline at end of file diff --git a/examples/scripts/complex.py b/examples/scripts/complex.py index be5095f..a3fee2b 100644 --- a/examples/scripts/complex.py +++ b/examples/scripts/complex.py @@ -2,21 +2,25 @@ model = liz.LizzyModel() model.read_mesh_file("../meshes/Complex_rotated.msh") -model.assign_simulation_parameters(mu=0.1, wo_delta_time=100, fill_tolerance=0.01) +model.assign_simulation_parameters(wo_delta_time=100, fill_tolerance=0.01) -model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "material_iso") -model.create_material(1E-10, 1E-11, 1E-11, 0.5, 1.0, "material_aniso") -model.create_material(1E-7, 1E-7, 1E-7, 0.5, 0.5, "material_racetrack") -rosette_ramp = model.create_rosette(model._mesh.nodes[12].coords, model._mesh.nodes[13].coords, "rosette_ramp") +model.create_resin("resin", 0.1) +model.assign_resin("resin") + +model.create_material("material_iso", (1E-10, 1E-10, 1E-10), 0.5, 1.0, ) +model.create_material("material_aniso", (1E-10, 1E-11, 1E-11), 0.5, 1.0, ) +model.create_material("material_racetrack", (1E-7, 1E-7, 1E-7), 0.5, 0.5, ) + +rosette_ramp = model.create_rosette("rosette_ramp", (model._mesh.nodes[12].coords - model._mesh.nodes[13].coords) ) model.assign_material("material_iso", 'Lshape') model.assign_material("material_aniso", 'ramp', rosette_ramp) model.assign_material("material_racetrack", 'racetrack') -model.create_inlet(100000, "inlet_left") +model.create_inlet("inlet_left", 100000) model.assign_inlet("inlet_left", "inlet") model.initialise_solver() -solution = model.solve() +model.solve() -model.save_results(solution, "Complex_rotated") +model.save_results() diff --git a/examples/scripts/multi_materials.py b/examples/scripts/multi_materials.py index 440fa58..6e9683a 100644 --- a/examples/scripts/multi_materials.py +++ b/examples/scripts/multi_materials.py @@ -2,18 +2,21 @@ model = liz.LizzyModel() model.read_mesh_file("../meshes/Triforce_R1.msh") -model.assign_simulation_parameters(mu=0.1, wo_delta_time=100) +model.assign_simulation_parameters(wo_delta_time=100) -model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "high_perm_material") -model.create_material(1E-13, 1E-13, 1E-13, 0.5, 1.0, "low_perm_material") +model.create_resin("resin", 0.1) +model.assign_resin("resin") + +model.create_material("high_perm_material", (1E-10, 1E-10, 1E-10), 0.5, 1.0, ) +model.create_material("low_perm_material", (1E-13, 1E-13, 1E-13), 0.5, 1.0, ) model.assign_material("high_perm_material", 'background') model.assign_material("low_perm_material", 'triforce') -model.create_inlet(1E+05, "inner_inlet") +model.create_inlet("inner_inlet", 1E+05) model.assign_inlet("inner_inlet", "inner_rim") model.initialise_solver() -solution = model.solve() +model.solve() -model.save_results(solution, "Triforce") \ No newline at end of file +model.save_results() \ No newline at end of file diff --git a/examples/scripts/radial_infusion.py b/examples/scripts/radial_infusion.py index 2400f5c..0615a55 100644 --- a/examples/scripts/radial_infusion.py +++ b/examples/scripts/radial_infusion.py @@ -2,16 +2,19 @@ model = liz.LizzyModel() model.read_mesh_file("../meshes/Radial.msh") -model.assign_simulation_parameters(mu=0.1, wo_delta_time=100) +model.assign_simulation_parameters(wo_delta_time=500) -rosette = model.create_rosette((1,1,0)) -model.create_material(1E-10, 1E-11, 1E-10, 0.5, 1.0, "aniso_material") -model.assign_material("aniso_material", 'domain', rosette) +model.create_resin("resin", 0.1) +model.assign_resin("resin") -model.create_inlet(1E+05, "inner_inlet") +rosette = model.create_rosette("rosette", (1,1,0)) +model.create_material("aniso_material", (1E-10, 1E-10, 1E-10), 0.5, 1.0) +model.assign_material("aniso_material", 'domain', "rosette") + +model.create_inlet("inner_inlet", 1E+05) model.assign_inlet("inner_inlet", "inner_rim") model.initialise_solver() -solution = model.solve() +model.solve() -model.save_results(solution, "Anisotropic_Radial") \ No newline at end of file +model.save_results() \ No newline at end of file diff --git a/examples/scripts/radial_infusion_annulus.py b/examples/scripts/radial_infusion_annulus.py new file mode 100644 index 0000000..cc74ce9 --- /dev/null +++ b/examples/scripts/radial_infusion_annulus.py @@ -0,0 +1,26 @@ +import lizzy as liz +from lizzy import SolverType + +model = liz.LizzyModel() + +model.read_mesh_file("../meshes/quarter_annulus.msh") +model.assign_simulation_parameters(wo_delta_time=100, fill_tolerance=0.00) + +model.create_resin("resin", 0.1) +model.assign_resin("resin") + +model.create_material("domain_material", (1E-10, 1E-10, 1E-10), 0.6, 1.0) +model.assign_material("domain_material", 'domain') + +model.create_inlet("inner_radius", 100000) +model.assign_inlet("inner_radius", "inlet") + +model.initialise_solver() +model.solve() + +model.save_results() + +# Lee, Y. M., et al. "Analysis of flow in the RTM Process." SAE transactions (1989): 65-75. +# Eq. (A-4) : +# r_0 = 1, r_m = 2, phi = 0.6, eta = 0.1, K = 1e-10, P0 = 1e5 +# Then the filling time = 3817.77 seconds \ No newline at end of file diff --git a/examples/scripts/solver_comparison.py b/examples/scripts/solver_comparison.py new file mode 100644 index 0000000..38c464e --- /dev/null +++ b/examples/scripts/solver_comparison.py @@ -0,0 +1,116 @@ +""" +Solver Comparison Example +========================= + +This example demonstrates how to use different solvers available in Lizzy. +The same filling simulation is run with three different solver types to +compare their results. + +Available solvers: +- DIRECT_SPARSE: Direct solver using sparse matrix factorization (recommended for small/medium meshes) +- DIRECT_DENSE: Direct solver using dense matrices (only for very small problems) +- ITERATIVE_PETSC: Iterative solver using PETSc library (recommended for large meshes) +""" + +import time +import lizzy as liz + +# ============================================================================= +# Model Setup (common for all solvers) +# ============================================================================= + +# Read mesh +model = liz.LizzyModel() +model.read_mesh_file("../meshes/Rect1M_R1.msh") + +# Set simulation parameters +model.assign_simulation_parameters( + wo_delta_time=100, # Initial time step guess [s] + fill_tolerance=0.00 # Fill fraction tolerance for CV to be considered full +) + +# Create and assign resin +model.create_resin("resin", 0.1) # Resin viscosity 0.1 [Pa.s] +model.assign_resin("resin") + +# Create and assign material +model.create_material( + "glass_fiber", # Material name + (1e-10, 1e-10, 1e-10), # Permeability (k1, k2, k3) [m^2] + 0.5, # Porosity [-] + 1.0 # Thickness [m] +) +model.assign_material("glass_fiber", "domain") + +# Create inlet with prescribed pressure +model.create_inlet("left_inlet", 1e5) # 1 bar +model.assign_inlet("left_inlet", "left_edge") + +# ============================================================================= +# Example 1: Direct Sparse Solver (default, recommended) +# ============================================================================= +print("\n" + "="*60) +print("Running with DIRECT_SPARSE solver") +print("="*60) + +model.initialise_solver(solver_type=liz.SolverType.DIRECT_SPARSE) +start_time = time.time() +solution_sparse = model.solve() +time_sparse = time.time() - start_time + +print(f"Fill time: {solution_sparse.time[-1]:.2f} s") +print(f"Number of time steps: {len(solution_sparse.time)}") +print(f"Solver runtime: {time_sparse:.2f} s") + +# ============================================================================= +# Example 2: Direct Dense Solver (for comparison only - slow for large meshes) +# ============================================================================= +print("\n" + "="*60) +print("Running with DIRECT_DENSE solver") +print("="*60) + +model.initialise_solver(solver_type=liz.SolverType.DIRECT_DENSE) +start_time = time.time() +solution_dense = model.solve() +time_dense = time.time() - start_time + +print(f"Fill time: {solution_dense.time[-1]:.2f} s") +print(f"Number of time steps: {len(solution_dense.time)}") +print(f"Solver runtime: {time_dense:.2f} s") + +# ============================================================================= +# Example 3: Iterative PETSc Solver (best for large meshes) +# ============================================================================= +print("\n" + "="*60) +print("Running with ITERATIVE_PETSC solver") +print("="*60) + +model.initialise_solver( + solver_type=liz.SolverType.ITERATIVE_PETSC, + solver_tol=1e-6, # Solver tolerance + solver_max_iter=1000, # Maximum iterations + ksp_type="cg", # Conjugate gradient method + pc_type="gamg" # Algebraic multigrid preconditioner +) +start_time = time.time() +solution_petsc = model.solve() +time_petsc = time.time() - start_time + +print(f"Fill time: {solution_petsc.time[-1]:.2f} s") +print(f"Number of time steps: {len(solution_petsc.time)}") +print(f"Solver runtime: {time_petsc:.2f} s") + +# ============================================================================= +# Compare Results +# ============================================================================= +print("\n" + "="*60) +print("Results Comparison") +print("="*60) +print(f"{'Solver':<20} {'Fill Time [s]':<15} {'Runtime [s]':<15}") +print("-"*50) +print(f"{'DIRECT_SPARSE':<20} {solution_sparse.time[-1]:<15.4f} {time_sparse:<15.2f}") +print(f"{'DIRECT_DENSE':<20} {solution_dense.time[-1]:<15.4f} {time_dense:<15.2f}") +print(f"{'ITERATIVE_PETSC':<20} {solution_petsc.time[-1]:<15.4f} {time_petsc:<15.2f}") + +# Save results from one of the solutions +model.save_results(solution_sparse, "solver_comparison") diff --git a/src/lizzy/_core/bcond/__init__.py b/src/lizzy/_core/bcond/__init__.py index 053b794..19f9a29 100644 --- a/src/lizzy/_core/bcond/__init__.py +++ b/src/lizzy/_core/bcond/__init__.py @@ -4,6 +4,5 @@ # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with this program. If not, see . -from .manager import BCManager +from .manager import GatesManager from .gates import Inlet -from .gates import SolverBCs \ No newline at end of file diff --git a/src/lizzy/_core/bcond/gates.py b/src/lizzy/_core/bcond/gates.py index 83cbaa9..7aed98e 100644 --- a/src/lizzy/_core/bcond/gates.py +++ b/src/lizzy/_core/bcond/gates.py @@ -4,21 +4,17 @@ # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with this program. If not, see . -from dataclasses import dataclass -import numpy as np - - class Inlet(): """A class representing an inlet boundary condition. Parameters ---------- + name : str + Name of the inlet. p_value : float Initial pressure value at the inlet [Pa]. - name : str, optional - Label assigned to the inlet. Will be used to select the inlet in future operations. If none assigned, a default 'unnamed_inlet' name is given. """ - def __init__(self, p_value:float, name:str = "unnamed_inlet"): + def __init__(self, name:str, p_value:float): super().__init__() self.p_value = p_value self._p0 = p_value @@ -68,8 +64,4 @@ def set_open(self, open: bool): -@dataclass() -class SolverBCs: - dirichlet_idx = np.array([], dtype=int) - dirichlet_vals = np.array([], dtype=float) - p0_idx = np.array([], dtype=int) + diff --git a/src/lizzy/_core/bcond/manager.py b/src/lizzy/_core/bcond/manager.py index 1207d37..5983cbd 100644 --- a/src/lizzy/_core/bcond/manager.py +++ b/src/lizzy/_core/bcond/manager.py @@ -1,12 +1,18 @@ +# Copyright 2025-2026 Simone Bancora, Paris Mulye +# +# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with this program. If not, see . + from .gates import Inlet from typing import Literal -class BCManager: +class GatesManager: """Manager for all boundary condition operations. """ def __init__(self): + self._created_inlets : dict[str, Inlet] = {} self._assigned_inlets : dict[str, Inlet] = {} - self._existing_inlets : dict[str, Inlet] = {} @property def assigned_inlets(self) -> dict[str, Inlet]: @@ -18,28 +24,11 @@ def assigned_inlets(self) -> dict[str, Inlet]: def existing_inlets(self) -> dict[str, Inlet]: """Dictionary of inlets that exist in the model (read-only). """ - return self._existing_inlets - - def create_inlet(self, initial_pressure_value:float, name:str = None) -> Inlet: - """Creates a new inlet and add it to the :attr:`~lizzy.bcond.bcond.BCManager.existing_inlets` dictionary. - - Parameters - ---------- - initial_pressure_value : float - Initial pressure value at the inlet. - name : str, optional - Label assigned to the inlet. Will be used to select the inlet in future operations. If none assigned, a default 'Inlet_{N}'name is given, where N is an incremental number of existing inlets. + return self._created_inlets - Returns - ------- - :class:`~lizzy.bcond.bcond.Inlet` - The created inlet object. - """ - if name is None: - inlet_count = len(self._existing_inlets) - name = f"Inlet_{inlet_count}" - new_inlet = Inlet(initial_pressure_value, name) - self._existing_inlets[name] = new_inlet + def create_inlet(self, name:str, initial_pressure_value:float) -> Inlet: + new_inlet = Inlet(name, initial_pressure_value) + self._created_inlets[name] = new_inlet return new_inlet def _fetch_inlet(self, inlet_selector:Inlet | str) -> Inlet: @@ -47,7 +36,7 @@ def _fetch_inlet(self, inlet_selector:Inlet | str) -> Inlet: return inlet_selector else: try: - selected_inlet = self._existing_inlets[inlet_selector] + selected_inlet = self._created_inlets[inlet_selector] except KeyError: raise KeyError(f"Inlet '{inlet_selector}' is not found in existing inlets. Check the name, or create the inlet first using `LizzyModel.create_inlet`.") return selected_inlet diff --git a/src/lizzy/_core/cvmesh/construction.py b/src/lizzy/_core/cvmesh/construction.py index 5d23179..8c46160 100644 --- a/src/lizzy/_core/cvmesh/construction.py +++ b/src/lizzy/_core/cvmesh/construction.py @@ -7,116 +7,213 @@ from __future__ import annotations from typing import TYPE_CHECKING if TYPE_CHECKING: - from lizzy._core.solver import FillSolver + from lizzy._core.solver.fillsolver import FillSolver import numpy as np from .collections import nodes, lines, elements -from lizzy._core.cvmesh.entities import Node, Line, Triangle, CV - - -def CreateNodes(mesh_data): - """ - Creates a Nodes. Returns a "nodes" list. - """ - nodes_coords = np.array(mesh_data['all_nodes_coords']) - all_nodes = nodes([]) - for n, coords in enumerate(nodes_coords): - node = Node(coords[0], coords[1], coords[2]) - node.idx = n - all_nodes.append(node) - all_nodes.XYZ = nodes_coords # redundant but might be useful - all_nodes.N = len(nodes_coords) - return all_nodes - -def CreateLines(mesh_data, triangles): - """ - Creates Lines. Returns a "lines" list. - """ - all_lines = lines([]) - line_counter = 0 - for n, triangle in enumerate(triangles): - #line 1 - line_1 = Line(triangle.nodes[0], triangle.nodes[1]) - line_1.idx = line_counter - line_1.triangles.append(triangle) - line_1.triangle_ids.append(triangle.idx) - line_counter += 1 - all_lines.append(line_1) - # line 2 - line_2 = Line(triangle.nodes[1], triangle.nodes[2]) - line_2.idx = line_counter - line_2.triangles.append(triangle) - line_2.triangle_ids.append(triangle.idx) - line_counter += 1 - all_lines.append(line_2) - # line 3 - line_3 = Line(triangle.nodes[2], triangle.nodes[0]) - line_3.idx = line_counter - line_3.triangles.append(triangle) - line_3.triangle_ids.append(triangle.idx) - line_counter += 1 - all_lines.append(line_3) - - # reference created lines to triangle - triangle.lines.append(line_1) - triangle.line_ids.append(line_1.idx) - triangle.lines.append(line_2) - triangle.line_ids.append(line_2.idx) - triangle.lines.append(line_3) - triangle.line_ids.append(line_3.idx) - - all_lines.N = len(all_lines) - return all_lines - -def CreateTriangles(mesh_data, nodes): - """ - Creates triangles. Returns a "triangles" list. - - Preliminary calculations (pre-processing) for tri elements. - Put triangles in planes and calculate Jacobians and areas: - A_el = A_xi * det(J) = 0.5 * abs(det(J)) - """ - conn = mesh_data['nodes_conn'] - all_triangles = elements([]) - for n, local_conn in enumerate(conn): - node_1 = nodes[local_conn[0]] - node_2 = nodes[local_conn[1]] - node_3 = nodes[local_conn[2]] - tri = Triangle(node_1, node_2, node_3) - tri.idx = n - tri.node_ids = [node_1.idx, node_2.idx, node_3.idx] - all_triangles.append(tri) - # also assign triangle to nodes - node_1.triangles.append(tri) - node_1.triangle_ids.append(tri.idx) - node_2.triangles.append(tri) - node_2.triangle_ids.append(tri.idx) - node_3.triangles.append(tri) - node_3.triangle_ids.append(tri.idx) - - # assign material_tag tag - for key in mesh_data['physical_domains']: - for i in mesh_data['physical_domains'][key]: - all_triangles[i].material_tag = key - - all_triangles.nodes_conn_table = mesh_data['nodes_conn'] # needed? - all_triangles.N = len(all_triangles) - return all_triangles - -def CreateControlVolumes(nodes : list[Node], fill_solver : FillSolver): - # for every nodes: - CVs : list[CV] = [] - for node in nodes: - CVs.append(CV(node)) - CVs = np.array(CVs) - # reference support CVs - for cv in CVs: - connected_nodes = cv.node.node_ids - cv.support_CVs = CVs[connected_nodes] - cv.GetCVLines() - cv.CheckFluxNormalOrientations() - cv.precompute_flux_terms() # this assignes cv.flux_terms, which is an array of variable size (len = n support triangles) - cv.support_triangle_ids = np.array([tri.idx for tri in cv.support_triangles]) #not needed anymore - fill_solver.map_cv_id_to_support_triangle_ids[cv.idx] = np.array([tri.idx for tri in cv.support_triangles]) #TODO should this be in Mesh - fill_solver.map_cv_id_to_flux_terms[cv.idx] = cv.flux_terms - return CVs \ No newline at end of file +from lizzy._core.cvmesh.entities import Node, Line, BoundaryLine, Triangle, CV + +import timeit + +class MeshView: + + def __init__(self): + self.n_nodes:int=0 + self.n_lines:int=0 + self.n_triangles:int=0 + self.node_idx_to_node_idxs: list[np.ndarray] = [] + self.node_idx_to_tri_idxs: list[np.ndarray] = [] + self.node_idx_to_flux_ndarray: list[np.ndarray] = [] + self.phys_boundary_name_to_node_idxs:dict = {} #for dirichlet mostly + self.phys_boundary_name_to_boundary_line_idxs:dict = {} + self.boundary_line_idx_to_node_idxs: np.ndarray = None + + +class MeshBuilder(): + def __init__(self): + self.n_nodes = 0 + self.n_triangles = 0 + self.n_lines = 0 + self.node_idx_to_node_idxs = None + self.node_idx_to_line_idxs = None + self.node_idx_to_tri_idxs_buffer = None + self.node_idx_to_tri_idxs = None + + self.line_idx_to_node_idxs = None + self.line_idx_to_line_idxs = None + self.line_idx_to_triangle_idxs = None + + self.triangle_idx_to_node_idxs = None + self.triangle_idx_to_triangle_idxs = None + self.triangle_idx_to_line_idxs = None + + self.node_idx_to_tri_idxs_for_fill_solver = None + + + def create_cross_referencing_maps(self, n_nodes, n_lines, n_triangles, tri_conn): + capacity_tris_per_node = 8 # initial buffer size + node_idx_to_tri_idxs_buffer = np.full((n_nodes, capacity_tris_per_node), -1, dtype=np.int32) + tri_idxs_local_pointer = np.zeros(n_nodes, dtype=np.uint8) + triangle_idx_to_line_idxs = np.empty((n_triangles, 3), dtype=np.int32) + line_idx_to_node_idxs = np.empty((n_lines, 2), dtype=np.int32) + + line_nodes_from_conn_selectors = [[0,1],[1,2],[2,0]] + for tri_id in range(n_triangles): + local_conn = tri_conn[tri_id] + + # populate `line_idx_to_node_idxs` + local_conn_pairs = [(local_conn[pair[0]], local_conn[pair[1]]) for pair in line_nodes_from_conn_selectors] + for j in range(3): + line_idx = tri_id*3+j + local_node_idxs = local_conn_pairs[j] + line_idx_to_node_idxs[line_idx] = local_node_idxs + + # populate `triangle_idx_to_line_idxs` + triangle_idx_to_line_idxs[tri_id, j] = line_idx + + # populate `node_idx_to_tri_idxs_buffer` + for local_node_id_selector in range(3): + node_id = local_conn[local_node_id_selector] + # check if we still have room for more ids + if tri_idxs_local_pointer[node_id] > capacity_tris_per_node - 1: + # increase bufefr size + capacity_tris_per_node *=2 + node_idx_to_tri_idxs_buffer = np.resize(node_idx_to_tri_idxs_buffer, (n_nodes, capacity_tris_per_node)) + # write the triangle id in the buffer (which is initially 5 tris per node) + node_idx_to_tri_idxs_buffer[node_id, tri_idxs_local_pointer[node_id]] = tri_id + # move local pointer + tri_idxs_local_pointer[node_id] +=1 + + # store + self.node_idx_to_tri_idxs_buffer = node_idx_to_tri_idxs_buffer + self.line_idx_to_node_idxs = line_idx_to_node_idxs + self.triangle_idx_to_node_idxs = tri_conn + self.triangle_idx_to_line_idxs = triangle_idx_to_line_idxs + + def create_entities(self, n_nodes, n_triangles, n_lines, node_coords, tri_conn, physical_lines_conn): + # preallocate lists + new_nodes = nodes([None]*n_nodes) + new_lines = lines([None]*n_lines) + new_boundary_lines = lines([None]*len(physical_lines_conn)) + new_triangles = elements([None]*n_triangles) + # create nodes + for i in range(n_nodes): + new_nodes[i] = Node(node_coords[i,0], node_coords[i,1], node_coords[i,2], i) + new_nodes.XYZ = node_coords + new_nodes.N = len(new_nodes) + # create lines + for i in range(n_lines): + local_conn = self.line_idx_to_node_idxs[i] + local_node_objs = [new_nodes[idx] for idx in local_conn] + new_lines[i] = Line(*local_node_objs, i) + new_lines.N = len(new_lines) + # create boundary lines: + for i in range(len(physical_lines_conn)): + local_conn = physical_lines_conn[i] + local_node_objs = [new_nodes[idx] for idx in local_conn] + new_boundary_lines[i] = BoundaryLine(*local_node_objs, i) + # create triangles + for i in range(n_triangles): + local_nodes_conn = self.triangle_idx_to_node_idxs[i] + local_node_objs = [new_nodes[idx] for idx in local_nodes_conn] + local_lines_conn = self.triangle_idx_to_line_idxs[i] + local_line_objs = [new_lines[idx] for idx in local_lines_conn] + + new_triangles[i] = Triangle(*local_node_objs, *local_line_objs, i) + new_triangles.nodes_conn_table = tri_conn + new_triangles.N = len(new_triangles) + + + + return new_nodes, new_lines, new_triangles, new_boundary_lines + + def assign_material_tags_to_elements(self, mesh_data, triangles:list[Triangle]): + # assign material_tag tag. key is a string (name of physical group) + for key in mesh_data['physical_domains']: + for i in mesh_data['physical_domains'][key]: + triangles[i].material_tag = key + + + + def assign_varying_number_references(self, nodes:list[Node], triangles, tri_conn): + node_idx_to_node_idxs = [None]*len(nodes) + node_idx_to_tri_idxs = [None]*len(nodes) + + for i in range(len(nodes)): + # assign triangles to nodes (varying number) + tri_ids_buffer = self.node_idx_to_tri_idxs_buffer[i] + tri_ids = tri_ids_buffer[tri_ids_buffer >= 0] + node_idx_to_tri_idxs[i] = np.array(tri_ids) + nodes[i].triangle_ids = tri_ids + triangle_objs = [triangles[idx] for idx in tri_ids] + nodes[i].triangles = triangle_objs + + connected_node_idxs = [idx for triangle in triangle_objs for idx in triangle.node_ids] + connected_node_idxs = list(set(connected_node_idxs)-{i}) + nodes[i].node_ids = connected_node_idxs + nodes[i].nodes = [nodes[idx] for idx in connected_node_idxs] + node_idx_to_node_idxs[i] = np.array(connected_node_idxs) + self.node_idx_to_node_idxs = node_idx_to_node_idxs + self.node_idx_to_tri_idxs = node_idx_to_tri_idxs + return node_idx_to_node_idxs, node_idx_to_tri_idxs + + + def build_mesh(self, mesh_data): + print("Creating Mesh...") + mesh_view = MeshView() + tri_conn:np.ndarray = mesh_data['nodes_conn'] + node_coords:np.ndarray = mesh_data['all_nodes_coords'] + n_nodes = node_coords.shape[0] + n_triangles = tri_conn.shape[0] + n_lines = n_triangles*3 + mesh_view.n_nodes = n_nodes + mesh_view.n_lines = n_lines + mesh_view.n_triangles = n_triangles + # time_create_cross_referencing(self, tri_conn) + self.create_cross_referencing_maps(n_nodes, n_lines, n_triangles, tri_conn) + # time_create_entities(self, node_coords, tri_conn) + physical_lines_conn = mesh_data["physical_lines_conn"] + phys_boundary_name_to_boundary_line_idxs = mesh_data["physical_lines"] + new_nodes, new_lines, new_triangles, new_boundary_lines = self.create_entities(n_nodes, n_triangles, n_lines, node_coords, tri_conn, physical_lines_conn) + # time_assign_varying_number_references(self, new_nodes, new_triangles, tri_conn) + node_idx_to_node_idxs, node_idx_to_tri_idxs = self.assign_varying_number_references(new_nodes, new_triangles, tri_conn) + mesh_view.node_idx_to_node_idxs = node_idx_to_node_idxs + mesh_view.node_idx_to_tri_idxs = node_idx_to_tri_idxs + mesh_view.phys_boundary_name_to_node_idxs = mesh_data['physical_nodes'] + mesh_view.phys_boundary_name_to_boundary_line_idxs = phys_boundary_name_to_boundary_line_idxs + mesh_view.boundary_line_idx_to_node_idxs = physical_lines_conn + cvs = self.create_control_volumes(new_nodes) + + self.assign_material_tags_to_elements(mesh_data, new_triangles) + return new_nodes, new_lines, new_boundary_lines, new_triangles, cvs, mesh_view + + + def create_control_volumes(self, nodes : list[Node]): + # for every nodes: + n_nodes = len(nodes) + CVs : list[CV] = [None]*n_nodes + for i in range(n_nodes): + CVs[i] = CV(nodes[i]) + # node_idx_to_flux_ndarray[i] = CVs[i].compute_flux_terms() + # reference support CVs + for cv in CVs: + connected_nodes = cv.node.node_ids + cv.support_CVs = [CVs[i] for i in connected_nodes] + return np.array(CVs) + + +def time_create_cross_referencing(self, tri_conn): + elapsed = timeit.timeit(lambda: self.create_cross_referencing_maps(tri_conn), number=1000) + print(f"CREATE CROSS REFERENCING: Average per run: {elapsed/1000:.10f} seconds") + +def time_create_entities(self, node_coords, tri_conn): + elapsed = timeit.timeit(lambda: self.create_entities(node_coords, tri_conn), number=1000) + print(f"CREATE ENTITIES: Average per run: {elapsed/1000:.10f} seconds") + +def time_assign_varying_number_references(self, new_nodes, new_triangles, tri_conn): + elapsed = timeit.timeit(lambda: self.assign_varying_number_references(new_nodes, new_triangles, tri_conn), number=1000) + print(f"ASSIGN VARYING NUMBER REFERENCES: Average per run: {elapsed/1000:.10f} seconds") + + + + diff --git a/src/lizzy/_core/cvmesh/entities.py b/src/lizzy/_core/cvmesh/entities.py index dd1a487..dbf1138 100644 --- a/src/lizzy/_core/cvmesh/entities.py +++ b/src/lizzy/_core/cvmesh/entities.py @@ -4,14 +4,33 @@ # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with this program. If not, see . + + +from __future__ import annotations +from typing import TYPE_CHECKING +if TYPE_CHECKING: + from lizzy._core.cvmesh.entities import Node, Line, Triangle, CV import numpy as np -from dataclasses import dataclass, field + class Node: """Class representing a mesh node.""" - def __init__(self, x: float, y: float, z: float): + + __slots__ = ( + "coords", + "idx", + "p", + "triangles", + "triangle_ids", + "lines", + "line_ids", + "nodes", + "node_ids" + ) + + def __init__(self, x: float, y: float, z: float, idx:int): self.coords = np.array([x, y, z]) - self.idx : int = 0 + self.idx : int = idx self.p : float = 0 self.triangles : list[Triangle] = [] self.triangle_ids : list[int] = [] @@ -21,9 +40,32 @@ def __init__(self, x: float, y: float, z: float): self.node_ids : list[int] = [] def __str__(self): return "Node ID: " + str(self.idx) + + def assign_triangle(self, triangle:Triangle): + self.triangles.append(triangle) + self.triangle_ids.append(triangle.idx) + class Element2D: + + __slots__ = ( + "idx", + "material_tag", + "A", + "h", + "k", + "grad_N", + "porosity", + "nodes", + "node_ids", + "nodes_coords", + "lines", + "line_ids", + "centroid", + "v" + ) + def __init__(self): self.idx : int = 0 self.material_tag : str = "" @@ -35,12 +77,11 @@ def __init__(self): self.nodes : tuple[Node] = () self.node_ids : list[int] = [] self.nodes_coords : np.ndarray = None - self.lines : list[Line] = [] + self.lines : tuple[Line] = [] self.line_ids : list[int] = [] self.centroid = np.zeros(3) self.v : np.ndarray = np.zeros((3,1)) - class Triangle(Element2D): """Class representing a triangular element. """ @@ -49,8 +90,13 @@ class Triangle(Element2D): dNdxi = np.array([[-1, -1], [1, 0], [0, 1]]) - def __init__(self, node_1:Node, node_2:Node, node_3:Node): + def __init__(self, node_1:Node, node_2:Node, node_3:Node, line_1:Line, line_2:Line, line_3:Line, n:int): super().__init__() + self.idx = n + self.nodes = (node_1, node_2, node_3) + self.node_ids = (node_1.idx, node_2.idx, node_3.idx) + self.lines = (line_1, line_2, line_3) + self.line_ids = (line_1.idx, line_2.idx, line_3.idx) x = np.array((node_1.coords, node_2.coords, node_3.coords)) J = np.array([ [x[1,0]-x[0,0], x[2,0]-x[0,0]], @@ -65,7 +111,6 @@ def __init__(self, node_1:Node, node_2:Node, node_3:Node): v = x[0,:] - x[2,:] n = np.cross(u, v) self.n = n / np.linalg.norm(n) - self.nodes = (node_1, node_2, node_3) self.grad_N = (Triangle.dNdxi @ dxidX).T self.A = 0.5 * detJ self.centroid = x.mean(0) @@ -139,29 +184,31 @@ def __str__(self): class Line: """Class representing a line between two nodes in the mesh. """ - def __init__(self, node_1:Node, node_2:Node): + + __slots__ = ( + "idx", + "nodes", + "triangles", + "triangle_ids", + "midpoint" + ) + def __init__(self, node_1:Node, node_2:Node, n): self.nodes = (node_1, node_2) - self.idx : int = 0 - self.midpoint : np.ndarray = self.ComputeMidPoint() - self.n : np.ndarray = self.ComputeNormal() + self.idx : int = n self.triangles = [] self.triangle_ids = [] + self.midpoint : np.ndarray = self._compute_midpoint() - def ComputeMidPoint(self): + def _compute_midpoint(self): x1 = self.nodes[0].coords x2 = self.nodes[1].coords return np.array((x1, x2)).mean(0) - #TODO: this is wrong in 3D but the line normal is not being used at the moment - def ComputeNormal(self): - x1 = self.nodes[0].coords - x2 = self.nodes[1].coords - DX = x1 - x2 - l = np.linalg.norm(DX) - nx = DX[1]/l - ny = -DX[0]/l - nz = 0 - return np.array((nx, ny, nz)) +class BoundaryLine(Line): + __slots__ = ("length") + def __init__(self, node_1:Node, node_2:Node, n): + super().__init__(node_1, node_2, n) + self.length = np.linalg.norm(self.nodes[0].coords - self.nodes[1].coords) class CV: @@ -176,12 +223,12 @@ def __init__(self, node:Node): self.support_CVs : list[CV] = [] self.support_lines : list[Line] = [] self.support_nodes : list[Node] = [] + self.support_node_ids = None self.support_triangles : list[Triangle] = node.triangles self.support_triangle_ids = None self.edges = [] - self.flux_terms = [] - self.cv_lines = [] - self.A, self.vol = self._calculate_area_and_volume() + self.cv_lines = self._create_cv_lines() + self._check_flux_normals() # The CV has this structure: @@ -189,8 +236,8 @@ def __init__(self, node:Node): # cv_lines (of each support triangle) = [ [line1, line2], [line1, lined], [line1, line2], ... ] # each line has normal, length - def GetCVLines(self): - self.cv_lines = [] + def _create_cv_lines(self): + cv_lines = [] for tri in self.support_triangles: elem_side_lines = [] for line in tri.lines: @@ -205,41 +252,21 @@ def GetCVLines(self): x2 = elem_side_lines[1].midpoint centroid = tri.centroid cv_lines_tri = [CVLine(x1, centroid, tri.n), CVLine(centroid, x2, tri.n)] - self.cv_lines.append(cv_lines_tri) - - def precompute_flux_terms(self): - for i, tri in enumerate(self.support_triangles): + cv_lines.append(cv_lines_tri) + return cv_lines + + def compute_flux_terms(self): # this computes flux_terms, which is an array of variable size (len = n support triangles) + n_support_triangles = len(self.support_triangles) + flux_terms = np.empty((n_support_triangles, 3)) + for i in range(n_support_triangles): + tri = self.support_triangles[i] line1 = self.cv_lines[i][0] # PSEUDO CODE ALL TO CHECK AND RE-WRITE line2 = self.cv_lines[i][1] n1 = line1.n n2 = line2.n flux_term = (-n1 * line1.l + -n2 * line2.l) * tri.h - self.flux_terms.append(flux_term) - self.flux_terms = np.array(self.flux_terms) - - - - # @staticmethod - # def polygon_area(points): - # """ - # Calculate the area of an irregular polygon using the Shoelace formula. - - # points: A list of (x, y) tuples representing the vertices of the polygon. - # The vertices should be provided in order (clockwise or counterclockwise). - # return: The area of the polygon. - # """ - # n = len(points) - # if n < 3: - # raise ValueError("A polygon must have at least 3 points.") - - # # Compute the Shoelace formula - # area = 0 - # for i in range(n): - # x1, y1 = points[i] - # x2, y2 = points[(i + 1) % n] # Next vertex (wrap around) - # area += x1 * y2 - y1 * x2 - - # return abs(area) / 2 + flux_terms[i] = flux_term + return flux_terms def _polygon_area_3d(self, points): """ @@ -272,10 +299,8 @@ def _polygon_area_3d(self, points): area = 0.5 * abs(np.dot(normal, cross_sum)) return area - - # TODO: check if this is correct : - def _calculate_area_and_volume(self) -> tuple[float, float]: + def calculate_area_and_volume(self) -> tuple[float, float]: area = 0 vol = 0 for tri in self.support_triangles: @@ -293,9 +318,11 @@ def _calculate_area_and_volume(self) -> tuple[float, float]: slice_vol = slice_area*tri.h*tri.porosity area += slice_area vol += slice_vol - return area, vol + self.A = area + self.vol = vol + - def CheckFluxNormalOrientations(self): + def _check_flux_normals(self): # by convention, normals are oriented outwards from the CV. This function checks and enforces that for i, tri in enumerate(self.support_triangles): for line in self.cv_lines[i]: @@ -307,20 +334,20 @@ def CheckFluxNormalOrientations(self): if dist_outer < dist_inner: line.n = -line.n - class CVLine: def __init__(self, p1, p2, tri_normal): self.p1 = p1 self.p2 = p2 - self.l = 0 - self.n = None self.tri_normal = tri_normal - self.ComputeLengthAndNormal() + self.l = self._compute_length() + self.n = self._compute_normal() - def ComputeLengthAndNormal(self): - DX = self.p1 - self.p2 + def _compute_length(self): + d = self.p1 - self.p2 self.midpoint = 0.5*(self.p1 + self.p2) - self.l = np.linalg.norm(DX) - n = np.cross(self.tri_normal, DX) - self.n = n / np.linalg.norm(n) - + return np.linalg.norm(d) + + def _compute_normal(self): + d = self.p1 - self.p2 + n = np.cross(self.tri_normal, d) + return n / np.linalg.norm(n) \ No newline at end of file diff --git a/src/lizzy/_core/cvmesh/mesh.py b/src/lizzy/_core/cvmesh/mesh.py index 28a80dd..1833b8d 100644 --- a/src/lizzy/_core/cvmesh/mesh.py +++ b/src/lizzy/_core/cvmesh/mesh.py @@ -7,15 +7,15 @@ from __future__ import annotations from typing import TYPE_CHECKING if TYPE_CHECKING: - from lizzy._core.solver import FillSolver from lizzy._core.io import Reader from lizzy._core.materials import MaterialManager, Rosette, PorousMaterial - from lizzy._core.cvmesh.entities import Node, Line, Triangle, CV + from lizzy._core.cvmesh.entities import Node, Line, BoundaryLine, Triangle, CV import numpy as np -from .construction import CreateNodes, CreateLines, CreateTriangles, CreateControlVolumes +from .construction import MeshBuilder, MeshView from .collections import nodes, lines, elements + class Mesh: r""" A class representing a FE/CV mesh. @@ -29,109 +29,18 @@ class Mesh: """ def __init__(self, mesh_reader:Reader): + self.mesh_view:MeshView = None self.mesh_data = mesh_reader.mesh_data self.nodes : list[Node] = nodes([]) + self.lines : list[Line] = lines([]) + self.boundary_lines : list[BoundaryLine] = lines([]) self.triangles : list[Triangle] = elements([]) self.tetras = elements([]) - self.lines : list[Line] = lines([]) self.CVs :list[CV] = [] - self.boundaries = mesh_reader.mesh_data['physical_nodes'] - self.preprocessed = False # Init methods: - self.PopulateFromMeshData(self.mesh_data) - self.CrossReferenceEntities() - - ################## ALL THIS BLOCK TO BE REMOVED - # create mesh of CVs for visualisation only - cv_mesh_nodes = [] - cv_mesh_conn = [] - nodes_counter = 0 - for cv in self.CVs: - for two_lines_tri in cv.cv_lines: - for line in two_lines_tri: - line_conn = [] - p1 = line.p1 - p2 = line.p2 - cv_mesh_nodes.append(p1) - line_conn.append(nodes_counter) - nodes_counter += 1 - cv_mesh_nodes.append(p2) - line_conn.append(nodes_counter) - nodes_counter += 1 - cv_mesh_conn.append(line_conn) - - self.cv_mesh_nodes = cv_mesh_nodes - self.cv_mesh_conn = cv_mesh_conn - - ################## ALL THIS BLOCK TO BE REMOVED - - def PopulateFromMeshData(self, mesh_data): - """ - Takes mesh data dictionary and initialises all mesh attributes: nodes, lines, triangles. - - Parameters - ---------- - mesh_data : dict - Dictionary of all mesh data read from mesh file - """ - self.nodes = CreateNodes(mesh_data) - self.triangles = CreateTriangles(mesh_data, self.nodes) - self.lines = CreateLines(mesh_data, self.triangles) - - def preprocess(self, material_manager: MaterialManager, fill_solver: FillSolver): - """ Pre-processes the mesh before simulation. Assigns material properties to elements, creates control volumes (CVs), and prepares data structures for simulation.""" - materials = material_manager.assigned_materials - rosettes = material_manager.assigned_rosettes - for tri in self.triangles: - try: - material : PorousMaterial = materials[tri.material_tag] - if material.is_isotropic: - tri.k = material.k_diag - else: - rosette : Rosette = rosettes[tri.material_tag] - u, v, w = rosette.project_along_normal(tri.n) - R = np.array([u, v, w]).T - tri.k = R @ material.k_diag @ R.T - tri.porosity = materials[tri.material_tag].porosity - tri.h = materials[tri.material_tag].thickness - except KeyError: - exit(f"Mesh contains unassigned material tag: {tri.material_tag}") - self.CVs = CreateControlVolumes(self.nodes, fill_solver) - # create a hashmap for CV id: [ids of supporting elements] - print("Mesh pre-processing completed\n") - - - - #TODO: for numba: - self.triangle_id_lists = [np.array(n.triangle_ids, dtype=np.int32) for n in self.nodes] - - - - - self.preprocessed = True - - def CrossReferenceEntities(self): - """ - Creates hierarchical connections between all objects that constitute the mesh: Nodes, Lines, Elements. - The purpose is to make any given object accessible from any other given object, to which it would be linked as an attribute. - Once this method is called on a mesh, all references should be created. - - Example - -------- - Given a mesh which has been cross referenced, fetch the nodes of the fourth element: - - >>> nodes = mesh.elements[3].nodes - """ - - # go through nodes and cross reference connected nodes, to help fetch support CVs - for node in self.nodes: - connected_nodes_ids = [] - for tri in node.triangles: - connected_nodes_ids.append(tri.node_ids) - connected_nodes_ids = np.unique(connected_nodes_ids).tolist() - connected_nodes_ids.remove(node.idx) - node.node_ids = connected_nodes_ids + self.mb = MeshBuilder() + self.nodes, self.lines, self.boundary_lines, self.triangles, self.CVs, self.mesh_view = self.mb.build_mesh(self.mesh_data) def empty_cvs(self): for cv in self.CVs: diff --git a/src/lizzy/_core/datatypes/simparams.py b/src/lizzy/_core/datatypes/simparams.py index 5af5947..e1ae7e0 100644 --- a/src/lizzy/_core/datatypes/simparams.py +++ b/src/lizzy/_core/datatypes/simparams.py @@ -12,8 +12,6 @@ class SimulationParameters: Attributes ---------- - mu : float - Viscosity [Pa s] wo_delta_time : float Interval of simulation time between solution write-outs [s]. Default: -1 (write-out every numerical time step) fill_tolerance : float @@ -23,7 +21,6 @@ class SimulationParameters: """ - mu: float = 0.1 wo_delta_time: float = -1 fill_tolerance: float = 0.01 has_been_assigned : bool = False @@ -36,7 +33,6 @@ def print_current(self): """Prints the currently assigned simulation parameters to the console.""" params = textwrap.dedent(rf""" Currently assigned simulation parameters: - - "mu": {self.mu} [Pa s], - "wo_delta_time": {self.wo_delta_time} [s], - "fill_tolerance": {self.fill_tolerance}, - "end_step_when_sensor_triggered": {self.end_step_when_sensor_triggered} @@ -54,7 +50,6 @@ def assign(self, **kwargs): Keyword arguments corresponding to parameter names and their new values. Each key must be a valid attribute of the `SimulationParameters` class, otherwise, an `AttributeError` is raised. Valid parameters are: - - ``mu``: viscosity [Pa s] - ``wo_delta_time``: interval of simulation time between solution write-outs [s]. Default: -1 (write-out every numerical time step) - ``fill_tolerance``: tolerance on the fill factor to consider a CV as filled. Default: 0.01 - ``end_step_when_sensor_triggered``: if True, ends current solution step and creates a write-out when a sensor changes state. Default: False diff --git a/src/lizzy/_core/datatypes/solution.py b/src/lizzy/_core/datatypes/solution.py index 3fb49d1..c28cb6a 100644 --- a/src/lizzy/_core/datatypes/solution.py +++ b/src/lizzy/_core/datatypes/solution.py @@ -9,24 +9,24 @@ class Solution: Attributes ---------- - time_steps_in_solution : int - The number of time steps stored in the solution. - time_step_idx : ndarray of int, shape (time_steps_in_solution,) - The indices of the time steps stored in the solution. The last index corresponds to the time step number at which this solution was saved. - p : np.ndarray of float, shape (time_steps_in_solution, N_nodes) + n_time_states : int + The number of time states stored in the solution. + time_step_idx : ndarray of int, shape (n_time_states,) + The indices of the time steps that were stored as time states in the solution. The last index corresponds to the time step number at which this solution was saved. + p : np.ndarray of float, shape (n_time_states, N_nodes) The pressure values at each step. - v : np.ndarray of float, shape (time_steps_in_solution, N_elements, 3) + v : np.ndarray of float, shape (n_time_states, N_elements, 3) The velocity values at each step. - v_nodal : np.ndarray of float, shape (time_steps_in_solution, N_nodes, 3) + v_nodal : np.ndarray of float, shape (n_time_states, N_nodes, 3) The nodal velocity values at each step. - time : np.ndarray of float, shape (time_steps_in_solution,) + time : np.ndarray of float, shape (n_time_states,) The simulation time values at each step. - fill_factor : np.ndarray of float, shape (time_steps_in_solution, N_nodes) + fill_factor : np.ndarray of float, shape (n_time_states, N_nodes) The fill factor values at each step. - free_surface : np.ndarray of int, shape (time_steps_in_solution, N_nodes) + free_surface : np.ndarray of int, shape (n_time_states, N_nodes) The free surface values at each step. """ - time_steps_in_solution : int + n_time_states : int time_step_idx : np.ndarray p : np.ndarray v : np.ndarray diff --git a/src/lizzy/_core/io/reader.py b/src/lizzy/_core/io/reader.py index f756701..53ee0a3 100644 --- a/src/lizzy/_core/io/reader.py +++ b/src/lizzy/_core/io/reader.py @@ -4,8 +4,7 @@ # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with this program. If not, see . -import os -import shutil + from pathlib import Path from enum import Enum, auto import numpy as np @@ -14,36 +13,6 @@ import numpy as np -# Extract lines: they will be REPEATED -def extract_lines(nodes_conn): - """ - Extract lines connectivity from a nodes connectivity table. - - Parameters - ---------- - nodes_conn: list - The nodes connectivity table of the mesh - - Returns - ------- - lines_conn : list - The lines connectivity table - """ - lines_conn = [] - n_geom = len(nodes_conn[0]) - for e in nodes_conn: - candidate_lines_conns = [] - for i in range(n_geom-1): - candidate_lines_conns.append([e[i], e[i+1]]) - candidate_lines_conns.append([e[i+1], e[0]]) - for line_conn in candidate_lines_conns: - line_test_1 = [line_conn[0], line_conn[1]] - line_test_2 = [line_conn[1], line_conn[0]] - # if normals give problems, might be necessary to not get repeated lines - if line_test_1 not in lines_conn and line_test_2 not in lines_conn: - lines_conn.append(line_conn) - return lines_conn - def extract_unique_nodes(node_ids_list): """ @@ -68,21 +37,21 @@ class Reader: The path to the mesh file. """ - def __init__(self, mesh_path:str): + def __init__(self, ): self.mesh_data:dict = {} # A dict containing all the mesh info from the gmsh file - self.mesh_path = Path(mesh_path) - self.case_name = self.__read_case_name() - self.__read_mesh_file() + self.case_name:str = None - def __read_mesh_file(self): - print(f"Reading mesh file: {self.mesh_path}") + def read_mesh_file(self, mesh_path:str): + mesh_path = Path(mesh_path) + self.case_name = self.__read_case_name(mesh_path) + print(f"Reading mesh file: {mesh_path}") _format = self._detect_format() match _format: case Format.MSH: - self.mesh_data = self._read_gmsh_file() + self.mesh_data = self._read_gmsh_file(mesh_path) - def __read_case_name(self): - case_name = self.mesh_path.stem + def __read_case_name(self, mesh_path:Path): + case_name = mesh_path.stem return case_name def _detect_format(self): @@ -90,14 +59,16 @@ def _detect_format(self): NOT IMPLEMENTED""" return Format.MSH - def _read_gmsh_file(self) -> dict: + def _read_gmsh_file(self, mesh_path:Path) -> dict: """ Reads a mesh file in .msh format (ASCII 4). Initialises all mesh attributes. """ try: - mesh_file = meshio.read(self.mesh_path, file_format="gmsh") + mesh_file = meshio.read(mesh_path, file_format="gmsh") except meshio._exceptions.ReadError: - raise FileNotFoundError(f"Mesh file not found: {self.mesh_path}") + raise FileNotFoundError(f"Mesh file not found: {mesh_path}") + + # TODO: this block until return is very slow all_nodes_coords : np.ndarray = mesh_file.points physical_domain_names = [] physical_line_names = [] @@ -120,12 +91,10 @@ def _read_gmsh_file(self) -> dict: # get node ids for nodes in the physical lines for key in physical_lines: physical_nodes_ids[key] = extract_unique_nodes(mesh_file.cells_dict["line"][physical_lines[key]]) - lines_conn = extract_lines(nodes_conn) mesh_data = { 'all_nodes_coords' : all_nodes_coords, 'nodes_conn' : nodes_conn, - 'lines_conn' : lines_conn, 'physical_lines_conn' : physical_lines_conn, 'physical_domains' : physical_domains, 'physical_lines' : physical_lines, diff --git a/src/lizzy/_core/io/writer.py b/src/lizzy/_core/io/writer.py index a2e000f..9f1fd77 100644 --- a/src/lizzy/_core/io/writer.py +++ b/src/lizzy/_core/io/writer.py @@ -15,23 +15,12 @@ class Writer: - """ - Handles writing results to output files. - - Attributes - ---------- - mesh : lizzy.Mesh - The Mesh object used in the simulation. - """ - def __init__(self, mesh): - """Class constructor - - Parameters - ---------- - mesh : lizzy.Mesh - The Mesh object of the simulation - """ - self.mesh = mesh + """Handles writing results to output files.""" + def __init__(self): + self._mesh = None + + def assign_mesh(self, mesh): + self._mesh = mesh def save_results(self, solution:Solution, result_name:str, **kwargs): """Save the results contained in the solution dictionary into an XDMF file. @@ -49,16 +38,16 @@ def save_results(self, solution:Solution, result_name:str, **kwargs): if os.path.isdir(destination_path): shutil.rmtree(destination_path) os.makedirs(destination_path, exist_ok=True) - points = self.mesh.nodes.XYZ # Node coordinates, assumed to be (N, 3) - cells = self.mesh.triangles.nodes_conn_table # Triangle connectivity (M, 3) + points = self._mesh.nodes.XYZ # Node coordinates, assumed to be (N, 3) + cells = self._mesh.triangles.nodes_conn_table # Triangle connectivity (M, 3) cells_list = [] for i in range(len(cells)) : cells_list.append(cells[i]) if save_cv_mesh: mesh_cv = meshio.Mesh( - points=self.mesh.cv_mesh_nodes, - cells=[("line", self.mesh.cv_mesh_conn)], # Triangle connectivity + points=self._mesh.cv_mesh_nodes, + cells=[("line", self._mesh.cv_mesh_conn)], # Triangle connectivity ) mesh_cv.write(destination_path / f"{result_name}_CV.vtk") @@ -66,7 +55,7 @@ def save_results(self, solution:Solution, result_name:str, **kwargs): filename = f"{result_name}.xdmf" with meshio.xdmf.TimeSeriesWriter(filename) as writer: writer.write_points_cells(points, [("triangle", cells_list)]) - for i in range(solution.time_steps_in_solution): + for i in range(solution.n_time_states): time = solution.time[i] point_data = { "Pressure" : solution.p[i], "FillFactor" : solution.fill_factor[i], diff --git a/src/lizzy/_core/materials/__init__.py b/src/lizzy/_core/materials/__init__.py index 990ccfd..943e8a7 100644 --- a/src/lizzy/_core/materials/__init__.py +++ b/src/lizzy/_core/materials/__init__.py @@ -1,3 +1,3 @@ from .manager import MaterialManager from .rosette import Rosette -from .materials import PorousMaterial \ No newline at end of file +from .materials import PorousMaterial, Resin \ No newline at end of file diff --git a/src/lizzy/_core/materials/manager.py b/src/lizzy/_core/materials/manager.py index 31008e6..2268992 100644 --- a/src/lizzy/_core/materials/manager.py +++ b/src/lizzy/_core/materials/manager.py @@ -4,8 +4,9 @@ # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. # You should have received a copy of the GNU General Public License along with this program. If not, see . +import sys import numpy as np -from .materials import PorousMaterial +from .materials import PorousMaterial, Resin from .rosette import Rosette class MaterialManager: @@ -15,6 +16,8 @@ def __init__(self): self._existing_materials : dict[str, PorousMaterial] = {} self._assigned_materials : dict[str, PorousMaterial] = {} self._assigned_rosettes : dict[str, Rosette] = {} + self._created_resins: dict[str, Resin] = {} + self._assigned_resin: Resin = None @property def assigned_materials(self) -> dict[str, PorousMaterial]: @@ -22,6 +25,12 @@ def assigned_materials(self) -> dict[str, PorousMaterial]: """ return self._assigned_materials + @property + def assigned_resin(self) -> Resin: + """Resin that has been assigned to the mdoel (read-only). + """ + return self._assigned_resin + @property def assigned_rosettes(self) -> dict[str, Rosette]: """Dictionary of rosettes assigned to materials in mesh regions (read-only). @@ -33,110 +42,52 @@ def existing_materials(self) -> dict[str, PorousMaterial]: """Dictionary of materials that exist in the model, but may have not been assigned yet (read-only). """ return self._existing_materials + + def _check_name_uniqueness_in_dict(self, name, dictionary): + if name in dictionary.keys(): + print(f"ERROR: the name '{name}' has been used more than once. Use unique names.") + sys.exit(1) def fetch_material(self, material_selector:str): - """Fetch an existing material by its label. - Parameters - ---------- - material_selector : str - Label of the material to fetch. - Returns - ------- - :class:`PorousMaterial` - Instance of the selected material. - """ try: selected_material = self._existing_materials[material_selector] except KeyError: - raise KeyError(f"Inlet '{material_selector}' is not found in existing inlets. Check the name, or create the inlet first using `LizzyModel.create_inlet`.") + raise KeyError(f"Inlet '{material_selector}' is not found in existing materials. Check the name, or create the material first using `LizzyModel.create_material`.") return selected_material - - def create_material(self, k1: float, k2: float, k3: float, porosity: float, thickness: float, name: str = None): - """Create a new material that can then be selected and used in the model. - - Parameters - ---------- - k1 : float - Permeability in the first principal direction. - k2 : float - Permeability in the second principal direction. - k3 : float - Permeability in the third principal direction. - porosity : float - Volumetric porosity of the material (porosity = 1 - fibre volume fraction). - thickness : float - Thickness of the material [mm]. - name : str, optional - Label assigned to the material. Necessary to select the material during assignment. If none assigned, a default 'Material_{N}'name is given, where N is an incremental number of existing materials. - - Returns - ------- - :class:`PorousMaterial` - Instance of the created material. - """ - if name is None: - material_count = len(self._existing_materials) - name = f"Material_{material_count}" - new_material = PorousMaterial(k1, k2, k3, porosity, thickness, name) - self._existing_materials[name] = new_material - return new_material - - def create_rosette(self, p1: tuple[float, float, float] = (1.0, 0, 0), p0: tuple[float, float, float] = (0.0, 0.0, 0.0), name: str = None): - """Create a new rosette that can then be selected and used in the model. - - Parameters - ---------- - p1 : tuple[float, float, float] - The first point defining the first axis of the rosette (k1 direction). - p0 : tuple[float, float, float] - The second point defining the first axis of the rosette (k1 direction). Default is (0,0,0). - name : str, optional - Label assigned to the rosette. Necessary to select the rosette during assignment. If none assigned, a default 'Rosette_{N}'name is given, where N is an incremental number of existing rosettes. - - Returns - ------- - :class:`Rosette` - Instance of the created rosette. - """ - if name is None: - rosette_count = len(self._assigned_rosettes) - name = f"Rosette_{rosette_count}" - new_rosette = Rosette(p1, p0, name) - self._assigned_rosettes[name] = new_rosette - return new_rosette + def _fetch_resin(self, resin_selector:str): + try: + selected_resin = self._created_resins[resin_selector] + except KeyError: + raise KeyError(f"Resin '{resin_selector}' was not found. Check the name, or create the resin first using `LizzyModel.create_resin`.") + return selected_resin def _fetch_rosette(self, rosette_selector: str): - """Fetch an existing rosette by its label. - - Parameters - ---------- - rosette_selector : str - Label of the rosette to fetch. - - Returns - ------- - :class:`Rosette` - Instance of the selected rosette. - """ try: selected_rosette = self._assigned_rosettes[rosette_selector] except KeyError: raise KeyError(f"Rosette '{rosette_selector}' is not found in assigned rosettes. Check the name, or create the rosette first using `LizzyModel.create_rosette`.") return selected_rosette + + def create_material(self, name:str, k_vals : tuple[float, float, float], porosity: float, thickness: float): + self._check_name_uniqueness_in_dict(name, self._existing_materials) + new_material = PorousMaterial(name, k_vals, porosity, thickness) + self._existing_materials[name] = new_material + return new_material + + def create_resin(self, name:str, viscosity:float): + self._check_name_uniqueness_in_dict(name, self._created_resins) + new_resin = Resin(name, viscosity) + self._created_resins[name] = new_resin + return new_resin + + def create_rosette(self, name:str, u: tuple[float, float, float] = (1.0, 0, 0)): + self._check_name_uniqueness_in_dict(name, self._assigned_rosettes) + new_rosette = Rosette(name, u) + self._assigned_rosettes[name] = new_rosette + return new_rosette def assign_material(self, material_selector:str, mesh_tag:str, rosette_selector:str | Rosette = None): - """Assign an existing material to a labeled mesh region. - - Parameters - ---------- - material_selector : str - Label of the material to assign. Must correspond to an existing material created with `LizzyModel.create_material`. - mesh_tag : str - Label of the mesh region where to assign the material. - rosette : Rosette, optional - Orientation rosette to apply to the material. If none provided, a default rosette with k1 aligned with the global X axis is assigned. - """ selected_material : PorousMaterial = self.fetch_material(material_selector) if rosette_selector is None: rosette = Rosette((1, 0, 0)) @@ -146,4 +97,8 @@ def assign_material(self, material_selector:str, mesh_tag:str, rosette_selector: rosette = rosette_selector selected_material.assigned = True self._assigned_materials[mesh_tag] = selected_material - self._assigned_rosettes[mesh_tag] = rosette \ No newline at end of file + self._assigned_rosettes[mesh_tag] = rosette + + def assign_resin(self, resin_selector:str): + selected_resin : Resin = self._fetch_resin(resin_selector) + self._assigned_resin = selected_resin diff --git a/src/lizzy/_core/materials/materials.py b/src/lizzy/_core/materials/materials.py index 01a728e..d870e5a 100644 --- a/src/lizzy/_core/materials/materials.py +++ b/src/lizzy/_core/materials/materials.py @@ -12,24 +12,35 @@ class PorousMaterial: Parameters ---------- - k1: float - Principal permeability in local direction e1. - k2: float - Principal permeability in local direction e2. - k3: float - Principal permeability in local direction e3. + name: str + Name/label of the material. + k_vals: tuple[float, float, float] + Permeability values in principal directions. porosity: float Porosity of the material (between 0 and 1). thickness: float Thickness of the material in the out-of-plane direction. - name: str - Name/label of the material. """ - def __init__(self, k1:float, k2:float, k3:float, porosity:float, thickness:float, name:str = "unnamed_material"): - self.k_diag = np.array([[k1, 0, 0],[0, k2, 0],[0, 0, k3]]) + def __init__(self, name:str, k_vals : tuple[float, float, float], porosity:float, thickness:float): + self.is_isotropic = np.allclose([k_vals[0], k_vals[1], k_vals[2]], k_vals[0], atol=1e-14, rtol=0) + self.k_princ = np.diag(k_vals) self.porosity = porosity self.thickness = thickness self.name = name - self.is_isotropic = np.allclose([k1, k2, k3], k1, atol=1e-14, rtol=0) self.assigned = False + +class Resin: + __slots__ = ("name", "mu") + """Resin defined by dynamic viscosity (constant). + + Parameters + ---------- + name: str + Name of the resin. + mu: float + Dynamic viscosity of the resin [Pa.s] + """ + def __init__(self, name:str, mu:float): + self.name = name + self.mu = mu \ No newline at end of file diff --git a/src/lizzy/_core/materials/rosette.py b/src/lizzy/_core/materials/rosette.py index 0f79cf7..1b761a5 100644 --- a/src/lizzy/_core/materials/rosette.py +++ b/src/lizzy/_core/materials/rosette.py @@ -11,16 +11,14 @@ class Rosette: Parameters ---------- + name : str + The unique name of the rosette. p1 : tuple[float, float, float] - The first point defining the first axis of the rosette (k1 direction). - p0 : tuple[float, float, float] - The second point defining the first axis of the rosette (k1 direction). Default is (0,0,0). + The vector defining the first axis of the rosette (k1 direction). """ - def __init__(self, p1=(1.0,0,0), p0=(0.0,0.0,0.0), name:str=None): - p1 = np.array(p1) - p0 = np.array(p0) - self.u = p1 - p0 + def __init__(self, name:str, u=(1.0,0,0)): self.name = name + self.u = np.array(u) #TODO: this needs reviewing def project_along_normal(self, normal): diff --git a/src/lizzy/_core/solver/__init__.py b/src/lizzy/_core/solver/__init__.py index 94bcfb1..81e480d 100644 --- a/src/lizzy/_core/solver/__init__.py +++ b/src/lizzy/_core/solver/__init__.py @@ -6,4 +6,5 @@ from . import fem as fe from .psolvers import SolverType -from .solver import Solver \ No newline at end of file +from .solver import Solver +from .fillsolver import FillSolver \ No newline at end of file diff --git a/src/lizzy/_core/solver/builtin/iter_solvers.py b/src/lizzy/_core/solver/builtin/iter_solvers.py index 263c7fb..2811862 100644 --- a/src/lizzy/_core/solver/builtin/iter_solvers.py +++ b/src/lizzy/_core/solver/builtin/iter_solvers.py @@ -236,7 +236,7 @@ def solve_pressure_petsc(k: np.ndarray, f: np.ndarray, tol: float = 1e-8, max_it ksp.solve(b, x) if verbose: - print(f"PETSc solver converged in {ksp.getIterationNumber()} iterations " + print(f"\nPETSc solver converged in {ksp.getIterationNumber()} iterations " f"to a tolerance of {ksp.getResidualNorm():.2e}") # Convert solution back to numpy array diff --git a/src/lizzy/_core/solver/fillsolver.py b/src/lizzy/_core/solver/fillsolver.py index 434cf78..b5a2c0d 100644 --- a/src/lizzy/_core/solver/fillsolver.py +++ b/src/lizzy/_core/solver/fillsolver.py @@ -13,7 +13,7 @@ def __init__(self): self.map_cv_id_to_support_triangle_ids = {} self.map_cv_id_to_flux_terms = {} - def find_free_surface_cvs(self, fill_factor_array : np.ndarray, cv_support_cvs_array : np.ndarray): + def find_free_surface_cvs(self, fill_factor_array : np.ndarray, cv_support_cvs_array): """ Finds the control volumes that are on the flow front. These cvs have a fill factor < 1. """ diff --git a/src/lizzy/_core/solver/preprocessor.py b/src/lizzy/_core/solver/preprocessor.py new file mode 100644 index 0000000..056aefc --- /dev/null +++ b/src/lizzy/_core/solver/preprocessor.py @@ -0,0 +1,106 @@ +# Copyright 2025-2026 Simone Bancora, Paris Mulye +# +# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with this program. If not, see . + +from __future__ import annotations +from typing import TYPE_CHECKING + +from lizzy._core.cvmesh import mesh +if TYPE_CHECKING: + from lizzy._core.sensors import SensorManager + from lizzy._core.bcond import GatesManager + from lizzy._core.cvmesh import Mesh + from lizzy._core.materials import MaterialManager, Rosette, PorousMaterial + from lizzy._core.datatypes import SimulationParameters + +import sys +import numpy as np +import time +from lizzy._core.solver import * +from .timestep_manager import TimeStepManager +from .vsolvers import VelocitySolver +from .fillsolver import FillSolver +from .psolvers import PressureSolver, SolverType + + + +class Preprocessor: + def __init__(self, mesh:Mesh, fill_solver:FillSolver, vsolver:VelocitySolver, material_manager:MaterialManager, simulation_parameters:SimulationParameters): + self.mesh = mesh + self.fill_solver = fill_solver + self.vsolver = vsolver + self.material_manager = material_manager + self.simulation_parameters = simulation_parameters + + # SOLVER has already instantiated fill solver, velocity solver (ONCE) + + # SEQUENCE: + + + + # 1. check things were assigned + def assignment_checks(self): + if not self.simulation_parameters.has_been_assigned: + print(f"Warning: Simulation parameters were not assigned. Running with default values: wo_delta_time={self.simulation_parameters.wo_delta_time}") + if self.material_manager._assigned_resin == None: + print(f"ERROR: No resin assigned to the model. Create a resin using `LizzyModel.create_resin` and assign it using `LizzyModel.assign_resin`") + sys.exit(1) + + # 2. assign materials to elements + def assign_materials_to_elements(self): + materials = self.material_manager.assigned_materials + rosettes = self.material_manager.assigned_rosettes + for tri in self.mesh.triangles: + try: + material : PorousMaterial = materials[tri.material_tag] + if material.is_isotropic: + tri.k = material.k_princ + else: + rosette : Rosette = rosettes[tri.material_tag] + u, v, w = rosette.project_along_normal(tri.n) + R = np.array([u, v, w]).T + tri.k = R @ material.k_princ @ R.T + tri.porosity = materials[tri.material_tag].porosity + tri.h = materials[tri.material_tag].thickness + except KeyError: + exit(f"Mesh contains unassigned material tag: {tri.material_tag}") + + # 3. setup control volumes + def setup_cvs(self): + cvs = self.mesh.CVs + n_cvs = len(cvs) + node_idx_to_flux_ndarray: list[np.ndarray] = [None]*n_cvs + for i in range(n_cvs): + cvs[i].calculate_area_and_volume() + node_idx_to_flux_ndarray[i] = cvs[i].compute_flux_terms() + self.mesh.mesh_view.node_idx_to_flux_ndarray = node_idx_to_flux_ndarray + + # 4. assign data to fill solver + def assign_fill_solver_maps(self): + self.fill_solver.map_cv_id_to_support_triangle_ids = self.mesh.mesh_view.node_idx_to_tri_idxs + self.fill_solver.map_cv_id_to_flux_terms = self.mesh.mesh_view.node_idx_to_flux_ndarray + + # 5. assemble global stiffness matrix (singular) + def assemble_global_stiffnes_matrix(self): + mu = self.material_manager.assigned_resin.mu + K_sing, f_orig = fe.Assembly(self.mesh, mu, sparse=True) + return K_sing, f_orig + + def run_preproc_sequence(self): + print("Preprocessing...") + self.assignment_checks() + self.assign_materials_to_elements() + self.setup_cvs() + self.assign_fill_solver_maps() + K_sing, f_orig = self.assemble_global_stiffnes_matrix() + self.vsolver.precalculate_darcy_operator(self.mesh.triangles) + return K_sing, f_orig + + + + + + + diff --git a/src/lizzy/_core/solver/psolvers.py b/src/lizzy/_core/solver/psolvers.py index 9d1bfc7..771f741 100644 --- a/src/lizzy/_core/solver/psolvers.py +++ b/src/lizzy/_core/solver/psolvers.py @@ -29,7 +29,7 @@ class SolverType(Enum): class PressureSolver: @staticmethod - def solve(k:np.ndarray, f:np.ndarray, method:SolverType = SolverType.DIRECT_SPARSE, + def solve(k:np.ndarray, f:np.ndarray, method:SolverType, tol:float = 1e-8, max_iter:int = 1000, verbose:bool = False, **solver_kwargs): """ Solve the system `K p = f`. diff --git a/src/lizzy/_core/solver/solver.py b/src/lizzy/_core/solver/solver.py index 68103d3..7dfd4f9 100644 --- a/src/lizzy/_core/solver/solver.py +++ b/src/lizzy/_core/solver/solver.py @@ -10,43 +10,61 @@ from lizzy._core.cvmesh import mesh if TYPE_CHECKING: from lizzy._core.sensors import SensorManager - from lizzy._core.bcond import BCManager + from lizzy._core.bcond import GatesManager from lizzy._core.cvmesh import Mesh + from lizzy._core.materials import MaterialManager import sys import numpy as np import time from lizzy._core.solver import * -from lizzy._core.bcond import SolverBCs from .timestep_manager import TimeStepManager from .vsolvers import VelocitySolver from .fillsolver import FillSolver from .psolvers import PressureSolver, SolverType +from .preprocessor import Preprocessor -# from scipy.sparse import lil_matrix -# import matplotlib.pyplot as plt +class SolverBCs: + __slots__ = ("dirichlet_idx", "dirichlet_vals", "p0_idx") + + def __init__(self): + self.dirichlet_idx = np.empty(0, dtype=np.uint32) + self.dirichlet_vals = np.empty(0, dtype=np.float64) + self.p0_idx = np.empty(0, dtype=np.uint32) class Solver: - def __init__(self, mesh, bc_manager, simulation_parameters, material_manager, sensor_manager:SensorManager, - solver_type=SolverType.DIRECT_SPARSE, solver_tol=1e-8, solver_max_iter=1000, + def __init__(self, mesh:Mesh, gates_manager, simulation_parameters, material_manager:MaterialManager, sensor_manager:SensorManager, + solver_type=SolverType.ITERATIVE_PETSC, solver_tol=1e-8, solver_max_iter=1000, solver_verbose=False, use_masked_solver=True, **solver_kwargs): + + # create / assign all core components self.mesh : Mesh = mesh - self.bc_manager : BCManager = bc_manager - self.simulation_parameters = simulation_parameters + self.fill_solver = FillSolver() self.material_manager = material_manager - self.time_step_manager = TimeStepManager(mesh.nodes.N, mesh.triangles.N) + self.simulation_parameters = simulation_parameters + self.vsolver = VelocitySolver(self.mesh.triangles) + self.preproc = Preprocessor(mesh, self.fill_solver, self.vsolver, material_manager, simulation_parameters) + + self.gates_manager : GatesManager = gates_manager + self.time_step_manager = TimeStepManager(mesh.mesh_view.n_nodes, mesh.mesh_view.n_triangles) self._sensor_manager = sensor_manager self.bcs = SolverBCs() - self.vsolver = None - self.fill_solver = None self.solver_type = solver_type + if solver_type == SolverType.ITERATIVE_PETSC: + try: + import petsc4py + petsc4py.init() + from petsc4py import PETSc + except ImportError: + print("Import Error: PETSc not available. Reverting to DIRECT_SPARSE builtin solver.") + self.solver_type = SolverType.DIRECT_SPARSE self.solver_tol = solver_tol self.solver_max_iter = solver_max_iter self.solver_verbose = solver_verbose self.use_masked_solver = use_masked_solver self.solver_kwargs = solver_kwargs - self.N_nodes = mesh.nodes.N + self.N_nodes = mesh.mesh_view.n_nodes self.K_sing = None self.f_orig = None self.K_sol = None @@ -59,58 +77,52 @@ def __init__(self, mesh, bc_manager, simulation_parameters, material_manager, se self.step_completed = False self.k_local_all = np.empty((self.mesh.triangles.N, 6)) self.f_local_all = np.zeros((self.mesh.triangles.N, 3)) - self.solver_vars = {"fill_factor_array" : np.empty(self.N_nodes), + self.solver_vars = {"fill_factor_array" : np.zeros(self.N_nodes, dtype=float), "filled_node_ids" : np.empty(self.N_nodes, dtype=int), "free_surface_array" : np.empty(self.N_nodes), "cv_volumes_array" : np.empty(self.N_nodes),} # self.cv_adj_matrix = lil_matrix((self.N_nodes, self.N_nodes), dtype=int) - self.cv_support_cvs_array = {} - self.perform_fe_precalcs() - self.compute_k_local() + self.cv_support_cvs_array = self.mesh.mesh_view.node_idx_to_node_idxs # TODO do cleaner + + self.perform_precalcs() + + # self.compute_k_local() # when a solver is instantiated, all simulation variables are initialised self.initialise_new_solution() + - def perform_fe_precalcs(self): - if not self.mesh.preprocessed: - self.fill_solver = FillSolver() - self.mesh.preprocess(self.material_manager, self.fill_solver) - if not self.simulation_parameters.has_been_assigned: - print(f"Warning: Simulation parameters were not assigned. Running with default values: mu={self.simulation_parameters.mu}, wo_delta_time={self.simulation_parameters.wo_delta_time}") - # assemble FE global matrix (singular) - use sparse for better performance - self.K_sing, self.f_orig = fe.Assembly(self.mesh, self.simulation_parameters.mu, sparse=True) - # TODO: reorder nodes here to reduce bandwidth - then reorder the whole mesh and objects + def perform_precalcs(self): + self.K_sing, self.f_orig = self.preproc.run_preproc_sequence() # TODO: reorder nodes here to reduce bandwidth - then reorder the whole mesh and objects + self.vectorize_solver_vars() + self.initialise_sensor_manager() # could move into preprocessor as this runs only once + + def initialise_sensor_manager(self): + # assign sensors + self._sensor_manager.initialise(self.mesh) - self.vsolver = VelocitySolver(self.mesh.triangles) + def vectorize_solver_vars(self): # precalculate vectorised version of all variables - fill_factor_list = [] + # TODO this part can go faster cv_volumes_list = [] for cv in self.mesh.CVs: - fill_factor_list.append(cv.fill) cv_volumes_list.append(cv.vol) - self.cv_support_cvs_array[cv.idx] = np.array([support_cv.idx for support_cv in cv.support_CVs]) - - # # construct a cv adjacency matrix - # for key, item in self.cv_support_cvs_array.items(): - # self.cv_adj_matrix[key, item] = 1 - # self.cv_adj_matrix = self.cv_adj_matrix.tocsr() - self.solver_vars["fill_factor_array"] = np.array(fill_factor_list, dtype=float) self.solver_vars["cv_volumes_array"] = np.array(cv_volumes_list, dtype=float) - # assign sensors - self._sensor_manager.initialise(self.mesh) + def update_dirichlet_bcs(self): + # TODO this is more "update inlet dirichlet bcs" since it only applies pressure (doesn't add empty 0 pressure). It can be faster, but it doesn't run often (only at beginning of time intervals) so it's not critical dirichlet_idx = [] dirichlet_vals = [] - for tag, inlet in self.bc_manager._assigned_inlets.items(): + for boundary_name, inlet in self.gates_manager._assigned_inlets.items(): try: - inlet_idx = self.mesh.boundaries[tag] + inlet_idx = self.mesh.mesh_view.phys_boundary_name_to_node_idxs[boundary_name] except KeyError: print("\nFatal error: The application has terminated.") - print(f"Mesh does not contain physical tag: {tag}") + print(f"Mesh does not contain physical tag: {boundary_name}") sys.exit(1) if inlet.is_open: dirichlet_idx.append(inlet_idx) - dirichlet_vals.append(np.ones(len(inlet_idx)) * inlet.p_value) + dirichlet_vals.append(np.full(len(inlet_idx), inlet.p_value, dtype=np.float64)) try: self.bcs.dirichlet_idx = np.concatenate(dirichlet_idx) self.bcs.dirichlet_vals = np.concatenate(dirichlet_vals) @@ -119,13 +131,11 @@ def update_dirichlet_bcs(self): print("No inlets are currently open. At least one inlet must be open at all times to allow resin to flow into the part.") sys.exit(1) - def update_empty_nodes_idx(self): + def get_empty_nodes_idx(self, fill_factor): """ Complementary to "update_dirichlet_bcs()", this updates the indices of all nodes with a fill factor < 1.0. These will be uses to assign an internal condition p=0. """ - # empty_node_ids = [cv.idx for cv in self.mesh.CVs if cv.fill < 1] # nodes with fill factor < 1 - empty_node_ids = np.where(self.solver_vars["fill_factor_array"] < 1.0)[0] - self.bcs.p0_idx = np.array(empty_node_ids) + return np.where(fill_factor < 1.0)[0] def fill_initial_cvs(self): """ @@ -136,11 +146,7 @@ def fill_initial_cvs(self): # for cv in initial_cvs: # cv.fill = 1.0 - def update_n_empty_cvs(self): - """ - Must be called AFTER calling "update_empty_nodes_idx()" - """ - self.n_empty_cvs = len(self.bcs.p0_idx) + def generate_initial_time_step(self): time_0 = 0 @@ -168,13 +174,13 @@ def initialise_new_solution(self): self.solver_vars["fill_factor_array"] = np.zeros(self.N_nodes) self.bcs = SolverBCs() self.mesh.empty_cvs() - self.bc_manager.reset_inlets() + self.gates_manager.reset_inlets() self.update_dirichlet_bcs() self.fill_initial_cvs() - self.update_empty_nodes_idx() - self.update_n_empty_cvs() + p0_idxs = self.get_empty_nodes_idx(self.solver_vars["fill_factor_array"]) + self.n_empty_cvs = len(p0_idxs) + self.bcs.p0_idx = p0_idxs # self.K_sol, self.f_sol = PressureSolver.apply_starting_bcs(self.K_sing, self.f_orig, self.bcs) - self.new_step_dofs = [] self.solver_vars["filled_node_ids"] = np.where(self.solver_vars["fill_factor_array"] >= 1)[0] active_cvs_ids, self.solver_vars["free_surface_array"] = self.fill_solver.find_free_surface_cvs( self.solver_vars["fill_factor_array"], self.cv_support_cvs_array) @@ -214,37 +220,11 @@ def handle_wo_by_sensor_triggered(self, current_write_out, fill_factor_array): # print("\nSensor triggered") return write_out - def compute_k_local(self): - for i, tri in enumerate(self.mesh.triangles): - mu = self.simulation_parameters.mu - k_el = tri.grad_N.T @ tri.k @ tri.grad_N * tri.A * tri.h / mu - self.k_local_all[i, 0] = k_el[0,0] - self.k_local_all[i, 1] = k_el[1,1] - self.k_local_all[i, 2] = k_el[2,2] - self.k_local_all[i, 3] = k_el[0,1] - self.k_local_all[i, 4] = k_el[0,2] - self.k_local_all[i, 5] = k_el[1,2] - - def update_and_collect_solver_input(self): - dirichlet_idx_full = np.concatenate((self.bcs.dirichlet_idx, self.bcs.p0_idx), axis=None) - dirichlet_vals_full = np.concatenate((self.bcs.dirichlet_vals, np.zeros((1, len(self.bcs.p0_idx)))), axis=None) - - mask_nodes = self.solver_vars["free_surface_array"].copy() - mask_nodes[self.solver_vars["filled_node_ids"]] = 1 - - elem_connectivity = self.mesh.mesh_data["nodes_conn"] - filled_node_ids = self.solver_vars["filled_node_ids"] - - node_is_filled = np.zeros(self.mesh.nodes.N, dtype=bool) - node_is_filled[filled_node_ids] = True - - elem_filled_status = node_is_filled[elem_connectivity] # shape (T, 3) - mask_elements = np.any(elem_filled_status, axis=1).astype(int) - return self.k_local_all, self.f_local_all, dirichlet_idx_full, dirichlet_vals_full, mask_nodes, mask_elements, self.new_step_dofs, elem_connectivity - - def solve_time_step(self): # Use masked solver (optimized) or traditional solver + fill_factor = self.solver_vars["fill_factor_array"] + free_surface = self.solver_vars["free_surface_array"] + cv_volumes = self.solver_vars["cv_volumes_array"] if self.use_masked_solver: p = PressureSolver.solve_with_mask( self.K_sing, self.f_orig, self.bcs, @@ -257,34 +237,29 @@ def solve_time_step(self): max_iter=self.solver_max_iter, verbose=self.solver_verbose, **self.solver_kwargs) - v_array = self.vsolver.calculate_elem_velocities(p, self.simulation_parameters.mu) + v_array = self.vsolver.calculate_elem_velocities(p, self.material_manager.assigned_resin.mu) # v_nodal_array = self.vsolver.calculate_nodal_velocities(self.mesh.nodes, v_array) v_nodal_array = np.zeros((self.N_nodes, 3)) - active_cvs_ids, self.solver_vars["free_surface_array"] = self.fill_solver.find_free_surface_cvs(self.solver_vars["fill_factor_array"], self.cv_support_cvs_array) - dt = self.fill_solver.calculate_time_step(active_cvs_ids, self.solver_vars["fill_factor_array"], self.solver_vars["cv_volumes_array"], v_array) + active_cvs_ids, free_surface = self.fill_solver.find_free_surface_cvs(fill_factor, self.cv_support_cvs_array) + dt = self.fill_solver.calculate_time_step(active_cvs_ids, fill_factor, cv_volumes, v_array) dt, write_out = self.handle_wo_criterion(dt) - self.solver_vars["fill_factor_array"] = self.fill_solver.fill_current_time_step(active_cvs_ids, self.solver_vars["fill_factor_array"], self.solver_vars["cv_volumes_array"], dt, self.simulation_parameters.fill_tolerance) - - # find the newly filled cv ids as difference from the previous step - current_filled_node_ids = np.where(self.solver_vars["fill_factor_array"] >= 1)[0] - self.new_step_dofs = [id for id in current_filled_node_ids if id not in self.solver_vars["filled_node_ids"]] - self.solver_vars["filled_node_ids"] = current_filled_node_ids + fill_factor = self.fill_solver.fill_current_time_step(active_cvs_ids, fill_factor, cv_volumes, dt, self.simulation_parameters.fill_tolerance) # Update the filling time self.current_time += dt - # save time step results + if self.simulation_parameters.end_step_when_sensor_triggered: - write_out = self.handle_wo_by_sensor_triggered(write_out, self.solver_vars["fill_factor_array"]) - self.time_step_manager.save_timestep(self.current_time, dt, p, v_array, v_nodal_array, self.solver_vars["fill_factor_array"], self.solver_vars["free_surface_array"], write_out) + write_out = self.handle_wo_by_sensor_triggered(write_out, fill_factor) + self.time_step_manager.save_timestep(self.current_time, dt, p, v_array, v_nodal_array, fill_factor, free_surface, write_out) if write_out: # update sensors - self._sensor_manager.probe_current_solution(p, v_nodal_array, self.solver_vars["fill_factor_array"], self.current_time) - # update the empty nodes for next step - self.update_empty_nodes_idx() - # Print number of empty cvs - self.update_n_empty_cvs() + self._sensor_manager.probe_current_solution(p, v_nodal_array, fill_factor, self.current_time) + # update the empty nodes idxs and count for next step + p0_idxs = self.get_empty_nodes_idx(fill_factor) + self.n_empty_cvs = len(p0_idxs) + self.bcs.p0_idx = p0_idxs def solve(self, log="on") -> dict: solve_time_start = time.time() @@ -296,7 +271,7 @@ def solve(self, log="on") -> dict: while self.n_empty_cvs > 0: self.solve_time_step() if log == "on": - print("\rFill time: {:.5f}".format(self.current_time) + ", Empty CVs: {:4}".format(self.n_empty_cvs), end='') + print("\rFill time: {:.2f}".format(self.current_time) + "s, Empty CVs: {:4}".format(self.n_empty_cvs), end='') solution = self.time_step_manager.pack_solution() # good night and good luck solve_time_end = time.time() @@ -313,7 +288,7 @@ def solve_time_interval(self, time_interval:float, log="off", lightweight=False) self.update_dirichlet_bcs() self.solve_time_step() if log == "on": - print("\rFill time: {:.5f}".format(self.current_time) + ", Empty CVs: {:4}".format(self.n_empty_cvs), + print("\rFill time: {:.2f}".format(self.current_time) + "s, Empty CVs: {:4}".format(self.n_empty_cvs), end='') if lightweight: solution = "Lightweight mode: no solution is saved" @@ -324,3 +299,10 @@ def solve_time_interval(self, time_interval:float, log="off", lightweight=False) # print("\nSTEP SOLVE COMPLETED in {:.2f} seconds".format(total_solve_time)) return solution + + +#was after `fill_current_time_step` +# find the newly filled cv ids as difference from the previous step +# current_filled_node_ids = np.where(self.solver_vars["fill_factor_array"] >= 1)[0] +# self.new_step_dofs = [id for id in current_filled_node_ids if id not in self.solver_vars["filled_node_ids"]] +# self.solver_vars["filled_node_ids"] = current_filled_node_ids \ No newline at end of file diff --git a/src/lizzy/_core/solver/vsolvers.py b/src/lizzy/_core/solver/vsolvers.py index e777c3f..c989bc2 100644 --- a/src/lizzy/_core/solver/vsolvers.py +++ b/src/lizzy/_core/solver/vsolvers.py @@ -12,7 +12,7 @@ class VelocitySolver: def __init__(self, triangles): self.darcy_operator = any self.nodes_conn = any - self.precalculate_darcy_operator(triangles) + # self.precalculate_darcy_operator(triangles) def precalculate_darcy_operator(self, triangles): """precalculate vectorised coefficient darcy_operator of shape function gradients for velocity: v = darcy_operator * p""" diff --git a/src/lizzy/materials/__init__.py b/src/lizzy/materials/__init__.py index 9aa219d..06648e3 100644 --- a/src/lizzy/materials/__init__.py +++ b/src/lizzy/materials/__init__.py @@ -1 +1 @@ -from lizzy._core.materials import PorousMaterial, Rosette \ No newline at end of file +from lizzy._core.materials import PorousMaterial, Rosette, Resin \ No newline at end of file diff --git a/src/lizzy/model/model.py b/src/lizzy/model/model.py index 70b5e1e..4aef759 100644 --- a/src/lizzy/model/model.py +++ b/src/lizzy/model/model.py @@ -8,21 +8,22 @@ from typing import TYPE_CHECKING if TYPE_CHECKING: from lizzy._core.sensors import Sensor - from lizzy._core.materials import PorousMaterial, Rosette + from lizzy._core.materials import PorousMaterial, Rosette, Resin from lizzy._core.bcond.gates import Inlet from lizzy.datatypes import Solution + from typing import Dict, Literal from types import MappingProxyType from lizzy._core.io import Reader, Writer from lizzy._core.cvmesh import Mesh -from lizzy._core.bcond import BCManager +from lizzy._core.bcond import GatesManager from lizzy._core.solver import Solver, SolverType from lizzy._core.sensors import SensorManager from lizzy._core.datatypes import SimulationParameters from lizzy._core.materials import MaterialManager from lizzy.utils.splash_logo import print_logo - +from lizzy.utils.decorators import State, preinit_only, postinit_only class LizzyModel: """ @@ -30,18 +31,28 @@ class LizzyModel: """ def __init__(self): print_logo() - self._model_name : str = None - self._reader : Reader = None - self._writer : Writer = None - self._mesh : Mesh = None - self._solver : Solver = None + self._model_name:str = None + self._reader:Reader = None + self._writer:Writer = None + self._simulation_parameters:SimulationParameters = None + self._material_manager:MaterialManager = None + self._gates_manager:GatesManager = None + self._sensor_manager:SensorManager = None + self._mesh:Mesh = None + self._solver:Solver = None self._renderer : any = None self._latest_solution: Solution = None + self._lightweight:bool = False + self._state:State = State.PRE_INIT + self._create_components() + + def _create_components(self): + self._reader = Reader() + self._writer = Writer() self._simulation_parameters = SimulationParameters() self._material_manager = MaterialManager() - self._bc_manager = BCManager() + self._gates_manager = GatesManager() self._sensor_manager = SensorManager() - self._lightweight = False @property def lightweight(self): @@ -86,14 +97,14 @@ def current_time(self) -> float: return self._solver.current_time @property - def latest_solution(self) -> dict: + def latest_solution(self) -> Solution: """The most recent solution from the model (read-only). This value is None is the model is run in `lightweight` mode. """ return self._latest_solution @property - def bc_manager(self) -> BCManager: - return self._bc_manager + def gates_manager(self) -> GatesManager: + return self._gates_manager @property def simulation_parameters(self) -> SimulationParameters: @@ -123,6 +134,7 @@ def get_current_time(self) -> float: """ return self.current_time + @postinit_only def get_latest_solution(self) -> Dict: """ Returns @@ -132,21 +144,25 @@ def get_latest_solution(self) -> Dict: """ return self.latest_solution - + @preinit_only def assign_simulation_parameters(self, **kwargs): r""" Assigns new values to one or more simulation parameters using keyword arguments. Parameters ---------- - **kwargs : dict + **kwargs Keyword arguments corresponding to parameter names and their new values. - Valid parameters are: + Valid keywords are: - - ``mu``: viscosity [Pa s] - - ``wo_delta_time``: interval of simulation time between solution write-outs [s]. Default: -1 (write-out every numerical time step) - - ``fill_tolerance``: tolerance on the fill factor to consider a CV as filled. Default: 0.01 - - ``end_step_when_sensor_triggered``: if True, ends current solution step and creates a write-out when a sensor changes state. Default: False + mu: float, optional + Viscosity [Pa s]. Default: 0.1 + wo_delta_time: float, optional + Interval of simulation time between solution write-outs [s]. Default: -1 (write-out every numerical time step) + fill_tolerance: float, optional + Tolerance on the fill factor to consider a CV as filled. Default: 0.01 + end_step_when_sensor_triggered: bool, optional + If True, ends current solution step and creates a write-out when a sensor changes state. Default: False Examples -------- @@ -160,6 +176,7 @@ def assign_simulation_parameters(self, **kwargs): """ self._simulation_parameters.assign(**kwargs) + @preinit_only def read_mesh_file(self, mesh_file_path:str): r""" Reads a mesh file and initialises the mesh. Currently only .MSH format is supported (Version 4 ASCII). @@ -169,10 +186,10 @@ def read_mesh_file(self, mesh_file_path:str): mesh_file_path : str Path to the mesh file from the current working folder. """ - self._reader = Reader(mesh_file_path) + self._reader.read_mesh_file(mesh_file_path) self._model_name = self._reader.case_name self._mesh = Mesh(self._reader) - self._writer = Writer(self._mesh) + def print_mesh_info(self) -> None: """Prints some information about the mesh. @@ -182,11 +199,14 @@ def print_mesh_info(self) -> None: return self._reader.print_mesh_info() - def create_material(self, k1: float, k2: float, k3: float, porosity: float, thickness: float, name:str= None) -> PorousMaterial: + @preinit_only + def create_material(self, name : str, k_vals : tuple[float, float, float], porosity: float, thickness: float) -> PorousMaterial: """Create a new material that can then be selected and used in the model. Parameters ---------- + name : str + Unique name of the new material. k1 : float Permeability in the first principal direction. k2 : float @@ -197,17 +217,35 @@ def create_material(self, k1: float, k2: float, k3: float, porosity: float, thic Volumetric porosity of the material (porosity = 1 - fibre volume fraction). thickness : float Thickness of the material [mm]. - name : str, optional - Label assigned to the material. Necessary to select the material during assignment. If none assigned, a default 'Material_{N}'name is given, where N is an incremental number of existing materials. - + Returns ------- :class:`~lizzy.core.materials.PorousMaterial` - Instance of the created material. + Reference to the created material. """ - new_material = self._material_manager.create_material(k1, k2, k3, porosity, thickness, name) + new_material = self._material_manager.create_material(name, k_vals, porosity, thickness) return new_material + + @preinit_only + def create_resin(self, name : str, viscosity : float) -> Resin: + """Create a new resin that can then be selected and used in the model. + Parameters + ---------- + name : str + Unique name of the resin. + viscosity : float + Dynamic viscosity of the resin [Pa.s] + + Returns + ------- + :class:`~lizzy.core.materials.Resin` + Reference to the created resin. + """ + new_resin = self._material_manager.create_resin(name, viscosity) + return new_resin + + @preinit_only def assign_material(self, material_selector, mesh_tag:str, rosette:Rosette = None): """Assign an existing material to a labeled mesh region. @@ -222,44 +260,56 @@ def assign_material(self, material_selector, mesh_tag:str, rosette:Rosette = Non """ self._material_manager.assign_material(material_selector, mesh_tag, rosette) - def create_rosette(self, p1:tuple[float, float, float], p0:tuple[float, float, float]=(0.0,0.0,0.0), name:str=None) -> Rosette: + @preinit_only + def assign_resin(self, resin_selector): + """Assign an existing resin to the model. + + Parameters + ---------- + material_selector : str + Name of the resin to assign. Must correspond to an existing resin created with `LizzyModel.create_resin`. + """ + self._material_manager.assign_resin(resin_selector) + + @preinit_only + def create_rosette(self, name:str, u:tuple[float, float, float]) -> Rosette: """Create a new rosette that can then be selected and used in the model. Parameters ---------- - p1 : tuple[float, float, float] + name : str + Name unique name of the new rosette. + u : tuple[float, float, float] The first point defining the first axis of the rosette (k1 direction). - p0 : tuple[float, float, float] - The second point defining the first axis of the rosette (k1 direction). Default is (0,0,0). - name : str, optional - Label assigned to the rosette. Necessary to select the rosette during assignment. If none assigned, a default 'Rosette_{N}'name is given, where N is an incremental number of existing rosettes. Returns ------- :class:`~lizzy.core.materials.Rosette` - Instance of the created rosette. + Reference to the created rosette. """ - new_rosette = self._material_manager.create_rosette(p1, p0, name) + new_rosette = self._material_manager.create_rosette(name, u) return new_rosette - def create_inlet(self, initial_pressure_value:float, name:str = None) -> Inlet: + @preinit_only + def create_inlet(self, name:str, initial_pressure_value:float) -> Inlet: """Creates a new inlet and add it to model existing inlets. Parameters ---------- + name : str + Name of the inlet. initial_pressure_value : float Initial pressure value at the inlet. - name : str, optional - Label assigned to the inlet. Will be used to select the inlet in future operations. If none assigned, a default 'Inlet_{N}'name is given, where N is an incremental number of existing inlets. Returns ------- :class:`~lizzy.core.bcond.Inlet` The created inlet object. """ - new_inlet = self._bc_manager.create_inlet(initial_pressure_value, name) + new_inlet = self._gates_manager.create_inlet(name, initial_pressure_value) return new_inlet + @preinit_only def assign_inlet(self, inlet_selector:Inlet | str, boundary_tag:str): """Selects an inlet from existing ones and assigns it to the indicated mesh boundary. @@ -270,7 +320,7 @@ def assign_inlet(self, inlet_selector:Inlet | str, boundary_tag:str): boundary_tag : str An existing mesh boundary tag where to assign the inlet. """ - self._bc_manager.assign_inlet(inlet_selector, boundary_tag) + self._gates_manager.assign_inlet(inlet_selector, boundary_tag) def fetch_inlet_by_name(self, inlet_name: str) -> Inlet: """Fetches an inlet from existing ones in the model. @@ -285,9 +335,10 @@ def fetch_inlet_by_name(self, inlet_name: str) -> Inlet: :class:`~lizzy.gates.Inlet` The fetched inlet object. """ - selected_inlet = self._bc_manager._fetch_inlet(inlet_name) + selected_inlet = self._gates_manager._fetch_inlet(inlet_name) return selected_inlet + def change_inlet_pressure(self, inlet_selector:Inlet | str, pressure_value:float, mode: Literal["set", "delta"] = "set"): """Changes the pressure value at the selected inlet to a new value, according to the selected mode. @@ -307,7 +358,7 @@ def change_inlet_pressure(self, inlet_selector:Inlet | str, pressure_value:float KeyError If the `mode` is not one of the allowed values. """ - self._bc_manager.change_inlet_pressure(inlet_selector, pressure_value, mode) + self._gates_manager.change_inlet_pressure(inlet_selector, pressure_value, mode) def open_inlet(self, inlet_selector:Inlet | str): """Sets the selected inlet state to `open`. When open, the inlet applies its p_value as a Dirichlet boundary condition. An inlet can be opened at any time during the simulation. @@ -317,7 +368,7 @@ def open_inlet(self, inlet_selector:Inlet | str): inlet_selector : Inlet | str Either the inlet object reference, or the name of an existing inlet. """ - self._bc_manager.open_inlet(inlet_selector) + self._gates_manager.open_inlet(inlet_selector) def close_inlet(self, inlet_selector:Inlet | str): """Sets the selected inlet state to `closed`. When closed, the inlet acts as a Neumann natural boundary condition (no flux). An inlet can be closed at any time during the simulation. The stored p_value is preserved when the inlet is closed. @@ -327,9 +378,10 @@ def close_inlet(self, inlet_selector:Inlet | str): inlet_selector : Inlet | str Either the inlet object reference, or the name of an existing inlet. """ - self._bc_manager.close_inlet(inlet_selector) + self._gates_manager.close_inlet(inlet_selector) #TODO: get coords arg as tuple or np array, then ids as int or string + @preinit_only def create_sensor(self, x:float, y:float, z:float): """Create a virtual sensor at the specified position and add it to the model. @@ -344,6 +396,7 @@ def create_sensor(self, x:float, y:float, z:float): """ self._sensor_manager.add_sensor(x, y, z) + def print_sensor_readings(self): """Prints to the console the current values of :attr:`~lizzy.sensors.Sensor.time`, :attr:`~lizzy.sensors.Sensor.pressure`, :attr:`~lizzy.sensors.Sensor.fill_factor` and :attr:`~lizzy.sensors.Sensor.velocity` of each sensor. """ @@ -358,7 +411,8 @@ def get_sensor_by_id(self, idx) -> Sensor: """ return self._sensor_manager.get_sensor_by_id(idx) - def initialise_solver(self, solver_type:SolverType = SolverType.DIRECT_SPARSE, + + def initialise_solver(self, solver_type:SolverType = SolverType.ITERATIVE_PETSC, solver_tol:float = 1e-8, solver_max_iter:int = 1000, solver_verbose:bool = False, use_masked_solver:bool = True, **solver_kwargs): @@ -369,7 +423,7 @@ def initialise_solver(self, solver_type:SolverType = SolverType.DIRECT_SPARSE, ---------- solver_type : SolverType Type of linear solver (DIRECT_DENSE, DIRECT_SPARSE, ITERATIVE_PETSC). - Default is DIRECT_SPARSE which is more efficient for sparse matrices. + Default is ITERATIVE_PETSC and will revert to DIRECT_SPARSE is PETSc is not installed. solver_tol : float Convergence tolerance for iterative solvers solver_max_iter : int @@ -378,16 +432,20 @@ def initialise_solver(self, solver_type:SolverType = SolverType.DIRECT_SPARSE, Print solver convergence information use_masked_solver : bool Use optimized masked solver (only solves for free DOFs). - Default is True for better performance (10-400x speedup in early timesteps). + Default is True for better performance. **solver_kwargs Additional solver-specific keyword arguments """ - self._solver = Solver(self._mesh, self._bc_manager, self._simulation_parameters, + + self._solver = Solver(self._mesh, self._gates_manager, self._simulation_parameters, self._material_manager, self._sensor_manager, solver_type, solver_tol, solver_max_iter, solver_verbose, use_masked_solver, **solver_kwargs) + + self._state = State.POST_INIT - def solve(self, log="on") -> dict: + @postinit_only + def solve(self, log="on") -> Solution: """Advance the filling simulation from the current time until the part is filled. Parameters @@ -397,13 +455,14 @@ def solve(self, log="on") -> dict: Returns ------- - solution : dict - A dictionary containing the solution. + solution : :class:`~lizzy.datatypes.Solution` + A Solution object storing the solution fields up to the time step reached """ self._latest_solution = self._solver.solve(log=log) return self._latest_solution - def solve_time_interval(self, time_interval:float, log="off") -> dict: + @postinit_only + def solve_time_interval(self, time_interval:float, log="off") -> Solution: """Advance the filling simulation from the current time for the specified time interval. Parameters @@ -415,18 +474,20 @@ def solve_time_interval(self, time_interval:float, log="off") -> dict: Returns ------- - solution : dict - A dictionary containing the solution up to the time step reached. + solution : :class:`~lizzy.datatypes.Solution` + A Solution object storing the solution fields up to the time step reached. """ self._latest_solution = self._solver.solve_time_interval(time_interval, log=log, lightweight=self._lightweight) return self._latest_solution + @postinit_only def initialise_new_solution(self): """ Initialises a new solution, resetting all simulation variables. The part will be emptied and initial boundary conditions restored. This method can be called to reset a simulation and run a new one, without resetting the model. """ self._solver.initialise_new_solution() + @postinit_only def save_results(self, solution: Solution = None, result_name:str = None, **kwargs): """Save the results contained in the solution dictionary into an XDMF file. @@ -441,6 +502,7 @@ def save_results(self, solution: Solution = None, result_name:str = None, **kwar solution = self._latest_solution if result_name == None: result_name = self._model_name + '_RES' + self._writer.assign_mesh(self._mesh) self._writer.save_results(solution, result_name, **kwargs) def get_node_by_id(self, node_id:int): diff --git a/src/lizzy/utils/decorators.py b/src/lizzy/utils/decorators.py new file mode 100644 index 0000000..d68b65f --- /dev/null +++ b/src/lizzy/utils/decorators.py @@ -0,0 +1,31 @@ +from __future__ import annotations +from typing import TYPE_CHECKING +if TYPE_CHECKING: + from lizzy.model import LizzyModel + +import sys +from functools import wraps +from enum import Enum + + +class State(Enum): + PRE_INIT = 0 + POST_INIT = 1 + +def preinit_only(method): + @wraps(method) + def wrapper(self:LizzyModel, *args, **kwargs): + if self._state != State.PRE_INIT: + print(f"ERROR: Method '{method.__name__}' must be called before `initialise_solver()`.") + sys.exit(1) + return method(self, *args, **kwargs) + return wrapper + +def postinit_only(method): + @wraps(method) + def wrapper(self:LizzyModel, *args, **kwargs): + if self._state != State.POST_INIT: + print(f"ERROR: Method '{method.__name__}' must be called after initialise_solver().") + sys.exit(1) + return method(self, *args, **kwargs) + return wrapper \ No newline at end of file diff --git a/tests/test_rect.py b/tests/test_rect.py index e40cd84..be01cca 100644 --- a/tests/test_rect.py +++ b/tests/test_rect.py @@ -12,14 +12,16 @@ def model(): model = liz.LizzyModel() model.read_mesh_file("tests/test_meshes/Rect_1M_R1.msh") - model.assign_simulation_parameters(mu=0.1, wo_delta_time=100, fill_tolerance=0.00) - model.create_material(1E-10, 1E-10, 1E-10, 0.5, 1.0, "test_material") + model.assign_simulation_parameters(wo_delta_time=100, fill_tolerance=0.00) + model.create_resin("resin", 0.1) + model.assign_resin("resin") + model.create_material("test_material", (1E-10, 1E-10, 1E-10), 0.5, 1.0) model.assign_material("test_material", 'domain') return model def test_fill_1bar(model: liz.LizzyModel): analytical_solution = 2500 - model.create_inlet(1E+05, "inlet_left") + model.create_inlet("inlet_left", 1E+05) model.assign_inlet("inlet_left", "left_edge") model.initialise_solver() solution = model.solve() @@ -28,7 +30,7 @@ def test_fill_1bar(model: liz.LizzyModel): def test_fill_01bar(model: liz.LizzyModel): analytical_solution = 25000 - model.create_inlet(1E+04, "inlet_left") + model.create_inlet("inlet_left", 1E+04) model.assign_inlet("inlet_left", "left_edge") model.initialise_solver() solution = model.solve()