From 0e0bd475cda848fac7d58192be44b33106d9c9e0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 13:16:46 -0400 Subject: [PATCH 1/9] normal as property --- festim/mesh/mesh.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py index bb2aa0654..e88b1b53f 100644 --- a/festim/mesh/mesh.py +++ b/festim/mesh/mesh.py @@ -1,3 +1,6 @@ +import ufl + + class Mesh: """ Mesh class @@ -32,3 +35,7 @@ def vdim(self): @property def fdim(self): return self.mesh.topology.dim - 1 + + @property + def n(self): + return ufl.FacetNormal(self.mesh) From 0920e568d0825e9468b97c57de5da8519189ce88 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 13:17:02 -0400 Subject: [PATCH 2/9] modify permeation test for new property --- test/test_permeation_problem.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index b8534107e..ab50725c9 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -8,7 +8,7 @@ form, assemble_scalar, ) -from ufl import dot, grad, exp, FacetNormal +from ufl import dot, grad, exp import numpy as np import tqdm.autonotebook @@ -44,9 +44,6 @@ def test_permeation_problem(): 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) return S * pressure**0.5 @@ -100,7 +97,7 @@ def siverts_law(T, S_0, E_S, pressure): mobile_xdmf.write_function(u, t) - surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) + surface_flux = form(D * dot(grad(u), my_mesh.n) * my_model.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) @@ -174,9 +171,6 @@ def test_permeation_problem_multi_volume(): 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) return S * pressure**0.5 @@ -230,7 +224,7 @@ def siverts_law(T, S_0, E_S, pressure): mobile_xdmf.write_function(u, t) - surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) + surface_flux = form(D * dot(grad(u), my_mesh.n) * my_model.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) From 3ad9c2706c99cb43d6085df8ef04f0285a09857c Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 08:21:52 -0400 Subject: [PATCH 3/9] updated doc strings --- festim/mesh/mesh.py | 1 + 1 file changed, 1 insertion(+) diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py index e88b1b53f..fc0980317 100644 --- a/festim/mesh/mesh.py +++ b/festim/mesh/mesh.py @@ -12,6 +12,7 @@ class Mesh: mesh (dolfinx.mesh.Mesh): the mesh vdim (int): the dimension of the mesh cells fdim (int): the dimension of the mesh facets + n (ufl.FacetNormal): the normal vector to the facets """ def __init__(self, mesh=None): From 2e8e011a5ce3a51f7bdc6e44a99167bd179c12b9 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 08:22:13 -0400 Subject: [PATCH 4/9] permeation test to use fdim property --- test/test_permeation_problem.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index ab50725c9..0cbc2adc5 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -48,11 +48,10 @@ def siverts_law(T, S_0, E_S, pressure): S = S_0 * exp(-E_S / F.k_B / T) return S * pressure**0.5 - fdim = my_mesh.mesh.topology.dim - 1 left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, fdim, left_facets) + left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, fdim, right_facets) + right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) S_0 = 4.02e21 E_S = 1.04 @@ -175,11 +174,10 @@ def siverts_law(T, S_0, E_S, pressure): S = S_0 * exp(-E_S / F.k_B / T) return S * pressure**0.5 - fdim = my_mesh.mesh.topology.dim - 1 left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, fdim, left_facets) + left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, fdim, right_facets) + right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) S_0 = 4.02e21 E_S = 1.04 From 5c52539f9d8cdb1bdcfc47fb6038054a5c6d9ae7 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:53:19 -0400 Subject: [PATCH 5/9] remove file --- test/test_permeation_problem.py | 256 -------------------------------- 1 file changed, 256 deletions(-) delete mode 100644 test/test_permeation_problem.py diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py deleted file mode 100644 index 0cbc2adc5..000000000 --- a/test/test_permeation_problem.py +++ /dev/null @@ -1,256 +0,0 @@ -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 -import numpy as np -import tqdm.autonotebook - - -import festim as F - - -def test_permeation_problem(): - L = 3e-04 - vertices = np.linspace(0, L, num=1001) - - my_mesh = F.Mesh1D(vertices) - - my_model = F.HydrogenTransportProblem() - my_model.mesh = my_mesh - - my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") - my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat) - left_surface = F.SurfaceSubdomain1D(id=1, x=0) - right_surface = F.SurfaceSubdomain1D(id=2, x=L) - my_model.subdomains = [my_subdomain, left_surface, right_surface] - - mobile_H = F.Species("H") - my_model.species = [mobile_H] - - temperature = Constant(my_mesh.mesh, 500.0) - my_model.temperature = temperature - - my_model.initialise() - - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) - - V = my_model.function_space - u = mobile_H.solution - - def siverts_law(T, S_0, E_S, pressure): - S = S_0 * exp(-E_S / F.k_B / T) - return S * pressure**0.5 - - left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) - right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) - - S_0 = 4.02e21 - E_S = 1.04 - P_up = 100 - surface_conc = siverts_law(T=temperature, S_0=S_0, E_S=E_S, pressure=P_up) - bc_sieverts = dirichletbc( - Constant(my_mesh.mesh, PETSc.ScalarType(surface_conc)), left_dofs, V - ) - bc_outgas = dirichletbc(Constant(my_mesh.mesh, PETSc.ScalarType(0)), right_dofs, V) - my_model.boundary_conditions = [bc_sieverts, bc_outgas] - my_model.create_solver() - - my_model.solver.convergence_criterion = "incremental" - my_model.solver.rtol = 1e-10 - my_model.solver.atol = 1e10 - - my_model.solver.report = True - ksp = my_model.solver.krylov_solver - opts = PETSc.Options() - option_prefix = ksp.getOptionsPrefix() - opts[f"{option_prefix}ksp_type"] = "cg" - opts[f"{option_prefix}pc_type"] = "gamg" - 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), my_mesh.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() - - # analytical solution - S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability = float(D) * S - times = np.array(times) - - n_array = np.arange(1, 10000)[:, np.newaxis] - summation = np.sum( - (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), - axis=0, - ) - analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) - - analytical_flux = np.abs(analytical_flux) - flux_values = np.array(np.abs(flux_values)) - - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - - relative_error = relative_error[ - np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - ] - error = relative_error.mean() - - assert error < 0.01 - - -def test_permeation_problem_multi_volume(): - L = 3e-04 - vertices = np.linspace(0, L, num=1001) - - my_mesh = F.Mesh1D(vertices) - - my_model = F.HydrogenTransportProblem() - my_model.mesh = my_mesh - - my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat") - my_subdomain_1 = F.VolumeSubdomain1D(id=1, borders=[0, L / 4], material=my_mat) - my_subdomain_2 = F.VolumeSubdomain1D(id=2, borders=[L / 4, L / 2], material=my_mat) - my_subdomain_3 = F.VolumeSubdomain1D( - id=3, borders=[L / 2, 3 * L / 4], material=my_mat - ) - my_subdomain_4 = F.VolumeSubdomain1D(id=4, borders=[3 * L / 4, L], material=my_mat) - left_surface = F.SurfaceSubdomain1D(id=1, x=0) - right_surface = F.SurfaceSubdomain1D(id=2, x=L) - my_model.subdomains = [ - my_subdomain_1, - my_subdomain_2, - my_subdomain_3, - my_subdomain_4, - left_surface, - right_surface, - ] - - mobile_H = F.Species("H") - my_model.species = [mobile_H] - - temperature = Constant(my_mesh.mesh, 500.0) - my_model.temperature = temperature - - my_model.initialise() - - D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature) - - V = my_model.function_space - u = mobile_H.solution - - def siverts_law(T, S_0, E_S, pressure): - S = S_0 * exp(-E_S / F.k_B / T) - return S * pressure**0.5 - - left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) - right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) - - S_0 = 4.02e21 - E_S = 1.04 - P_up = 100 - surface_conc = siverts_law(T=temperature, S_0=S_0, E_S=E_S, pressure=P_up) - bc_sieverts = dirichletbc( - Constant(my_mesh.mesh, PETSc.ScalarType(surface_conc)), left_dofs, V - ) - bc_outgas = dirichletbc(Constant(my_mesh.mesh, PETSc.ScalarType(0)), right_dofs, V) - my_model.boundary_conditions = [bc_sieverts, bc_outgas] - my_model.create_solver() - - my_model.solver.convergence_criterion = "incremental" - my_model.solver.rtol = 1e-10 - my_model.solver.atol = 1e10 - - my_model.solver.report = True - ksp = my_model.solver.krylov_solver - opts = PETSc.Options() - option_prefix = ksp.getOptionsPrefix() - opts[f"{option_prefix}ksp_type"] = "cg" - opts[f"{option_prefix}pc_type"] = "gamg" - 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), my_mesh.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() - - # analytical solution - S = S_0 * exp(-E_S / F.k_B / float(temperature)) - permeability = float(D) * S - times = np.array(times) - - n_array = np.arange(1, 10000)[:, np.newaxis] - summation = np.sum( - (-1) ** n_array * np.exp(-((np.pi * n_array) ** 2) * float(D) / L**2 * times), - axis=0, - ) - analytical_flux = P_up**0.5 * permeability / L * (2 * summation + 1) - - analytical_flux = np.abs(analytical_flux) - flux_values = np.array(np.abs(flux_values)) - - relative_error = np.abs((flux_values - analytical_flux) / analytical_flux) - - relative_error = relative_error[ - np.where(analytical_flux > 0.01 * np.max(analytical_flux)) - ] - error = relative_error.mean() - - assert error < 0.01 From b0901a0308116bb3c330d9ce1ee5f3401c9dca81 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:59:00 -0400 Subject: [PATCH 6/9] update permeation test --- test/test_permeation_problem.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 8ceddcda5..62b792b30 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -4,7 +4,7 @@ dirichletbc, locate_dofs_topological, ) -from ufl import exp, FacetNormal +from ufl import exp import numpy as np @@ -44,11 +44,10 @@ def siverts_law(T, S_0, E_S, pressure): S = S_0 * exp(-E_S / F.k_B / T) return S * pressure**0.5 - fdim = my_mesh.mesh.topology.dim - 1 left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, fdim, left_facets) + left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, fdim, right_facets) + right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) S_0 = 4.02e21 E_S = 1.04 @@ -143,20 +142,15 @@ def test_permeation_problem_multi_volume(): 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) return S * pressure**0.5 - fdim = my_mesh.mesh.topology.dim - 1 left_facets = my_model.facet_meshtags.find(1) - left_dofs = locate_dofs_topological(V, fdim, left_facets) + left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets) right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, fdim, right_facets) + right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets) S_0 = 4.02e21 E_S = 1.04 From aca8ebb425a49ab558994fb166020e2536d8cd50 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 12:41:15 -0400 Subject: [PATCH 7/9] use mesh property n --- festim/hydrogen_transport_problem.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3b9f2c549..d719cf4c1 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -267,11 +267,10 @@ def run(self, final_time: float): 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)) + surface_flux = form(D * dot(grad(cm), self.mesh.n) * self.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) From b1f529f299b469c3f9ad08344d16e2b7113f2c4a Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 13:02:32 -0400 Subject: [PATCH 8/9] remove conda workflow --- .github/workflows/ci.yml | 66 ++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 33036bf90..682bdeebc 100644 --- a/.github/workflows/ci.yml +++ b/.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 From 6b863f44011be8ced71f341c1f692c584b083d20 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Tue, 17 Oct 2023 15:19:08 -0400 Subject: [PATCH 9/9] Update festim/hydrogen_transport_problem.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- festim/hydrogen_transport_problem.py | 1 - 1 file changed, 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index d719cf4c1..97f547818 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -266,7 +266,6 @@ def run(self, final_time: float): mobile_xdmf.write_function(self.u, t) cm = self.species[0].solution - # TODO this should be a property of Mesh D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature )