From d1a041a449ba8652493cd74dbf0d3584b2907137 Mon Sep 17 00:00:00 2001 From: PhilipHildebrand <99822700+PhilipHildebrand@users.noreply.github.com> Date: Mon, 23 Jan 2023 10:58:05 +0100 Subject: [PATCH] Adding FEniCSx solver to partitioned-heat-equation tutorial (#317) * adding fenicsx solver to partitioned-heat-equation tutorial * Delete printStats.py * cleanup work * changed README.md * output reference solution * Changed output folder --- partitioned-heat-conduction/README.md | 14 +- partitioned-heat-conduction/fenicsx/clean.sh | 6 + .../fenicsx/errorcomputation.py | 16 ++ partitioned-heat-conduction/fenicsx/heat.py | 252 ++++++++++++++++++ .../fenicsx/my_enums.py | 17 ++ .../fenicsx/precice-adapter-config-D.json | 9 + .../fenicsx/precice-adapter-config-N.json | 9 + .../fenicsx/problem_setup.py | 46 ++++ partitioned-heat-conduction/fenicsx/run.sh | 16 ++ 9 files changed, 380 insertions(+), 5 deletions(-) create mode 100755 partitioned-heat-conduction/fenicsx/clean.sh create mode 100644 partitioned-heat-conduction/fenicsx/errorcomputation.py create mode 100644 partitioned-heat-conduction/fenicsx/heat.py create mode 100644 partitioned-heat-conduction/fenicsx/my_enums.py create mode 100644 partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json create mode 100644 partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json create mode 100644 partitioned-heat-conduction/fenicsx/problem_setup.py create mode 100755 partitioned-heat-conduction/fenicsx/run.sh diff --git a/partitioned-heat-conduction/README.md b/partitioned-heat-conduction/README.md index 616ea02dc..8b1713d36 100644 --- a/partitioned-heat-conduction/README.md +++ b/partitioned-heat-conduction/README.md @@ -1,7 +1,7 @@ --- title: Partitioned heat conduction permalink: tutorials-partitioned-heat-conduction.html -keywords: FEniCS, Nutils, Heat conduction +keywords: FEniCSx, FEniCS, Nutils, Heat conduction summary: We solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion. --- @@ -25,6 +25,8 @@ This simple case allows us to compare the solution for the partitioned case to a You can either couple a solver with itself or different solvers with each other. In any case you will need to have preCICE and the python bindings installed on your system. +* FEniCSx. Install [FEniCSx](https://fenicsproject.org/download/) and the [FEniCSx-adapter](https://github.com/precice/fenicsx-adapter). The code is largely based on this [fenics-tutorial](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) from [1] and has been adapted to FEniCSx. + * FEniCS. Install [FEniCS](https://fenicsproject.org/download/) and the [FEniCS-adapter](https://github.com/precice/fenics-adapter). The code is largely based on this [fenics-tutorial](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) from [1]. * Nutils. Install [Nutils](https://nutils.org/install-nutils.html). @@ -43,18 +45,18 @@ For choosing whether you want to run the Dirichlet-kind and a Neumann-kind parti For running the case, open two terminals run: ```bash -cd fenics +cd fenicsx ./run.sh -d ``` and ```bash -cd fenics +cd fenicsx ./run.sh -n ``` -If you want to use Nutils for one or both sides of the setup, just `cd nutils`. The FEniCS case also supports parallel runs. Here, you cannot use the `run.sh` script, but must simply execute +If you want to use FEniCS or Nutils, use `cd fenics`/ `cd nutils` instead of `cd fenicsx`. The FEniCS case also supports parallel runs. Here, you cannot use the `run.sh` script, but must simply execute ```bash mpirun -n heat.py -d @@ -66,7 +68,9 @@ You can mix the Nutils and FEniCS solver, if you like. Note that the error for a ## Visualization -Output is written into the folders `fenics/out` and `nutils`. +Output is written into the folders `fenicsx/out`, `fenics/out` and `nutils`. + +For FEniCSx you can visualize the content with paraview by opening the `*.xdmf` files. The files `Dirichlet.xdmf` and `Neumann.xdmf` correspond to the numerical solution of the Dirichlet, respectively Neumann, problem, while the files with the prefix `ref` correspond to the analytical reference solution. For FEniCS you can visualize the content with paraview by opening the `*.pvd` files. The files `Dirichlet.pvd` and `Neumann.pvd` correspond to the numerical solution of the Dirichlet, respectively Neumann, problem, while the files with the prefix `ref` correspond to the analytical reference solution, the files with `error` show the error and the files with `ranks` the ranks of the solvers (if executed in parallel). diff --git a/partitioned-heat-conduction/fenicsx/clean.sh b/partitioned-heat-conduction/fenicsx/clean.sh new file mode 100755 index 000000000..3a8b4619d --- /dev/null +++ b/partitioned-heat-conduction/fenicsx/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_fenics . diff --git a/partitioned-heat-conduction/fenicsx/errorcomputation.py b/partitioned-heat-conduction/fenicsx/errorcomputation.py new file mode 100644 index 000000000..069dd879f --- /dev/null +++ b/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/partitioned-heat-conduction/fenicsx/heat.py b/partitioned-heat-conduction/fenicsx/heat.py new file mode 100644 index 000000000..ccaf70622 --- /dev/null +++ b/partitioned-heat-conduction/fenicsx/heat.py @@ -0,0 +1,252 @@ +""" +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.3 # 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/partitioned-heat-conduction/fenicsx/my_enums.py b/partitioned-heat-conduction/fenicsx/my_enums.py new file mode 100644 index 000000000..5708badb3 --- /dev/null +++ b/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/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json b/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json new file mode 100644 index 000000000..7989ba340 --- /dev/null +++ b/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/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json b/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json new file mode 100644 index 000000000..2afd52ce3 --- /dev/null +++ b/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/partitioned-heat-conduction/fenicsx/problem_setup.py b/partitioned-heat-conduction/fenicsx/problem_setup.py new file mode 100644 index 000000000..9dbbf7437 --- /dev/null +++ b/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/partitioned-heat-conduction/fenicsx/run.sh b/partitioned-heat-conduction/fenicsx/run.sh new file mode 100755 index 000000000..e31f07a10 --- /dev/null +++ b/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