diff --git a/.github/workflows/run-tutorials.yml b/.github/workflows/run-tutorials.yml index 1c514a1..1cff7ce 100644 --- a/.github/workflows/run-tutorials.yml +++ b/.github/workflows/run-tutorials.yml @@ -12,7 +12,7 @@ jobs: run_ht_simple: name: Run HT, simple runs-on: ubuntu-latest - container: benjaminrodenberg/precice:2.5.0 + container: precice/precice:2.5.1 steps: - name: Checkout Repository uses: actions/checkout@v2 @@ -27,8 +27,9 @@ jobs: run: python3 -c "import dolfinx; print('FEniCS-X version '+dolfinx.__version__)" - name: Install adapter run: pip3 install --user . - - name: Get tutorials - run: git clone -b master https://github.com/PhilipHildebrand/tutorials.git + # Tutorial is currently hosted under https://github.com/precice/fenicsx-adapter, should be moved to https://github.com/precice/tutorials when working. See https://github.com/precice/tutorials/issues/491 + # - name: Get tutorials + # run: git clone -b develop https://github.com/precice/tutorials.git - name: Run tutorial run: | cd tutorials/partitioned-heat-conduction/fenicsx diff --git a/README.md b/README.md index 89a2695..8a1f5b9 100644 --- a/README.md +++ b/README.md @@ -82,7 +82,7 @@ If you are using FEniCSx, please also consider the information on [the official ## Development history -2018: The initial version of the [fenics-adapter](https://github.com/precice/fenics-adapter) was developed by [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) during his research stay at Lund University in the group for [Numerical Analysis](http://www.maths.lu.se/english/research/research-divisions/numerical-analysis/) in close collaboration with [Peter Meisrimel](https://portal.research.lu.se/en/persons/peter-meisrimel). +2018: The initial version of the [fenics-adapter](https://github.com/precice/fenics-adapter) was developed by [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) during his research stay at Lund University in the group for [Numerical Analysis](http://www.maths.lu.se/english/research/research-divisions/numerical-analysis/) in close collaboration with Peter Meisrimel. 2019: [Richard Hertrich](https://github.com/richahert) contributed the possibility to perform FSI simulations using the adapter in his [Bachelor thesis](https://mediatum.ub.tum.de/node?id=1520579). @@ -90,4 +90,4 @@ If you are using FEniCSx, please also consider the information on [the official 2021: For development of FEniCSx support, `precice/fenics-adapter@v1.2.0` was forked as `precice/fenicsx-adapter`. The required modifications were carried out by [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). -2023: [Philip Hildebrand](https://github.com/PhilipHildebrand) updated the adapter to a first minimal working version (https://github.com/precice/fenicsx-adapter/pull/15) and contributed a first tutorial (https://github.com/precice/tutorials/pull/317) in the scope of his Bachelor's thesis ["Extending the FEniCSx Adapter for the Coupling Library preCICE"](https://mediatum.ub.tum.de/node?id=1706280) under supervision of [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). +2023: [Philip Hildebrand](https://github.com/PhilipHildebrand) updated the adapter to a [first minimal working version](https://github.com/precice/fenicsx-adapter/pull/15) and contributed a [first tutorial](https://github.com/precice/tutorials/pull/317) in the scope of his Bachelor's thesis ["Extending the FEniCSx Adapter for the Coupling Library preCICE"](https://mediatum.ub.tum.de/node?id=1706280) under supervision of [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). diff --git a/setup.py b/setup.py index b063af3..0a24b23 100644 --- a/setup.py +++ b/setup.py @@ -39,7 +39,7 @@ author_email='info@precice.org', license='LGPL-3.0', packages=['fenicsxprecice'], - install_requires=['pyprecice>=2.0.0', 'scipy', 'numpy>=1.13.3', 'mpi4py'], + install_requires=['pyprecice~=2.0.0', 'scipy', 'numpy>=1.13.3', 'mpi4py'], tests_require=['sympy'], test_suite='tests', zip_safe=False) diff --git a/tutorials/partitioned-heat-conduction/fenicsx/clean.sh b/tutorials/partitioned-heat-conduction/fenicsx/clean.sh new file mode 100755 index 0000000..3a8b461 --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_fenics . diff --git a/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py b/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py new file mode 100644 index 0000000..069dd87 --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py @@ -0,0 +1,16 @@ +from ufl import dx +from dolfinx.fem import assemble_scalar, form +import numpy as np +from mpi4py import MPI + + +def compute_errors(u_approx, u_ref, total_error_tol=10 ** -4): + mesh = u_ref.function_space.mesh + + # compute total L2 error between reference and calculated solution + error_pointwise = form(((u_approx - u_ref) / u_ref) ** 2 * dx) + error_total = np.sqrt(mesh.comm.allreduce(assemble_scalar(error_pointwise), MPI.SUM)) + + assert (error_total < total_error_tol) + + return error_total diff --git a/tutorials/partitioned-heat-conduction/fenicsx/heat.py b/tutorials/partitioned-heat-conduction/fenicsx/heat.py new file mode 100644 index 0000000..7ecef9f --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/heat.py @@ -0,0 +1,253 @@ +""" +The basic example is taken from "Langtangen, Hans Petter, and Anders Logg. Solving PDEs in Python: The FEniCS +Tutorial I. Springer International Publishing, 2016." + +The example code has been extended with preCICE API calls and mixed boundary conditions to allow for a Dirichlet-Neumann +coupling of two separate heat equations. It also has been adapted to be compatible with FEniCSx. + +The original source code can be found on https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html. + +Heat equation with Dirichlet conditions. (Dirichlet problem) + u'= Laplace(u) + f in the unit square [0,1] x [0,1] + u = u_C on the coupling boundary at x = 1 + u = u_D on the remaining boundary + u = u_0 at t = 0 + u = 1 + x^2 + alpha*y^2 + \beta*t + f = beta - 2 - 2*alpha + +Heat equation with mixed boundary conditions. (Neumann problem) + u'= Laplace(u) + f in the shifted unit square [1,2] x [0,1] + du/dn = f_N on the coupling boundary at x = 1 + u = u_D on the remaining boundary + u = u_0 at t = 0 + u = 1 + x^2 + alpha*y^2 + \beta*t + f = beta - 2 - 2*alpha +""" + +from __future__ import print_function, division +from mpi4py import MPI +from dolfinx.fem import Function, FunctionSpace, Expression, Constant, dirichletbc, locate_dofs_geometrical +from dolfinx.fem.petsc import LinearProblem +from dolfinx.io import XDMFFile +from ufl import TrialFunction, TestFunction, dx, ds, dot, grad, inner, lhs, rhs, FiniteElement, VectorElement +from fenicsxprecice import Adapter +from errorcomputation import compute_errors +from my_enums import ProblemType, DomainPart +import argparse +import numpy as np +from problem_setup import get_geometry + + +def determine_gradient(V_g, u): + """ + compute flux following http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html#tut-poisson-gradu + :param V_g: Vector function space + :param u: solution where gradient is to be determined + """ + + w = TrialFunction(V_g) + v = TestFunction(V_g) + + a = inner(w, v) * dx + L = inner(grad(u), v) * dx + problem = LinearProblem(a, L) + return problem.solve() + + +parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") +command_group = parser.add_mutually_exclusive_group(required=True) +command_group.add_argument("-d", "--dirichlet", help="create a dirichlet problem", dest="dirichlet", + action="store_true") +command_group.add_argument("-n", "--neumann", help="create a neumann problem", dest="neumann", action="store_true") +parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-6,) + +args = parser.parse_args() + +fenics_dt = .1 # time step size +# Error is bounded by coupling accuracy. In theory we would obtain the analytical solution. +error_tol = args.error_tol + +alpha = 3 # parameter alpha +beta = 1.2 # parameter beta + +if args.dirichlet and not args.neumann: + problem = ProblemType.DIRICHLET + domain_part = DomainPart.LEFT +elif args.neumann and not args.dirichlet: + problem = ProblemType.NEUMANN + domain_part = DomainPart.RIGHT + +mesh, coupling_boundary, remaining_boundary = get_geometry(MPI.COMM_WORLD, domain_part) + +# Define function space using mesh +scalar_element = FiniteElement("P", mesh.ufl_cell(), 2) +vector_element = VectorElement("P", mesh.ufl_cell(), 1) +V = FunctionSpace(mesh, scalar_element) +V_g = FunctionSpace(mesh, vector_element) +W = V_g.sub(0).collapse()[0] + +# Define boundary conditions + + +class Expression_u_D: + def __init__(self): + self.t = 0.0 + + def eval(self, x): + return np.full(x.shape[1], 1 + x[0] * x[0] + alpha * x[1] * x[1] + beta * self.t) + + +u_D = Expression_u_D() +u_D_function = Function(V) +u_D_function.interpolate(u_D.eval) + +if problem is ProblemType.DIRICHLET: + # Define flux in x direction + class Expression_f_N: + def __init__(self): + pass + + def eval(self, x): + return np.full(x.shape[1], 2 * x[0]) + + f_N = Expression_f_N() + f_N_function = Function(V) + f_N_function.interpolate(f_N.eval) + +# Define initial value +u_n = Function(V, name="Temperature") +u_n.interpolate(u_D.eval) + +precice, precice_dt, initial_data = None, 0.0, None + +# Initialize the adapter according to the specific participant +if problem is ProblemType.DIRICHLET: + precice = Adapter(MPI.COMM_WORLD, adapter_config_filename="precice-adapter-config-D.json") + precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function) +elif problem is ProblemType.NEUMANN: + precice = Adapter(MPI.COMM_WORLD, adapter_config_filename="precice-adapter-config-N.json") + precice_dt = precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D_function) + +dt = Constant(mesh, 0.0) +dt.value = np.min([fenics_dt, precice_dt]) + +# Define variational problem +u = TrialFunction(V) +v = TestFunction(V) + + +class Expression_f: + def __init__(self): + self.t = 0.0 + + def eval(self, x): + return np.full(x.shape[1], beta - 2 - 2 * alpha) + + +f = Expression_f() +f_function = Function(V) +F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f_function) * v * dx + +dofs_remaining = locate_dofs_geometrical(V, remaining_boundary) +bcs = [dirichletbc(u_D_function, dofs_remaining)] + +coupling_expression = precice.create_coupling_expression() +read_data = precice.read_data() +precice.update_coupling_expression(coupling_expression, read_data) +if problem is ProblemType.DIRICHLET: + # modify Dirichlet boundary condition on coupling interface + dofs_coupling = locate_dofs_geometrical(V, coupling_boundary) + bcs.append(dirichletbc(coupling_expression, dofs_coupling)) +if problem is ProblemType.NEUMANN: + # modify Neumann boundary condition on coupling interface, modify weak + # form correspondingly + F += v * coupling_expression * ds + +a, L = lhs(F), rhs(F) + +# Time-stepping +u_np1 = Function(V, name="Temperature") +t = 0 + +# reference solution at t=0 +u_ref = Function(V, name="reference") +u_ref.interpolate(u_D_function) + +# Generating output files +temperature_out = XDMFFile(MPI.COMM_WORLD, f"./output/{precice.get_participant_name()}.xdmf", "w") +temperature_out.write_mesh(mesh) +ref_out = XDMFFile(MPI.COMM_WORLD, f"./output/ref{precice.get_participant_name()}.xdmf", "w") +ref_out.write_mesh(mesh) + +# output solution at t=0, n=0 +n = 0 + +temperature_out.write_function(u_n, t) +ref_out.write_function(u_ref, t) + +u_D.t = t + dt.value +u_D_function.interpolate(u_D.eval) +f.t = t + dt.value +f_function.interpolate(f.eval) + +if problem is ProblemType.DIRICHLET: + flux = Function(V_g, name="Heat-Flux") + +while precice.is_coupling_ongoing(): + + # write checkpoint + if precice.is_action_required(precice.action_write_iteration_checkpoint()): + precice.store_checkpoint(u_n, t, n) + + read_data = precice.read_data() + + # Update the coupling expression with the new read data + precice.update_coupling_expression(coupling_expression, read_data) + + dt.value = np.min([fenics_dt, precice_dt]) + + linear_problem = LinearProblem(a, L, bcs=bcs) + u_np1 = linear_problem.solve() + + # Write data to preCICE according to which problem is being solved + if problem is ProblemType.DIRICHLET: + # Dirichlet problem reads temperature and writes flux on boundary to Neumann problem + flux = determine_gradient(V_g, u_np1) + flux_x = Function(W) + flux_x.interpolate(flux.sub(0)) + precice.write_data(flux_x) + elif problem is ProblemType.NEUMANN: + # Neumann problem reads flux and writes temperature on boundary to Dirichlet problem + precice.write_data(u_np1) + + precice_dt = precice.advance(dt.value) + + # roll back to checkpoint + if precice.is_action_required(precice.action_read_iteration_checkpoint()): + u_cp, t_cp, n_cp = precice.retrieve_checkpoint() + u_n.interpolate(u_cp) + t = t_cp + n = n_cp + else: # update solution + u_n.interpolate(u_np1) + t += dt.value + n += 1 + + if precice.is_time_window_complete(): + u_ref.interpolate(u_D_function) + error = compute_errors(u_n, u_ref, total_error_tol=error_tol) + print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error)) + print('output u^%d and u_ref^%d' % (n, n)) + + temperature_out.write_function(u_n, t) + ref_out.write_function(u_ref, t) + + # Update Dirichlet BC + u_D.t = t + dt.value + u_D_function.interpolate(u_D.eval) + f.t = t + dt.value + f_function.interpolate(f.eval) + + +# Hold plot +precice.finalize() diff --git a/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py b/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py new file mode 100644 index 0000000..5708bad --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py @@ -0,0 +1,17 @@ +from enum import Enum + + +class ProblemType(Enum): + """ + Enum defines problem type. Details see above. + """ + DIRICHLET = 1 # Dirichlet problem + NEUMANN = 2 # Neumann problem + + +class DomainPart(Enum): + """ + Enum defines which part of the domain [x_left, x_right] x [y_bottom, y_top] we compute. + """ + LEFT = 1 # left part of domain in simple interface case + RIGHT = 2 # right part of domain in simple interface case diff --git a/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json b/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json new file mode 100644 index 0000000..7989ba3 --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json @@ -0,0 +1,9 @@ +{ + "participant_name": "Dirichlet", + "config_file_name": "../precice-config.xml", + "interface": { + "coupling_mesh_name": "Dirichlet-Mesh", + "write_data_name": "Heat-Flux", + "read_data_name": "Temperature" + } +} diff --git a/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json b/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json new file mode 100644 index 0000000..2afd52c --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json @@ -0,0 +1,9 @@ +{ + "participant_name": "Neumann", + "config_file_name": "../precice-config.xml", + "interface": { + "coupling_mesh_name": "Neumann-Mesh", + "write_data_name": "Temperature", + "read_data_name": "Heat-Flux" + } +} diff --git a/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py b/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py new file mode 100644 index 0000000..9dbbf74 --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py @@ -0,0 +1,46 @@ +""" +Problem setup for partitioned-heat-conduction/fenicsx tutorial +""" +from dolfinx.mesh import DiagonalType, create_rectangle +from my_enums import DomainPart +import numpy as np + + +y_bottom, y_top = 0, 1 +x_left, x_right = 0, 2 +x_coupling = 1.0 # x coordinate of coupling interface + + +def exclude_straight_boundary(x): + tol = 1E-14 + return np.logical_or( + np.logical_or(~np.isclose(x[0], x_coupling, tol), np.isclose(x[1], y_top, tol)), + np.isclose(x[1], y_bottom, tol) + ) + + +def straight_boundary(x): + tol = 1E-14 + return np.isclose(x[0], x_coupling, tol) + + +def get_geometry(mpi_comm, domain_part): + nx = ny = 9 + low_resolution = 5 + high_resolution = 5 + n_vertices = 20 + + if domain_part is DomainPart.LEFT: + p0 = (x_left, y_bottom) + p1 = (x_coupling, y_top) + elif domain_part is DomainPart.RIGHT: + p0 = (x_coupling, y_bottom) + p1 = (x_right, y_top) + else: + raise Exception("invalid domain_part: {}".format(domain_part)) + + mesh = create_rectangle(mpi_comm, [p0, p1], [nx, ny], diagonal=DiagonalType.left) + coupling_boundary = straight_boundary + remaining_boundary = exclude_straight_boundary + + return mesh, coupling_boundary, remaining_boundary diff --git a/tutorials/partitioned-heat-conduction/fenicsx/run.sh b/tutorials/partitioned-heat-conduction/fenicsx/run.sh new file mode 100755 index 0000000..e31f07a --- /dev/null +++ b/tutorials/partitioned-heat-conduction/fenicsx/run.sh @@ -0,0 +1,16 @@ +#!/bin/sh +set -e -u + +while getopts ":dn" opt; do + case ${opt} in + d) + python3 heat.py -d --error-tol 10e-3 + ;; + n) + python3 heat.py -n --error-tol 10e-3 + ;; + \?) + echo "Usage: cmd [-d] [-n]" + ;; + esac +done diff --git a/tutorials/partitioned-heat-conduction/precice-config.xml b/tutorials/partitioned-heat-conduction/precice-config.xml new file mode 100644 index 0000000..baab5d5 --- /dev/null +++ b/tutorials/partitioned-heat-conduction/precice-config.xml @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +