From f00fb5708b4fccf8808ccbc80f2046b5ae03ad7d Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 15:50:15 -0400 Subject: [PATCH 01/25] benchmarking script --- .github/workflows/ct.yml | 48 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 .github/workflows/ct.yml diff --git a/.github/workflows/ct.yml b/.github/workflows/ct.yml new file mode 100644 index 000000000..39157d31a --- /dev/null +++ b/.github/workflows/ct.yml @@ -0,0 +1,48 @@ +name: CT +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: Download previous benchmark data + uses: actions/cache@v1 + with: + path: ./cache + key: ${{ runner.os }}-benchmark + + - name: Run benchmark + shell: bash -l {0} + run: | + pytest benchmark.py --benchmark-json output.json + + - name: Store benchmark result + uses: benchmark-action/github-action-benchmark@v1 + with: + tool: 'pytest' + output-file-path: output.json + external-data-json-path: ./cache/benchmark-data.json + fail-on-alert: true + github-token: ${{ secrets.GITHUB_TOKEN }} + comment-on-alert: true From 624375e3afc1aa8e469af041183f63972797cd3e Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 15:55:40 -0400 Subject: [PATCH 02/25] rename all --- .github/workflows/{ct.yml => benchmark.yml} | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename .github/workflows/{ct.yml => benchmark.yml} (97%) diff --git a/.github/workflows/ct.yml b/.github/workflows/benchmark.yml similarity index 97% rename from .github/workflows/ct.yml rename to .github/workflows/benchmark.yml index 39157d31a..db24b9f89 100644 --- a/.github/workflows/ct.yml +++ b/.github/workflows/benchmark.yml @@ -1,8 +1,8 @@ -name: CT +name: benchmark on: [pull_request, push] jobs: - run-on-conda: + benchmark: runs-on: ubuntu-latest steps: From 1383fdd439e1cef74c70fab30666aca53adcfb7f Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:03:16 -0400 Subject: [PATCH 03/25] indentation and add benchmark script --- .github/workflows/benchmark.yml | 8 +- benchmark.py | 293 ++++++++++++++++++++++++++++++++ 2 files changed, 297 insertions(+), 4 deletions(-) create mode 100644 benchmark.py diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index db24b9f89..4a4a10051 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -27,10 +27,10 @@ jobs: pip install .[test] - name: Download previous benchmark data - uses: actions/cache@v1 - with: - path: ./cache - key: ${{ runner.os }}-benchmark + uses: actions/cache@v1 + with: + path: ./cache + key: ${{ runner.os }}-benchmark - name: Run benchmark shell: bash -l {0} diff --git a/benchmark.py b/benchmark.py new file mode 100644 index 000000000..61de1d48a --- /dev/null +++ b/benchmark.py @@ -0,0 +1,293 @@ +import numpy as np +from mpi4py import MPI +from petsc4py import PETSc +from dolfinx.io import XDMFFile +from dolfinx.fem import ( + Constant, + dirichletbc, + Function, + FunctionSpace, + locate_dofs_topological, + locate_dofs_geometrical, + form, + assemble_scalar, +) +from dolfinx.fem.petsc import ( + NonlinearProblem, +) +from dolfinx.nls.petsc import NewtonSolver +from ufl import ( + dot, + FiniteElement, + grad, + TestFunction, + exp, + FacetNormal, + dx, + ds, + Cell, + Mesh, + VectorElement, + Measure, + MixedElement, +) +from dolfinx.mesh import create_mesh, meshtags, locate_entities +from dolfinx import log +import numpy as np +import tqdm.autonotebook +import festim as F +import pytest + + +def pure_fenics(): + # mesh nodes + indices = np.linspace(0, 3e-4, num=5001) + + gdim, shape, degree = 1, "interval", 1 + cell = Cell(shape, geometric_dimension=gdim) + domain = Mesh(VectorElement("Lagrange", cell, degree)) + mesh_points = np.reshape(indices, (len(indices), 1)) + indexes = np.arange(mesh_points.shape[0]) + cells = np.stack((indexes[:-1], indexes[1:]), axis=-1) + my_mesh = create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain) + fdim = my_mesh.topology.dim - 1 + vdim = my_mesh.topology.dim + + elements = FiniteElement("CG", my_mesh.ufl_cell(), 1) + V = FunctionSpace(my_mesh, elements) + u = Function(V) + u_n = Function(V) + v = TestFunction(V) + + borders = [0, 3e-04] + + num_cells = my_mesh.topology.index_map(vdim).size_local + cells = np.arange(num_cells, dtype=np.int32) + markers = np.full(num_cells, 1, dtype=np.int32) + markers[ + locate_entities( + my_mesh, + vdim, + lambda x: np.logical_and(x[0] >= borders[0], x[0] <= borders[1]), + ) + ] = 1 + mesh_tags_volumes = meshtags(my_mesh, vdim, cells, markers) + + tags_volumes = np.array([1], dtype=np.int32) + + dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0)) + dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], indices[-1])) + + dofs_facets = np.array([dofs_L[0], dofs_R[0]], dtype=np.int32) + tags_facets = np.array([1, 2], dtype=np.int32) + + facet_dimension = my_mesh.topology.dim - 1 + volume_dimension = my_mesh.topology.dim + + mesh_tags_facets = meshtags(my_mesh, facet_dimension, dofs_facets, tags_facets) + # mesh_tags_volumes = meshtags(my_mesh, volume_dimension, dofs_volumes, tags_volumes) + ds = Measure("ds", domain=my_mesh, subdomain_data=mesh_tags_facets) + dx = Measure("dx", domain=my_mesh, subdomain_data=mesh_tags_volumes) + + temperature = 500 + k_B = 8.6173303e-5 + n = FacetNormal(my_mesh) + + def siverts_law(T, S_0, E_S, pressure): + S = S_0 * exp(-E_S / k_B / T) + return S * pressure**0.5 + + fdim = my_mesh.topology.dim - 1 + left_facets = mesh_tags_facets.find(1) + left_dofs = locate_dofs_topological(V, fdim, left_facets) + right_facets = mesh_tags_facets.find(2) + right_dofs = locate_dofs_topological(V, fdim, right_facets) + + surface_conc = siverts_law(T=temperature, S_0=4.02e21, E_S=1.04, pressure=100) + bc_sieverts = dirichletbc( + Constant(my_mesh, PETSc.ScalarType(surface_conc)), left_dofs, V + ) + bc_outgas = dirichletbc(Constant(my_mesh, PETSc.ScalarType(0)), right_dofs, V) + bcs = [bc_sieverts, bc_outgas] + + D_0 = 1.9e-7 + E_D = 0.2 + + D = D_0 * exp(-E_D / k_B / temperature) + + dt = 1 / 20 + final_time = 50 + num_steps = int(final_time / dt) + + f = Constant(my_mesh, (PETSc.ScalarType(0))) + F = dot(D * grad(u), grad(v)) * dx(1) + F += ((u - u_n) / dt) * v * dx(1) + + problem = NonlinearProblem(F, u, bcs=bcs) + solver = NewtonSolver(MPI.COMM_WORLD, problem) + solver.convergence_criterion = "incremental" + solver.rtol = 1e-10 + solver.atol = 1e10 + solver.report = True + ksp = 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() + # log.set_log_level(log.LogLevel.INFO) + + mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") + mobile_xdmf.write_mesh(my_mesh) + + flux_values = [] + times = [] + t = 0 + progress = tqdm.autonotebook.tqdm( + desc="Solving H transport problem", total=num_steps + ) + for i in range(num_steps): + progress.update(1) + t += dt + + no_iters, converged = solver.solve(u) + + # surface_flux = form(D * dot(grad(u), n) * ds(2)) + # flux = assemble_scalar(surface_flux) + # flux_values.append(flux) + # times.append(t) + # np.savetxt("outgassing_flux.txt", np.array(flux_values)) + # np.savetxt("times.txt", np.array(times)) + + mobile_xdmf.write_function(u, t) + + u_n.x.array[:] = u.x.array[:] + + mobile_xdmf.close() + + +def festim_script(): + 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 + + # 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) + right_facets = my_model.facet_meshtags.find(2) + right_dofs = locate_dofs_topological(V, 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), 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_fenicsx_benchmark(benchmark): + benchmark(pure_fenics) + + +def test_festim_benchmark(benchmark): + benchmark(festim_script) From 86db081114b1cc26403283aa8ec7052149d48edb Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:16:28 -0400 Subject: [PATCH 04/25] use docker --- .github/workflows/benchmark.yml | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 4a4a10051..e347fd935 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -4,25 +4,13 @@ on: [pull_request, push] jobs: benchmark: runs-on: ubuntu-latest + container: dolfinx/dolfinx:stable 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] @@ -33,7 +21,6 @@ jobs: key: ${{ runner.os }}-benchmark - name: Run benchmark - shell: bash -l {0} run: | pytest benchmark.py --benchmark-json output.json From 202ee8d8363055cd0e25309120e0382b3391a97b Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:20:46 -0400 Subject: [PATCH 05/25] no options --- .github/workflows/benchmark.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index e347fd935..ee680da5a 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -22,7 +22,7 @@ jobs: - name: Run benchmark run: | - pytest benchmark.py --benchmark-json output.json + python3 -m pytest benchmark.py - name: Store benchmark result uses: benchmark-action/github-action-benchmark@v1 From 9acb0057c861cfe02b4193cc2ceecd0dd50818b2 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:27:08 -0400 Subject: [PATCH 06/25] install pytest benchmark --- setup.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.cfg b/setup.cfg index 5d988731b..b24c50ecd 100644 --- a/setup.cfg +++ b/setup.cfg @@ -34,3 +34,4 @@ install_requires = test = pytest >= 5.4.3 pytest-cov + pytest-benchmark From 8688c7bdcde94c5ed8bcfc8b00ef91d97b3493be Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:39:33 -0400 Subject: [PATCH 07/25] no token, use json --- .github/workflows/benchmark.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index ee680da5a..30c90d3a9 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -22,7 +22,7 @@ jobs: - name: Run benchmark run: | - python3 -m pytest benchmark.py + pytest benchmark.py --benchmark-json output.json - name: Store benchmark result uses: benchmark-action/github-action-benchmark@v1 @@ -31,5 +31,4 @@ jobs: output-file-path: output.json external-data-json-path: ./cache/benchmark-data.json fail-on-alert: true - github-token: ${{ secrets.GITHUB_TOKEN }} comment-on-alert: true From 19edfa66b6290c4399e5d95afe580da2bd07f06d Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:44:27 -0400 Subject: [PATCH 08/25] use token --- .github/workflows/benchmark.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 30c90d3a9..e347fd935 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -31,4 +31,5 @@ jobs: output-file-path: output.json external-data-json-path: ./cache/benchmark-data.json fail-on-alert: true + github-token: ${{ secrets.GITHUB_TOKEN }} comment-on-alert: true From 9a66db699061ec1b65af0ec9671ea866558fd847 Mon Sep 17 00:00:00 2001 From: J Dark Date: Mon, 16 Oct 2023 16:59:26 -0400 Subject: [PATCH 09/25] alert on threshold --- .github/workflows/benchmark.yml | 1 + benchmark.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index e347fd935..84fc032d4 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -33,3 +33,4 @@ jobs: fail-on-alert: true github-token: ${{ secrets.GITHUB_TOKEN }} comment-on-alert: true + alert-threshold: "120%" diff --git a/benchmark.py b/benchmark.py index 61de1d48a..97621c64a 100644 --- a/benchmark.py +++ b/benchmark.py @@ -41,7 +41,7 @@ def pure_fenics(): # mesh nodes - indices = np.linspace(0, 3e-4, num=5001) + indices = np.linspace(0, 3e-4, num=6001) gdim, shape, degree = 1, "interval", 1 cell = Cell(shape, geometric_dimension=gdim) From 132bfa2ce799d1a49fa0cdb0ed18cdf0289d3094 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 08:28:47 -0400 Subject: [PATCH 10/25] compare to previous run --- .github/workflows/benchmark.yml | 2 +- benchmark.py | 46 ++++++++------------------------- 2 files changed, 12 insertions(+), 36 deletions(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 84fc032d4..6fab470be 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -22,7 +22,7 @@ jobs: - name: Run benchmark run: | - pytest benchmark.py --benchmark-json output.json + pytest benchmark.py --benchmark-json output.json --benchmark-compare - name: Store benchmark result uses: benchmark-action/github-action-benchmark@v1 diff --git a/benchmark.py b/benchmark.py index 97621c64a..729fe0c33 100644 --- a/benchmark.py +++ b/benchmark.py @@ -41,7 +41,7 @@ def pure_fenics(): # mesh nodes - indices = np.linspace(0, 3e-4, num=6001) + indices = np.linspace(0, 3e-4, num=1001) gdim, shape, degree = 1, "interval", 1 cell = Cell(shape, geometric_dimension=gdim) @@ -153,12 +153,12 @@ def siverts_law(T, S_0, E_S, pressure): no_iters, converged = solver.solve(u) - # surface_flux = form(D * dot(grad(u), n) * ds(2)) - # flux = assemble_scalar(surface_flux) - # flux_values.append(flux) - # times.append(t) - # np.savetxt("outgassing_flux.txt", np.array(flux_values)) - # np.savetxt("times.txt", np.array(times)) + surface_flux = form(D * dot(grad(u), n) * ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(t) + np.savetxt("outgassing_flux.txt", np.array(flux_values)) + np.savetxt("times.txt", np.array(times)) mobile_xdmf.write_function(u, t) @@ -232,8 +232,8 @@ 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) + mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") + mobile_xdmf.write_mesh(my_model.mesh.mesh) final_time = 50 @@ -249,7 +249,7 @@ def siverts_law(T, S_0, E_S, pressure): my_model.solver.solve(u) - # mobile_xdmf.write_function(u, t) + mobile_xdmf.write_function(u, t) surface_flux = form(D * dot(grad(u), n) * my_model.ds(2)) flux = assemble_scalar(surface_flux) @@ -258,31 +258,7 @@ def siverts_law(T, S_0, E_S, pressure): 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 + mobile_xdmf.close() def test_fenicsx_benchmark(benchmark): From d2bea58301a06d606915895f1ce951466f7e8920 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 09:44:09 -0400 Subject: [PATCH 11/25] just test festim, autosave --- .github/workflows/benchmark.yml | 2 +- benchmark.py | 268 +------------------------------- 2 files changed, 3 insertions(+), 267 deletions(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 6fab470be..e3441def5 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -22,7 +22,7 @@ jobs: - name: Run benchmark run: | - pytest benchmark.py --benchmark-json output.json --benchmark-compare + pytest benchmark.py --benchmark-json output.json --benchmark-compare --benchmark-autosave - name: Store benchmark result uses: benchmark-action/github-action-benchmark@v1 diff --git a/benchmark.py b/benchmark.py index 729fe0c33..e0b08a8f4 100644 --- a/benchmark.py +++ b/benchmark.py @@ -1,269 +1,5 @@ -import numpy as np -from mpi4py import MPI -from petsc4py import PETSc -from dolfinx.io import XDMFFile -from dolfinx.fem import ( - Constant, - dirichletbc, - Function, - FunctionSpace, - locate_dofs_topological, - locate_dofs_geometrical, - form, - assemble_scalar, -) -from dolfinx.fem.petsc import ( - NonlinearProblem, -) -from dolfinx.nls.petsc import NewtonSolver -from ufl import ( - dot, - FiniteElement, - grad, - TestFunction, - exp, - FacetNormal, - dx, - ds, - Cell, - Mesh, - VectorElement, - Measure, - MixedElement, -) -from dolfinx.mesh import create_mesh, meshtags, locate_entities -from dolfinx import log -import numpy as np -import tqdm.autonotebook -import festim as F -import pytest - - -def pure_fenics(): - # mesh nodes - indices = np.linspace(0, 3e-4, num=1001) - - gdim, shape, degree = 1, "interval", 1 - cell = Cell(shape, geometric_dimension=gdim) - domain = Mesh(VectorElement("Lagrange", cell, degree)) - mesh_points = np.reshape(indices, (len(indices), 1)) - indexes = np.arange(mesh_points.shape[0]) - cells = np.stack((indexes[:-1], indexes[1:]), axis=-1) - my_mesh = create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain) - fdim = my_mesh.topology.dim - 1 - vdim = my_mesh.topology.dim - - elements = FiniteElement("CG", my_mesh.ufl_cell(), 1) - V = FunctionSpace(my_mesh, elements) - u = Function(V) - u_n = Function(V) - v = TestFunction(V) - - borders = [0, 3e-04] - - num_cells = my_mesh.topology.index_map(vdim).size_local - cells = np.arange(num_cells, dtype=np.int32) - markers = np.full(num_cells, 1, dtype=np.int32) - markers[ - locate_entities( - my_mesh, - vdim, - lambda x: np.logical_and(x[0] >= borders[0], x[0] <= borders[1]), - ) - ] = 1 - mesh_tags_volumes = meshtags(my_mesh, vdim, cells, markers) - - tags_volumes = np.array([1], dtype=np.int32) - - dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0)) - dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], indices[-1])) - - dofs_facets = np.array([dofs_L[0], dofs_R[0]], dtype=np.int32) - tags_facets = np.array([1, 2], dtype=np.int32) - - facet_dimension = my_mesh.topology.dim - 1 - volume_dimension = my_mesh.topology.dim - - mesh_tags_facets = meshtags(my_mesh, facet_dimension, dofs_facets, tags_facets) - # mesh_tags_volumes = meshtags(my_mesh, volume_dimension, dofs_volumes, tags_volumes) - ds = Measure("ds", domain=my_mesh, subdomain_data=mesh_tags_facets) - dx = Measure("dx", domain=my_mesh, subdomain_data=mesh_tags_volumes) - - temperature = 500 - k_B = 8.6173303e-5 - n = FacetNormal(my_mesh) - - def siverts_law(T, S_0, E_S, pressure): - S = S_0 * exp(-E_S / k_B / T) - return S * pressure**0.5 - - fdim = my_mesh.topology.dim - 1 - left_facets = mesh_tags_facets.find(1) - left_dofs = locate_dofs_topological(V, fdim, left_facets) - right_facets = mesh_tags_facets.find(2) - right_dofs = locate_dofs_topological(V, fdim, right_facets) - - surface_conc = siverts_law(T=temperature, S_0=4.02e21, E_S=1.04, pressure=100) - bc_sieverts = dirichletbc( - Constant(my_mesh, PETSc.ScalarType(surface_conc)), left_dofs, V - ) - bc_outgas = dirichletbc(Constant(my_mesh, PETSc.ScalarType(0)), right_dofs, V) - bcs = [bc_sieverts, bc_outgas] - - D_0 = 1.9e-7 - E_D = 0.2 - - D = D_0 * exp(-E_D / k_B / temperature) - - dt = 1 / 20 - final_time = 50 - num_steps = int(final_time / dt) - - f = Constant(my_mesh, (PETSc.ScalarType(0))) - F = dot(D * grad(u), grad(v)) * dx(1) - F += ((u - u_n) / dt) * v * dx(1) - - problem = NonlinearProblem(F, u, bcs=bcs) - solver = NewtonSolver(MPI.COMM_WORLD, problem) - solver.convergence_criterion = "incremental" - solver.rtol = 1e-10 - solver.atol = 1e10 - solver.report = True - ksp = 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() - # log.set_log_level(log.LogLevel.INFO) - - mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") - mobile_xdmf.write_mesh(my_mesh) - - flux_values = [] - times = [] - t = 0 - progress = tqdm.autonotebook.tqdm( - desc="Solving H transport problem", total=num_steps - ) - for i in range(num_steps): - progress.update(1) - t += dt - - no_iters, converged = solver.solve(u) - - surface_flux = form(D * dot(grad(u), n) * ds(2)) - flux = assemble_scalar(surface_flux) - flux_values.append(flux) - times.append(t) - np.savetxt("outgassing_flux.txt", np.array(flux_values)) - np.savetxt("times.txt", np.array(times)) - - mobile_xdmf.write_function(u, t) - - u_n.x.array[:] = u.x.array[:] - - mobile_xdmf.close() - - -def festim_script(): - 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 - - # 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) - right_facets = my_model.facet_meshtags.find(2) - right_dofs = locate_dofs_topological(V, 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), 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() - - -def test_fenicsx_benchmark(benchmark): - benchmark(pure_fenics) +from test.test_permeation_problem import test_permeation_problem def test_festim_benchmark(benchmark): - benchmark(festim_script) + benchmark(test_permeation_problem) From dfe39efb761ba1c5a9883450bfe177b285ef91e2 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 09:49:10 -0400 Subject: [PATCH 12/25] test worse performance --- test/test_permeation_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index b8534107e..a4824e4bf 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -18,7 +18,7 @@ def test_permeation_problem(): L = 3e-04 - vertices = np.linspace(0, L, num=1001) + vertices = np.linspace(0, L, num=10001) my_mesh = F.Mesh1D(vertices) From 6ef1b4224bd463f309786cdaaf32bf3832c7ff14 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 09:58:17 -0400 Subject: [PATCH 13/25] permission to write benchmark results in cache --- .github/workflows/benchmark.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index e3441def5..305ad8cd8 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -4,6 +4,7 @@ on: [pull_request, push] jobs: benchmark: runs-on: ubuntu-latest + permissions: write-all container: dolfinx/dolfinx:stable steps: From b7d956c0ea3062cc1534b458f42d440ea4d3eba5 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:35:45 -0400 Subject: [PATCH 14/25] simplify benchmark --- .github/workflows/benchmark.yml | 20 +--- benchmark.py | 202 +++++++++++++++++++++++++++++++- test/test_permeation_problem.py | 2 +- 3 files changed, 202 insertions(+), 22 deletions(-) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 305ad8cd8..5ad6ed1c0 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -4,7 +4,6 @@ on: [pull_request, push] jobs: benchmark: runs-on: ubuntu-latest - permissions: write-all container: dolfinx/dolfinx:stable steps: @@ -15,23 +14,6 @@ jobs: run: | pip install .[test] - - name: Download previous benchmark data - uses: actions/cache@v1 - with: - path: ./cache - key: ${{ runner.os }}-benchmark - - name: Run benchmark run: | - pytest benchmark.py --benchmark-json output.json --benchmark-compare --benchmark-autosave - - - name: Store benchmark result - uses: benchmark-action/github-action-benchmark@v1 - with: - tool: 'pytest' - output-file-path: output.json - external-data-json-path: ./cache/benchmark-data.json - fail-on-alert: true - github-token: ${{ secrets.GITHUB_TOKEN }} - comment-on-alert: true - alert-threshold: "120%" + python3 benchmark.py diff --git a/benchmark.py b/benchmark.py index e0b08a8f4..12c7f676c 100644 --- a/benchmark.py +++ b/benchmark.py @@ -1,5 +1,203 @@ +from mpi4py import MPI +from petsc4py import PETSc +from dolfinx.io import XDMFFile +from dolfinx.fem import ( + Constant, + dirichletbc, + Function, + FunctionSpace, + locate_dofs_topological, + locate_dofs_geometrical, + form, + assemble_scalar, +) +from dolfinx.fem.petsc import ( + NonlinearProblem, +) +from dolfinx.nls.petsc import NewtonSolver +from ufl import ( + dot, + FiniteElement, + grad, + TestFunction, + exp, + FacetNormal, + dx, + ds, + Cell, + Mesh, + VectorElement, + Measure, +) +from dolfinx.mesh import create_mesh, meshtags, locate_entities +import numpy as np +import tqdm.autonotebook +import time from test.test_permeation_problem import test_permeation_problem -def test_festim_benchmark(benchmark): - benchmark(test_permeation_problem) +def fenics_test_permeation_problem(): + L = 3e-04 + indices = np.linspace(0, L, num=1001) + gdim, shape, degree = 1, "interval", 1 + cell = Cell(shape, geometric_dimension=gdim) + domain = Mesh(VectorElement("Lagrange", cell, degree)) + mesh_points = np.reshape(indices, (len(indices), 1)) + indexes = np.arange(mesh_points.shape[0]) + cells = np.stack((indexes[:-1], indexes[1:]), axis=-1) + my_mesh = create_mesh(MPI.COMM_WORLD, cells, mesh_points, domain) + fdim = my_mesh.topology.dim - 1 + vdim = my_mesh.topology.dim + n = FacetNormal(my_mesh) + + elements = FiniteElement("CG", my_mesh.ufl_cell(), 1) + V = FunctionSpace(my_mesh, elements) + u = Function(V) + u_n = Function(V) + v = TestFunction(V) + + num_cells = my_mesh.topology.index_map(vdim).size_local + cells = np.arange(num_cells, dtype=np.int32) + markers = np.full(num_cells, 1, dtype=np.int32) + markers[ + locate_entities( + my_mesh, + vdim, + lambda x: np.logical_and(x[0] >= indices[0], x[0] <= indices[1]), + ) + ] = 1 + mesh_tags_volumes = meshtags(my_mesh, vdim, cells, markers) + + dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0)) + dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], indices[-1])) + + dofs_facets = np.array([dofs_L[0], dofs_R[0]], dtype=np.int32) + tags_facets = np.array([1, 2], dtype=np.int32) + + mesh_tags_facets = meshtags(my_mesh, fdim, dofs_facets, tags_facets) + mesh_tags_volumes = meshtags(my_mesh, vdim, cells, markers) + ds = Measure("ds", domain=my_mesh, subdomain_data=mesh_tags_facets) + dx = Measure("dx", domain=my_mesh, subdomain_data=mesh_tags_volumes) + + temperature = 500 + k_B = 8.6173303e-5 + + def siverts_law(T, S_0, E_S, pressure): + S = S_0 * exp(-E_S / k_B / T) + return S * pressure**0.5 + + left_facets = mesh_tags_facets.find(1) + left_dofs = locate_dofs_topological(V, fdim, left_facets) + right_facets = mesh_tags_facets.find(2) + right_dofs = locate_dofs_topological(V, 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, PETSc.ScalarType(surface_conc)), left_dofs, V + ) + bc_outgas = dirichletbc(Constant(my_mesh, PETSc.ScalarType(0)), right_dofs, V) + bcs = [bc_sieverts, bc_outgas] + + D_0 = 1.9e-7 + E_D = 0.2 + D = D_0 * exp(-E_D / k_B / temperature) + + dt = 1 / 20 + final_time = 50 + num_steps = int(final_time / dt) + + F = dot(D * grad(u), grad(v)) * dx(1) + F += ((u - u_n) / dt) * v * dx(1) + + problem = NonlinearProblem(F, u, bcs=bcs) + solver = NewtonSolver(MPI.COMM_WORLD, problem) + solver.convergence_criterion = "incremental" + solver.rtol = 1e-10 + solver.atol = 1e10 + solver.report = True + ksp = 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_mesh) + + flux_values = [] + times = [] + t = 0 + progress = tqdm.autonotebook.tqdm( + desc="Solving H transport problem", total=num_steps + ) + for i in range(num_steps): + progress.update(1) + t += dt + + solver.solve(u) + + surface_flux = form(D * dot(grad(u), n) * ds(2)) + flux = assemble_scalar(surface_flux) + flux_values.append(flux) + times.append(t) + + mobile_xdmf.write_function(u, t) + + u_n.x.array[:] = u.x.array[:] + + mobile_xdmf.close() + + # analytical solution + S = S_0 * exp(-E_S / 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_festim_vs_fenics_permeation_benchmark(): + repetitions = 10 + + fenics_times = [] + for i in range(repetitions): + start = time.time() + fenics_test_permeation_problem() + fenics_times.append(time.time() - start) + fenics_time = np.mean(fenics_times) + + festim_times = [] + for i in range(repetitions): + start = time.time() + test_permeation_problem() + festim_times.append(time.time() - start) + festim_time = np.mean(festim_times) + + diff = (np.abs(fenics_time - festim_time) / ((fenics_time + festim_time) / 2)) * 100 + if diff > 20: + raise ValueError(f"festim is {diff:.1f}% slower than fenics") + + +if __name__ == "__main__": + test_festim_vs_fenics_permeation_benchmark() diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index a4824e4bf..b8534107e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -18,7 +18,7 @@ def test_permeation_problem(): L = 3e-04 - vertices = np.linspace(0, L, num=10001) + vertices = np.linspace(0, L, num=1001) my_mesh = F.Mesh1D(vertices) From de817ea415650c37602df59a01cd6f78a34cd277 Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:39:39 -0400 Subject: [PATCH 15/25] test benchmark fails --- test/test_permeation_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index b8534107e..88fad97d2 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -18,7 +18,7 @@ def test_permeation_problem(): L = 3e-04 - vertices = np.linspace(0, L, num=1001) + vertices = np.linspace(0, L, num=3001) my_mesh = F.Mesh1D(vertices) From b9947b6fb5cd7c139fa1a70d2a8c95f4c2c9ebad Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:48:09 -0400 Subject: [PATCH 16/25] back to normal --- test/test_permeation_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 88fad97d2..b8534107e 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -18,7 +18,7 @@ def test_permeation_problem(): L = 3e-04 - vertices = np.linspace(0, L, num=3001) + vertices = np.linspace(0, L, num=1001) my_mesh = F.Mesh1D(vertices) From 819b138859f814ffbfe07266e4b0113789b2e26c Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 11:59:40 -0400 Subject: [PATCH 17/25] dont install benchmark --- setup.cfg | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index b24c50ecd..5d988731b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -34,4 +34,3 @@ install_requires = test = pytest >= 5.4.3 pytest-cov - pytest-benchmark From 206558b6edd6233e12f97c1a4b1b0b6bba87ec6f Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 16:37:15 -0400 Subject: [PATCH 18/25] move benchmark to test folder --- .github/workflows/benchmark.yml | 2 +- benchmark.py => test/benchmark.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) rename benchmark.py => test/benchmark.py (98%) diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 5ad6ed1c0..66c6ab2f4 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -16,4 +16,4 @@ jobs: - name: Run benchmark run: | - python3 benchmark.py + python3 test/benchmark.py diff --git a/benchmark.py b/test/benchmark.py similarity index 98% rename from benchmark.py rename to test/benchmark.py index 12c7f676c..0c8894813 100644 --- a/benchmark.py +++ b/test/benchmark.py @@ -33,7 +33,7 @@ import numpy as np import tqdm.autonotebook import time -from test.test_permeation_problem import test_permeation_problem +from test_permeation_problem import test_permeation_problem def fenics_test_permeation_problem(): From 505fc374cad112c017d0f7a77aaccdda43d992ad Mon Sep 17 00:00:00 2001 From: J Dark Date: Tue, 17 Oct 2023 16:54:10 -0400 Subject: [PATCH 19/25] add doc strings, print difference --- test/benchmark.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/benchmark.py b/test/benchmark.py index 0c8894813..d81cec771 100644 --- a/test/benchmark.py +++ b/test/benchmark.py @@ -178,6 +178,8 @@ def siverts_law(T, S_0, E_S, pressure): def test_festim_vs_fenics_permeation_benchmark(): + """Runs a problem with pure fenicsx and the same problem with FESTIM and + raise ValueError if difference is too high""" repetitions = 10 fenics_times = [] @@ -196,7 +198,11 @@ def test_festim_vs_fenics_permeation_benchmark(): diff = (np.abs(fenics_time - festim_time) / ((fenics_time + festim_time) / 2)) * 100 if diff > 20: - raise ValueError(f"festim is {diff:.1f}% slower than fenics") + raise ValueError( + f"festim is {diff:.1f}% slower than fenics, current acceptble threshold of 20%" + ) + else: + print(f"festim is {diff:.1f}% slower than fenics") if __name__ == "__main__": From 39e0d2d6d004be13c202e64ea8cb75f8d14610a5 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 08:10:54 -0400 Subject: [PATCH 20/25] function takes arguments mesh size and time --- test/test_permeation_problem.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index 8ceddcda5..c1b003a39 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -11,9 +11,9 @@ import festim as F -def test_permeation_problem(): +def test_permeation_problem(final_time=50, mesh_size=1001): L = 3e-04 - vertices = np.linspace(0, L, num=1001) + vertices = np.linspace(0, L, num=mesh_size) my_mesh = F.Mesh1D(vertices) @@ -32,8 +32,6 @@ 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) @@ -74,8 +72,6 @@ def siverts_law(T, S_0, E_S, pressure): opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" ksp.setFromOptions() - final_time = 50 - times, flux_values = my_model.run(final_time=final_time) # analytical solution From 9887a96492c7380e3621d01829f6918ed865967e Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 08:11:11 -0400 Subject: [PATCH 21/25] shouldnt be done in the loop --- festim/hydrogen_transport_problem.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 3b9f2c549..c426aa288 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -253,7 +253,11 @@ def run(self, final_time: float): mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") mobile_xdmf.write_mesh(self.mesh.mesh) - + n = FacetNormal(self.mesh.mesh) + D = self.subdomains[0].material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature + ) + cm = self.species[0].solution progress = tqdm.autonotebook.tqdm( desc="Solving H transport problem", total=final_time ) @@ -265,12 +269,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 - 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) From 7cbb8245c838d4d0e1639b7a920bf359effa18e0 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 09:59:31 -0400 Subject: [PATCH 22/25] benchmark to use 20001 cells --- test/benchmark.py | 49 ++++++++++++++++++++--------------------------- 1 file changed, 21 insertions(+), 28 deletions(-) diff --git a/test/benchmark.py b/test/benchmark.py index d81cec771..5d07beb0b 100644 --- a/test/benchmark.py +++ b/test/benchmark.py @@ -36,9 +36,9 @@ from test_permeation_problem import test_permeation_problem -def fenics_test_permeation_problem(): +def fenics_test_permeation_problem(mesh_size=1001): L = 3e-04 - indices = np.linspace(0, L, num=1001) + indices = np.linspace(0, L, num=mesh_size) gdim, shape, degree = 1, "interval", 1 cell = Cell(shape, geometric_dimension=gdim) domain = Mesh(VectorElement("Lagrange", cell, degree)) @@ -133,21 +133,21 @@ def siverts_law(T, S_0, E_S, pressure): times = [] t = 0 progress = tqdm.autonotebook.tqdm( - desc="Solving H transport problem", total=num_steps + desc="Solving H transport problem", total=final_time ) - for i in range(num_steps): - progress.update(1) - t += dt + while t < final_time: + progress.update(float(dt)) + t += float(dt) solver.solve(u) + mobile_xdmf.write_function(u, t) + surface_flux = form(D * dot(grad(u), n) * ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) times.append(t) - mobile_xdmf.write_function(u, t) - u_n.x.array[:] = u.x.array[:] mobile_xdmf.close() @@ -180,29 +180,22 @@ def siverts_law(T, S_0, E_S, pressure): def test_festim_vs_fenics_permeation_benchmark(): """Runs a problem with pure fenicsx and the same problem with FESTIM and raise ValueError if difference is too high""" - repetitions = 10 - - fenics_times = [] - for i in range(repetitions): - start = time.time() - fenics_test_permeation_problem() - fenics_times.append(time.time() - start) - fenics_time = np.mean(fenics_times) - - festim_times = [] - for i in range(repetitions): - start = time.time() - test_permeation_problem() - festim_times.append(time.time() - start) - festim_time = np.mean(festim_times) - - diff = (np.abs(fenics_time - festim_time) / ((fenics_time + festim_time) / 2)) * 100 - if diff > 20: + + start = time.time() + fenics_test_permeation_problem(mesh_size=20001) + fenics_time = time.time() - start + + start = time.time() + test_permeation_problem(mesh_size=20001) + festim_time = time.time() - start + + diff = (fenics_time - festim_time) / fenics_time + if diff < -0.1: raise ValueError( - f"festim is {diff:.1f}% slower than fenics, current acceptble threshold of 20%" + f"festim is {np.abs(diff):.1%} slower than fenics, current acceptable threshold of 10%" ) else: - print(f"festim is {diff:.1f}% slower than fenics") + print(f"avg relative diff between festim and fenics {diff:.1%}") if __name__ == "__main__": From cb045ccce9cc33e149439a7701173f10172b1115 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 09:59:43 -0400 Subject: [PATCH 23/25] dont need total time for now --- test/test_permeation_problem.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py index c1b003a39..aad1821a0 100644 --- a/test/test_permeation_problem.py +++ b/test/test_permeation_problem.py @@ -11,7 +11,7 @@ import festim as F -def test_permeation_problem(final_time=50, mesh_size=1001): +def test_permeation_problem(mesh_size=1001): L = 3e-04 vertices = np.linspace(0, L, num=mesh_size) @@ -72,6 +72,8 @@ def siverts_law(T, S_0, E_S, pressure): opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps" ksp.setFromOptions() + final_time = 50 + times, flux_values = my_model.run(final_time=final_time) # analytical solution From 62353ce2bdd0049610acab2430918ddd02f8e666 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 10:22:37 -0400 Subject: [PATCH 24/25] fix coverage --- festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 81bfffd09..3fc9311a0 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -253,7 +253,7 @@ def run(self, final_time: float): mobile_xdmf = XDMFFile(MPI.COMM_WORLD, "mobile_concentration.xdmf", "w") mobile_xdmf.write_mesh(self.mesh.mesh) - n = FacetNormal(self.mesh.mesh) + n = self.mesh.n D = self.subdomains[0].material.get_diffusion_coefficient( self.mesh.mesh, self.temperature ) From b9aab59d502bb79d49b0cbd493f1fda3edbe28a3 Mon Sep 17 00:00:00 2001 From: J Dark Date: Wed, 18 Oct 2023 10:41:15 -0400 Subject: [PATCH 25/25] refactoring --- test/benchmark.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/benchmark.py b/test/benchmark.py index 5d07beb0b..e69f2a568 100644 --- a/test/benchmark.py +++ b/test/benchmark.py @@ -190,9 +190,10 @@ def test_festim_vs_fenics_permeation_benchmark(): festim_time = time.time() - start diff = (fenics_time - festim_time) / fenics_time - if diff < -0.1: + threshold = -0.1 + if diff < threshold: raise ValueError( - f"festim is {np.abs(diff):.1%} slower than fenics, current acceptable threshold of 10%" + f"festim is {np.abs(diff):.1%} slower than fenics, current acceptable threshold of {np.abs(threshold):.1%}" ) else: print(f"avg relative diff between festim and fenics {diff:.1%}")