diff --git a/tutorials/partitioned-heat-conduction/README.md b/tutorials/partitioned-heat-conduction/README.md deleted file mode 100644 index 8f98175..0000000 --- a/tutorials/partitioned-heat-conduction/README.md +++ /dev/null @@ -1,83 +0,0 @@ ---- -title: Partitioned heat conduction -permalink: tutorials-partitioned-heat-conduction.html -keywords: 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. ---- - -{% note %} -Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/partitioned-heat-conduction). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html). -{% endnote %} - -## Setup - -We solve a partitioned heat equation. For information on the non-partitioned case, please refer to [1, p.37ff]. In this tutorial the computational domain is partitioned and coupled via preCICE. The coupling roughly follows the approach described in [2]. - -![Case setup of partitioned-heat-conduction case](images/tutorials-partitioned-heat-conduction-setup.png) - -Case setup from [3]. `D` denotes the Dirichlet participant and `N` denotes the Neumann participant. - -The heat equation is solved on a rectangular domain `Omega = [0,2] x [0,1]` with given Dirichlet boundary conditions. We split the domain at `x_c = 1` using a straight vertical line, the coupling interface. The left part of the domain will be referred to as the Dirichlet partition and the right part as the Neumann partition. To couple the two participants we use Dirichlet-Neumann coupling. Here, the Dirichlet participant receives Dirichlet boundary conditions (`Temperature`) at the coupling interface and solves the heat equation using these boundary conditions on the left part of the domain. Then the Dirichlet participant computes the resulting heat flux (`Flux`) from the solution and sends it to the Neumann participant. The Neumann participant uses the flux as a Neumann boundary condition to solve the heat equation on the right part of the domain. We then extract the temperature from the solution and send it back to the Dirichlet participant. This establishes the coupling between the two participants. - -This simple case allows us to compare the solution for the partitioned case to a known analytical solution (method of manufactures solutions, see [1, p.37ff]). For more usage examples and details, please refer to [3, sect. 4.1]. - -## Available solvers and dependencies - -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. - -* 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). - -* OpenFOAM. This case also requires [funkySetFields](https://openfoamwiki.net/index.php/Contrib/funkySetFields) (part of [swak4Foam](https://openfoamwiki.net/index.php/Contrib/swak4Foam)) and uses the custom [heatTransfer](https://github.com/precice/tutorials/blob/master/partitioned-heat-conduction/openfoam-solver/heatTransfer.C) solver (find it in `openfoam-solver` and build with `wmake`). Read more details in the [OpenFOAM adapter](https://precice.org/adapter-openfoam-overview.html). - -## Running the simulation - -This tutorial is for FEniCS and Nutils. You can find the corresponding `run.sh` script in the folders `fenics` and `nutils`. - -For choosing whether you want to run the Dirichlet-kind and a Neumann-kind participant, please provide the following commandline input: - -* `-d` flag will enforce Dirichlet boundary conditions on the coupling interface. -* `-n` flag will enforce Neumann boundary conditions on the coupling interface. - -For running the case, open two terminals run: - -```bash -cd fenics -./run.sh -d -``` - -and - -```bash -cd fenics -./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 - -```bash -mpirun -n heat.py -d -``` - -### Note on the combination of Nutils & FEniCS - -You can mix the Nutils and FEniCS solver, if you like. Note that the error for a pure FEniCS simulation is lower than for a mixed one. We did not yet study the origin of this error, but assume that this is due to the fact that Nutils uses Gauss points as coupling mesh and therefore entails extrapolation in the data mapping at the top and bottom corners. - -## Visualization - -Output is written into the folders `fenics/out` and `nutils`. - -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). - -For Nutils, please use the files `Dirichlet-*.vtk` or `Neumann-*.vtk`. Please note that these files contain the temperature as well as the reference solution. - -![Animation of the partitioned heat equation](images/tutorials-partitioned-heat-conduction-FEniCS-movie.gif) - -Visualization in paraview for `x_c = 1.5`. - -## References - -[1] Hans Petter Langtangen and Anders Logg. "Solving PDEs in Minutes-The FEniCS Tutorial Volume I." (2016). [pdf](https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf) -[2] Azahar Monge and Philipp Birken. "Convergence Analysis of the Dirichlet-Neumann Iteration for Finite Element Discretizations." (2016). Proceedings in Applied Mathematics and Mechanics. [doi](https://doi.org/10.1002/pamm.201610355) -[3] Benjamin Rüth, Benjamin Uekermann, Miriam Mehl, Philipp Birken, Azahar Monge, and Hans Joachim Bungartz. "Quasi-Newton waveform iteration for partitioned surface-coupled multiphysics applications." (2020). International Journal for Numerical Methods in Engineering. [doi](https://doi.org/10.1002/nme.6443) diff --git a/tutorials/partitioned-heat-conduction/clean-tutorial.sh b/tutorials/partitioned-heat-conduction/clean-tutorial.sh deleted file mode 120000 index 4713f50..0000000 --- a/tutorials/partitioned-heat-conduction/clean-tutorial.sh +++ /dev/null @@ -1 +0,0 @@ -../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/tutorials/partitioned-heat-conduction/fenics/.gitignore b/tutorials/partitioned-heat-conduction/fenics/.gitignore deleted file mode 100644 index 7c6571e..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -venv -*.pyc diff --git a/tutorials/partitioned-heat-conduction/fenics/clean.sh b/tutorials/partitioned-heat-conduction/fenics/clean.sh deleted file mode 100755 index 3a8b461..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_fenics . diff --git a/tutorials/partitioned-heat-conduction/fenics/errorcomputation.py b/tutorials/partitioned-heat-conduction/fenics/errorcomputation.py deleted file mode 100644 index 2c9ec73..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/errorcomputation.py +++ /dev/null @@ -1,15 +0,0 @@ -from fenics import inner, assemble, dx, project, sqrt - - -def compute_errors(u_approx, u_ref, v, total_error_tol=10 ** -4): - # compute pointwise L2 error - error_normalized = (u_ref - u_approx) / u_ref - # project onto function space - error_pointwise = project(abs(error_normalized), v) - # determine L2 norm to estimate total error - error_total = sqrt(assemble(inner(error_pointwise, error_pointwise) * dx)) - error_pointwise.rename("error", " ") - - assert (error_total < total_error_tol) - - return error_total, error_pointwise diff --git a/tutorials/partitioned-heat-conduction/fenics/heat.py b/tutorials/partitioned-heat-conduction/fenics/heat.py deleted file mode 100644 index e7113f5..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/heat.py +++ /dev/null @@ -1,230 +0,0 @@ -""" -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. - -The original source code can be found on https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py. - -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 fenics import Function, FunctionSpace, Expression, Constant, DirichletBC, TrialFunction, TestFunction, \ - File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace, MeshFunction, MPI -from fenicsprecice 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 -import dolfin -from dolfin import FacetNormal, dot - - -def determine_gradient(V_g, u, flux): - """ - 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 - :param flux: returns calculated flux into this value - """ - - w = TrialFunction(V_g) - v = TestFunction(V_g) - - a = inner(w, v) * dx - L = inner(grad(u), v) * dx - solve(a == L, flux) - - -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(domain_part) - -# Define function space using mesh -V = FunctionSpace(mesh, 'P', 2) -V_g = VectorFunctionSpace(mesh, 'P', 1) -W = V_g.sub(0).collapse() - -# Define boundary conditions -u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t', degree=2, alpha=alpha, beta=beta, t=0) -u_D_function = interpolate(u_D, V) - -if problem is ProblemType.DIRICHLET: - # Define flux in x direction - f_N = Expression("2 * x[0]", degree=1, alpha=alpha, t=0) - f_N_function = interpolate(f_N, W) - -# Define initial value -u_n = interpolate(u_D, V) -u_n.rename("Temperature", "") - -precice, precice_dt, initial_data = None, 0.0, None - -# Initialize the adapter according to the specific participant -if problem is ProblemType.DIRICHLET: - precice = Adapter(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(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(0) -dt.assign(np.min([fenics_dt, precice_dt])) - -# Define variational problem -u = TrialFunction(V) -v = TestFunction(V) -f = Expression('beta - 2 - 2*alpha', degree=2, alpha=alpha, beta=beta, t=0) -F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f) * v * dx - -bcs = [DirichletBC(V, u_D, remaining_boundary)] - -# Set boundary conditions at coupling interface once wrt to the coupling -# expression -coupling_expression = precice.create_coupling_expression() -if problem is ProblemType.DIRICHLET: - # modify Dirichlet boundary condition on coupling interface - bcs.append(DirichletBC(V, coupling_expression, coupling_boundary)) -if problem is ProblemType.NEUMANN: - # modify Neumann boundary condition on coupling interface, modify weak - # form correspondingly - F += v * coupling_expression * dolfin.ds - -a, L = lhs(F), rhs(F) - -# Time-stepping -u_np1 = Function(V) -u_np1.rename("Temperature", "") -t = 0 - -# reference solution at t=0 -u_ref = interpolate(u_D, V) -u_ref.rename("reference", " ") - -# mark mesh w.r.t ranks -mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim()) -if problem is ProblemType.NEUMANN: - mesh_rank.set_all(MPI.rank(MPI.comm_world) + 4) -else: - mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0) -mesh_rank.rename("myRank", "") - -# Generating output files -temperature_out = File("out/%s.pvd" % precice.get_participant_name()) -ref_out = File("out/ref%s.pvd" % precice.get_participant_name()) -error_out = File("out/error%s.pvd" % precice.get_participant_name()) -ranks = File("out/ranks%s.pvd" % precice.get_participant_name()) - -# output solution and reference solution at t=0, n=0 -n = 0 -print('output u^%d and u_ref^%d' % (n, n)) -temperature_out << u_n -ref_out << u_ref -ranks << mesh_rank - -error_total, error_pointwise = compute_errors(u_n, u_ref, V) -error_out << error_pointwise - -# set t_1 = t_0 + dt, this gives u_D^1 -# call dt(0) to evaluate FEniCS Constant. Todo: is there a better way? -u_D.t = t + dt(0) -f.t = t + dt(0) - -if problem is ProblemType.DIRICHLET: - flux = Function(V_g) - flux.rename("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.assign(np.min([fenics_dt, precice_dt])) - - # Compute solution u^n+1, use bcs u_D^n+1, u^n and coupling bcs - solve(a == L, u_np1, bcs) - - # 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 - determine_gradient(V_g, u_np1, flux) - flux_x = interpolate(flux.sub(0), W) - 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(0)) - - # 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.assign(u_cp) - t = t_cp - n = n_cp - else: # update solution - u_n.assign(u_np1) - t += float(dt) - n += 1 - - if precice.is_time_window_complete(): - u_ref = interpolate(u_D, V) - u_ref.rename("reference", " ") - error, error_pointwise = compute_errors(u_n, u_ref, V, total_error_tol=error_tol) - print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error)) - # output solution and reference solution at t_n+1 - print('output u^%d and u_ref^%d' % (n, n)) - temperature_out << u_n - ref_out << u_ref - error_out << error_pointwise - - # Update Dirichlet BC - u_D.t = t + float(dt) - f.t = t + float(dt) - -# Hold plot -precice.finalize() diff --git a/tutorials/partitioned-heat-conduction/fenics/my_enums.py b/tutorials/partitioned-heat-conduction/fenics/my_enums.py deleted file mode 100644 index 27fc5d9..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/my_enums.py +++ /dev/null @@ -1,19 +0,0 @@ -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 - CIRCULAR = 3 # circular part of domain in complex interface case - RECTANGLE = 4 # domain excluding circular part of complex interface case diff --git a/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-D.json b/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-D.json deleted file mode 100644 index c2a7788..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-D.json +++ /dev/null @@ -1,9 +0,0 @@ -{ - "participant_name": "Dirichlet", - "config_file_name": "../precice-config.xml", - "interface": { - "coupling_mesh_name": "Dirichlet-Mesh", - "write_data_name": "Flux", - "read_data_name": "Temperature" - } -} diff --git a/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-N.json b/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-N.json deleted file mode 100644 index 36ea15d..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/precice-adapter-config-N.json +++ /dev/null @@ -1,9 +0,0 @@ -{ - "participant_name": "Neumann", - "config_file_name": "../precice-config.xml", - "interface": { - "coupling_mesh_name": "Neumann-Mesh", - "write_data_name": "Temperature", - "read_data_name": "Flux" - } -} diff --git a/tutorials/partitioned-heat-conduction/fenics/printStats.py b/tutorials/partitioned-heat-conduction/fenics/printStats.py deleted file mode 100644 index baa0f83..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/printStats.py +++ /dev/null @@ -1,4 +0,0 @@ -import pstats -from pstats import SortKey -p = pstats.Stats('Dirichlet.prof') -p.strip_dirs().sort_stats(1).print_stats(20) diff --git a/tutorials/partitioned-heat-conduction/fenics/problem_setup.py b/tutorials/partitioned-heat-conduction/fenics/problem_setup.py deleted file mode 100644 index cb1b164..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/problem_setup.py +++ /dev/null @@ -1,56 +0,0 @@ -""" -Problem setup for partitioned-heat-conduction/fenics-fenics tutorial -""" - -from fenics import SubDomain, Point, RectangleMesh, near, Function, VectorFunctionSpace, Expression -from my_enums import DomainPart - - -y_bottom, y_top = 0, 1 -x_left, x_right = 0, 2 -x_coupling = 1.0 # x coordinate of coupling interface -radius = 0.2 -midpoint = Point(0.5, 0.5) - - -class ExcludeStraightBoundary(SubDomain): - def get_user_input_args(self, args): - self._interface = args.interface - - def inside(self, x, on_boundary): - tol = 1E-14 - if on_boundary and not near(x[0], x_coupling, tol) or near(x[1], y_top, tol) or near(x[1], y_bottom, tol): - return True - else: - return False - - -class StraightBoundary(SubDomain): - def inside(self, x, on_boundary): - tol = 1E-14 - if on_boundary and near(x[0], x_coupling, tol): - return True - else: - return False - - -def get_geometry(domain_part): - nx = ny = 9 - low_resolution = 5 - high_resolution = 5 - n_vertices = 20 - - if domain_part is DomainPart.LEFT: - p0 = Point(x_left, y_bottom) - p1 = Point(x_coupling, y_top) - elif domain_part is DomainPart.RIGHT: - p0 = Point(x_coupling, y_bottom) - p1 = Point(x_right, y_top) - else: - raise Exception("invalid domain_part: {}".format(domain_part)) - - mesh = RectangleMesh(p0, p1, nx, ny, diagonal="left") - coupling_boundary = StraightBoundary() - remaining_boundary = ExcludeStraightBoundary() - - return mesh, coupling_boundary, remaining_boundary diff --git a/tutorials/partitioned-heat-conduction/fenics/run.sh b/tutorials/partitioned-heat-conduction/fenics/run.sh deleted file mode 100755 index e31f07a..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/run.sh +++ /dev/null @@ -1,16 +0,0 @@ -#!/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/fenics/tutorials.patch b/tutorials/partitioned-heat-conduction/fenics/tutorials.patch deleted file mode 100644 index 8dc831f..0000000 --- a/tutorials/partitioned-heat-conduction/fenics/tutorials.patch +++ /dev/null @@ -1,33 +0,0 @@ -diff --git a/partitioned-heat-conduction/fenics/problem_setup.py b/partitioned-heat-conduction/fenics/problem_setup.py -index cb1b164..e8073db 100644 ---- a/partitioned-heat-conduction/fenics/problem_setup.py -+++ b/partitioned-heat-conduction/fenics/problem_setup.py -@@ -35,7 +35,7 @@ class StraightBoundary(SubDomain): - - - def get_geometry(domain_part): -- nx = ny = 9 -+ nx = ny = 3 - low_resolution = 5 - high_resolution = 5 - n_vertices = 20 -diff --git a/partitioned-heat-conduction/precice-config.xml b/partitioned-heat-conduction/precice-config.xml -index a810f28..b0262fa 100644 ---- a/partitioned-heat-conduction/precice-config.xml -+++ b/partitioned-heat-conduction/precice-config.xml -@@ -41,13 +41,14 @@ - - - -- -+ - - - - - - -+ - - - diff --git a/tutorials/partitioned-heat-conduction/fenicsx/.gitignore b/tutorials/partitioned-heat-conduction/fenicsx/.gitignore deleted file mode 100644 index b48121a..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/.gitignore +++ /dev/null @@ -1,5 +0,0 @@ -venv -*.pyc -*.log -*.json -out/ diff --git a/tutorials/partitioned-heat-conduction/fenicsx/clean.sh b/tutorials/partitioned-heat-conduction/fenicsx/clean.sh deleted file mode 100755 index 3a8b461..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/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 deleted file mode 100644 index d2d6ff3..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py +++ /dev/null @@ -1,17 +0,0 @@ -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, error_pointwise - return error_total diff --git a/tutorials/partitioned-heat-conduction/fenicsx/heat.py b/tutorials/partitioned-heat-conduction/fenicsx/heat.py deleted file mode 100644 index cc03160..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/heat.py +++ /dev/null @@ -1,287 +0,0 @@ -""" -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. - -The original source code can be found on https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py. - -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 # TODO update do dolfinx -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 - self.alpha = alpha - self.beta = beta - - def eval(self, x): - return np.full(x.shape[1], 1 + x[0] * x[0] + self.alpha * x[1] * x[1] + self.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): - self.alpha = alpha - self.t = 0.0 - - 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.alpha = alpha - self.beta = beta - self.t = 0.0 - - def eval(self, x): - return np.full(x.shape[1], self.beta - 2 - 2 * self.alpha) - - -f = Expression_f() -f_function = Function(V) -# f_function.interpolate(f.eval) -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)] - -# Set boundary conditions at coupling interface once wrt to the coupling -# expression -# TODO hide that in coupling_expression -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) - -''' -# TODO -# mark mesh w.r.t ranks -mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim()) -if problem is ProblemType.NEUMANN: - mesh_rank.set_all(MPI.rank(MPI.comm_world) + 4) -else: - mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0) -mesh_rank.rename("myRank", "") -''' - -with XDMFFile(MPI.COMM_WORLD, f"./out/{precice.get_participant_name()}.xdmf", "w") as xdmf: - xdmf.write_mesh(mesh) - - # output solution and reference solution at t=0, n=0 - n = 0 - xdmf.write_function(u_n, t) - ''' - print('output u^%d and u_ref^%d' % (n, n)) - temperature_out << u_n - ref_out << u_ref - ranks << mesh_rank - ''' - - # error_total, error_pointwise = compute_errors(u_n, u_ref, V) # TODO - '' - # TODO - ''' - error_out << error_pointwise - ''' - 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="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]) - - # Compute solution u^n+1, use bcs u_D^n+1, u^n and coupling bcs - # , petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) # TODO is it possible to do that only once (before th coupling-loop)? - 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, error_pointwise = compute_errors(mesh, u_n, u_ref, total_error_tol=error_tol) - 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)) - - xdmf.write_function(u_n, t) - - # output solution and reference solution at t_n+1 - # temperature_out << u_n - # ref_out << u_ref - # error_out << error_pointwise - - # 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 deleted file mode 100644 index 27fc5d9..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py +++ /dev/null @@ -1,19 +0,0 @@ -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 - CIRCULAR = 3 # circular part of domain in complex interface case - RECTANGLE = 4 # domain excluding circular part of complex 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 deleted file mode 100644 index c2a7788..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-D.json +++ /dev/null @@ -1,9 +0,0 @@ -{ - "participant_name": "Dirichlet", - "config_file_name": "../precice-config.xml", - "interface": { - "coupling_mesh_name": "Dirichlet-Mesh", - "write_data_name": "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 deleted file mode 100644 index 36ea15d..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/precice-adapter-config-N.json +++ /dev/null @@ -1,9 +0,0 @@ -{ - "participant_name": "Neumann", - "config_file_name": "../precice-config.xml", - "interface": { - "coupling_mesh_name": "Neumann-Mesh", - "write_data_name": "Temperature", - "read_data_name": "Flux" - } -} diff --git a/tutorials/partitioned-heat-conduction/fenicsx/printStats.py b/tutorials/partitioned-heat-conduction/fenicsx/printStats.py deleted file mode 100644 index baa0f83..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/printStats.py +++ /dev/null @@ -1,4 +0,0 @@ -import pstats -from pstats import SortKey -p = pstats.Stats('Dirichlet.prof') -p.strip_dirs().sort_stats(1).print_stats(20) diff --git a/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py b/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py deleted file mode 100644 index 7146bb5..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py +++ /dev/null @@ -1,65 +0,0 @@ -""" -Problem setup for partitioned-heat-conduction/fenics-fenics 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 -radius = 0.2 -midpoint = (0.5, 0.5, 0) - - -# class ExcludeStraightBoundary(SubDomain): -# def get_user_input_args(self, args): -# self._interface = args.interface - -# def inside(self, x, on_boundary): -# tol = 1E-14 -# if on_boundary and not near(x[0], x_coupling, tol) or near(x[1], y_top, tol) or near(x[1], y_bottom, tol): -# return True -# else: -# return False -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) - ) - - -# class StraightBoundary(SubDomain): -# def inside(self, x, on_boundary): -# tol = 1E-14 -# if on_boundary and near(x[0], x_coupling, tol): -# return True -# else: -# return False -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 deleted file mode 100755 index e31f07a..0000000 --- a/tutorials/partitioned-heat-conduction/fenicsx/run.sh +++ /dev/null @@ -1,16 +0,0 @@ -#!/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/images/tutorials-partitioned-heat-conduction-FEniCS-movie.gif b/tutorials/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-FEniCS-movie.gif deleted file mode 100644 index 8726ae1..0000000 Binary files a/tutorials/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-FEniCS-movie.gif and /dev/null differ diff --git a/tutorials/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-setup.png b/tutorials/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-setup.png deleted file mode 100644 index b494c6a..0000000 Binary files a/tutorials/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-setup.png and /dev/null differ diff --git a/tutorials/partitioned-heat-conduction/precice-config.xml b/tutorials/partitioned-heat-conduction/precice-config.xml deleted file mode 100644 index a810f28..0000000 --- a/tutorials/partitioned-heat-conduction/precice-config.xml +++ /dev/null @@ -1,61 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -