Merge pull request #10 from jhdark/fenicsx

Fenicsx

.github/workflows/ci.yml

@@ -2,40 +2,40 @@ name: CI
 on: [pull_request, push]
 
 jobs:
-  run-on-conda:
-    runs-on: ubuntu-latest
-
-    steps:
-    - name: Checkout code
-      uses: actions/checkout@v2
-
-    - name: Set up Conda
-      uses: conda-incubator/setup-miniconda@v2
-      with:
-        activate-environment: myenv
-        mamba-version: "*"
-        channels: conda-forge, defaults
-
-    - name: Create Conda environment
-      shell: bash -l {0}
-      run: |
-        mamba install -c conda-forge fenics-dolfinx
-
-    - name: Install local package and dependencies
-      shell: bash -l {0}
-      run: |
-        pip install .[test]
-
-    - name: Run tests
-      shell: bash -l {0}
-      run: |
-        pytest test/ --cov festim --cov-report xml --cov-report term
+  # run-on-conda:
+  #   runs-on: ubuntu-latest
+
+  #   steps:
+  #   - name: Checkout code
+  #     uses: actions/checkout@v2
+
+  #   - name: Set up Conda
+  #     uses: conda-incubator/setup-miniconda@v2
+  #     with:
+  #       activate-environment: myenv
+  #       mamba-version: "*"
+  #       channels: conda-forge, defaults
+
+  #   - name: Create Conda environment
+  #     shell: bash -l {0}
+  #     run: |
+  #       mamba install -c conda-forge fenics-dolfinx
+
+  #   - name: Install local package and dependencies
+  #     shell: bash -l {0}
+  #     run: |
+  #       pip install .[test]
+
+  #   - name: Run tests
+  #     shell: bash -l {0}
+  #     run: |
+  #       pytest test/ --cov festim --cov-report xml --cov-report term
 
-    - name: Upload to codecov
-      uses: codecov/codecov-action@v3
-      with:
-        token: ${{ secrets.CODECOV_TOKEN }}
-        files: ./coverage.xml
+  #   - name: Upload to codecov
+  #     uses: codecov/codecov-action@v3
+  #     with:
+  #       token: ${{ secrets.CODECOV_TOKEN }}
+  #       files: ./coverage.xml
 
   run-on-docker:
     runs-on: ubuntu-latest

festim/__init__.py

@@ -29,3 +29,5 @@
 
 from .subdomain.surface_subdomain import SurfaceSubdomain1D
 from .subdomain.volume_subdomain import VolumeSubdomain1D
+
+from .exports.vtx import VTXExport

festim/exports/vtx.py

@@ -0,0 +1,83 @@
+from dolfinx.io import VTXWriter
+import mpi4py
+
+import festim as F
+
+
+class VTXExport:
+    """Export functions to VTX file
+
+    Args:
+        filename (str): the name of the output file
+        field (int): the field index to export
+
+    Attributes:
+        filename (str): the name of the output file
+        writer (dolfinx.io.VTXWriter): the VTX writer
+        field (festim.Species, list of festim.Species): the field index to export
+
+    Usage:
+        >>> u = dolfinx.fem.Function(V)
+        >>> my_export = festim.VTXExport("my_export.bp")
+        >>> my_export.define_writer(mesh.comm, [u])
+        >>> for t in range(10):
+        ...    u.interpolate(lambda x: t * (x[0] ** 2 + x[1] ** 2 + x[2] ** 2))
+        ...    my_export.write(t)
+    """
+
+    def __init__(self, filename: str, field) -> None:
+        self.filename = filename
+        self.field = field
+
+    @property
+    def filename(self):
+        return self._filename
+
+    @filename.setter
+    def filename(self, value):
+        if not isinstance(value, str):
+            raise TypeError("filename must be of type str")
+        if not value.endswith(".bp"):
+            raise ValueError("filename must end with .bp")
+        self._filename = value
+
+    @property
+    def field(self):
+        return self._field
+
+    @field.setter
+    def field(self, value):
+        # check that field is festim.Species or list of festim.Species
+        if not isinstance(value, F.Species) and not isinstance(value, list):
+            raise TypeError(
+                "field must be of type festim.Species or list of festim.Species"
+            )
+        # check that all elements of list are festim.Species
+        if isinstance(value, list):
+            for element in value:
+                if not isinstance(element, F.Species):
+                    raise TypeError(
+                        "field must be of type festim.Species or list of festim.Species"
+                    )
+        # if field is festim.Species, convert to list
+        if not isinstance(value, list):
+            value = [value]
+
+        self._field = value
+
+    def define_writer(self, comm: mpi4py.MPI.Intracomm, functions: list) -> None:
+        """Define the writer
+
+        Args:
+            comm (mpi4py.MPI.Intracomm): the MPI communicator
+            functions (list): the list of functions to export
+        """
+        self.writer = VTXWriter(comm, self.filename, functions, "BP4")
+
+    def write(self, t: float):
+        """Write functions to VTX file
+
+        Args:
+            t (float): the time of export
+        """
+        self.writer.write(t)

festim/hydrogen_transport_problem.py

@@ -1,11 +1,14 @@
 from dolfinx import fem
 from dolfinx.nls.petsc import NewtonSolver
+from dolfinx.io import XDMFFile
 import ufl
 from mpi4py import MPI
-from dolfinx.fem import Function
+from dolfinx.fem import Function, form, assemble_scalar
 from dolfinx.mesh import meshtags
-from ufl import TestFunction, dot, grad, Measure
+from ufl import TestFunction, dot, grad, Measure, FacetNormal
 import numpy as np
+import tqdm.autonotebook
+
 
 import festim as F
 
@@ -110,6 +113,18 @@ def initialise(self):
         self.define_markers_and_measures()
         self.assign_functions_to_species()
         self.create_formulation()
+        self.defing_export_writers()
+
+    def defing_export_writers(self):
+        """Defines the export writers of the model"""
+        for export in self.exports:
+            # TODO implement when export.field is an int or str
+            # then find solution from index of species
+
+            if isinstance(export, F.VTXExport):
+                export.define_writer(
+                    MPI.COMM_WORLD, [field.solution for field in export.field]
+                )
 
     def define_function_space(self):
         elements = ufl.FiniteElement("CG", self.mesh.mesh.ufl_cell(), 1)
@@ -171,6 +186,10 @@ def assign_functions_to_species(self):
             spe.prev_solution = Function(self.function_space)
             spe.test_function = TestFunction(self.function_space)
 
+        # TODO remove this
+        self.u = self.species[0].solution
+        self.u_n = self.species[0].prev_solution
+
     def create_formulation(self):
         """Creates the formulation of the model"""
         if len(self.sources) > 1:
@@ -218,3 +237,47 @@ def create_solver(self):
         )
         solver = NewtonSolver(MPI.COMM_WORLD, problem)
         self.solver = solver
+
+    def run(self, final_time: float):
+        """Runs the model for a given time
+
+        Args:
+            final_time (float): the final time of the simulation
+
+        Returns:
+            list of float: the times of the simulation
+            list of float: the fluxes of the simulation
+        """
+        t = 0
+        times, flux_values = [], []
+
+        mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w")
+        mobile_xdmf.write_mesh(self.mesh.mesh)
+
+        progress = tqdm.autonotebook.tqdm(
+            desc="Solving H transport problem", total=final_time
+        )
+        while t < final_time:
+            progress.update(float(self.dt))
+            t += float(self.dt)
+
+            self.solver.solve(self.u)
+
+            mobile_xdmf.write_function(self.u, t)
+
+            cm = self.species[0].solution
+            # TODO this should be a property of Mesh
+            n = FacetNormal(self.mesh.mesh)
+            D = self.subdomains[0].material.get_diffusion_coefficient(
+                self.mesh.mesh, self.temperature
+            )
+            surface_flux = form(D * dot(grad(cm), n) * self.ds(2))
+            flux = assemble_scalar(surface_flux)
+            flux_values.append(flux)
+            times.append(t)
+
+            # update previous solution
+            self.u_n.x.array[:] = self.u.x.array[:]
+
+        mobile_xdmf.close()
+        return times, flux_values

test/test_permeation_problem.py

@@ -1,16 +1,11 @@
-from mpi4py import MPI
 from petsc4py import PETSc
-from dolfinx.io import XDMFFile
 from dolfinx.fem import (
     Constant,
     dirichletbc,
     locate_dofs_topological,
-    form,
-    assemble_scalar,
 )
-from ufl import dot, grad, exp, FacetNormal
+from ufl import exp, FacetNormal
 import numpy as np
-import tqdm.autonotebook
 
 
 import festim as F
@@ -37,15 +32,13 @@ def test_permeation_problem():
     temperature = Constant(my_mesh.mesh, 500.0)
     my_model.temperature = temperature
 
+    my_model.exports = [F.VTXExport("test.bp", field=mobile_H)]
+
     my_model.initialise()
 
     D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)
 
     V = my_model.function_space
-    u = mobile_H.solution
-
-    # TODO this should be a property of Mesh
-    n = FacetNormal(my_mesh.mesh)
 
     def siverts_law(T, S_0, E_S, pressure):
         S = S_0 * exp(-E_S / F.k_B / T)
@@ -81,33 +74,9 @@ def siverts_law(T, S_0, E_S, pressure):
     opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
     ksp.setFromOptions()
 
-    mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w")
-    mobile_xdmf.write_mesh(my_model.mesh.mesh)
-
     final_time = 50
 
-    flux_values = []
-    times = []
-    t = 0
-    progress = tqdm.autonotebook.tqdm(
-        desc="Solving H transport problem", total=final_time
-    )
-    while t < final_time:
-        progress.update(float(my_model.dt))
-        t += float(my_model.dt)
-
-        my_model.solver.solve(u)
-
-        mobile_xdmf.write_function(u, t)
-
-        surface_flux = form(D * dot(grad(u), n) * my_model.ds(2))
-        flux = assemble_scalar(surface_flux)
-        flux_values.append(flux)
-        times.append(t)
-
-        mobile_H.prev_solution.x.array[:] = u.x.array[:]
-
-    mobile_xdmf.close()
+    times, flux_values = my_model.run(final_time=final_time)
 
     # analytical solution
     S = S_0 * exp(-E_S / F.k_B / float(temperature))
@@ -167,6 +136,8 @@ def test_permeation_problem_multi_volume():
     temperature = Constant(my_mesh.mesh, 500.0)
     my_model.temperature = temperature
 
+    my_model.exports = [F.VTXExport("test.bp", field=mobile_H)]
+
     my_model.initialise()
 
     D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)
@@ -211,33 +182,9 @@ def siverts_law(T, S_0, E_S, pressure):
     opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
     ksp.setFromOptions()
 
-    mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w")
-    mobile_xdmf.write_mesh(my_model.mesh.mesh)
-
     final_time = 50
 
-    flux_values = []
-    times = []
-    t = 0
-    progress = tqdm.autonotebook.tqdm(
-        desc="Solving H transport problem", total=final_time
-    )
-    while t < final_time:
-        progress.update(float(my_model.dt))
-        t += float(my_model.dt)
-
-        my_model.solver.solve(u)
-
-        mobile_xdmf.write_function(u, t)
-
-        surface_flux = form(D * dot(grad(u), n) * my_model.ds(2))
-        flux = assemble_scalar(surface_flux)
-        flux_values.append(flux)
-        times.append(t)
-
-        mobile_H.prev_solution.x.array[:] = u.x.array[:]
-
-    mobile_xdmf.close()
+    times, flux_values = my_model.run(final_time=final_time)
 
     # analytical solution
     S = S_0 * exp(-E_S / F.k_B / float(temperature))