Add nutils solver for heat conduction problem (#153)

* Strip case of everything besides simple heat conduction setup.

* Adds precice-config.xml supporting parallel runs.

* Add documentation and nutils case.

* Add complex case. Based on ae77ae5.

* Cleanup and Update to newest version of python bindings.

* Use scalar-valued flux.

* Apply new structure to complex case.

* Remove subcycling.

* Simplify heat nutils to one mesh, make config consistent

* Fix names.

* BCs working.

* Minor restructuring

* Fix initial condition and output.

* Mark nutils case as under construction. See #152.

* Add nutils prototype from 679fe8d.

* Add instructions on running nutils case to README

* Cleanup

* Fix geometry.

* Add convenience scripts and improve output.

* Make output more consistent.

* Update convenience scripts.

* Fix sign error.

* Add monolithic solver for comparison.

* Scale flux by element size.

* Use identical mesh for FEniCS and Nutils.

* Add argument for setting error tolerance.

* Update README.md

* Remove monolithic case.

* Smooth readme of partitioned heat case

* More information on error.

* Provide cleaning scripts according to #169.

* Streamline run and clean scripts.

* Cleanup.

* Align Nutils vtk output format with Paraview

* Make explanation of fenics-nutils error more precise

* Comment, format, and rewrite nutils heat code

Co-authored-by: uekerman <[email protected]>
Co-authored-by: Benjamin Uekermann <[email protected]>

Author: BenjaminRodenberg

partitioned-heat-conduction/README.md

@@ -2,14 +2,14 @@
 title: Partitioned heat conduction
 permalink: tutorials-partitioned-heat-conduction.html
 keywords: FEniCS, Nutils, Heat conduction
-summary: In this tutorial we solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion.
+summary: We solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion.
 ---
 
 {% include important.html content="We have not yet ported the documentation of the preCICE tutorials from the preCICE wiki to here. Please go to the [preCICE wiki](https://github.com/precice/precice/wiki#2-getting-started---tutorials)" %}
 
 ## Setup
 
-In this tutorial 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].
+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)
 
@@ -26,38 +26,48 @@ You can either couple a solver with itself or different solvers with each other.
 * `fenics`, requires you to 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].
 
 
-* :construction: This case is still under construction. See https://github.com/precice/tutorials/issues/152. :construction: `nutils`, requires you to install [Nutils](http://www.nutils.org/en/latest/).
+* `nutils`, requires you to install [Nutils](http://www.nutils.org/en/latest/).
 
 ## 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 and:
+For running the case, open two terminals run:
 
 ```
 cd fenics
-python3 heat.py -d
+python3 ./run.sh -d
 ```
 
 and
 
 ```
 cd fenics
-python3 heat.py -n
+python3 ./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, simply execute
+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
 
 ```
 mpirun -n <N_PROC> python3 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 folder `fenics/out`. 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).
+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)
 

partitioned-heat-conduction/clean-tutorial.sh

@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../tools/cleaning-tools.sh
+
+clean_tutorial .

partitioned-heat-conduction/fenics/Allclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+rm -rf out 
+rm -f *.log
+rm -f *-events.json

partitioned-heat-conduction/fenics/clean.sh

@@ -61,12 +61,14 @@ def determine_gradient(V_g, u, flux):
         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 = 10 ** -6
+error_tol = args.error_tol
 
 alpha = 3  # parameter alpha
 beta = 1.3  # parameter beta

partitioned-heat-conduction/fenics/heat.py

@@ -35,24 +35,21 @@ def inside(self, x, on_boundary):
 
 
 def get_geometry(domain_part):
-    nx = 5
-    ny = 10
+    nx = ny = 9
     low_resolution = 5
     high_resolution = 5
     n_vertices = 20
 
     if domain_part is DomainPart.LEFT:
-        nx = nx*3
         p0 = Point(x_left, y_bottom)
         p1 = Point(x_coupling, y_top)        
     elif domain_part is DomainPart.RIGHT:
-        ny = ny*2
         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)
+    mesh = RectangleMesh(p0, p1, nx, ny, diagonal="left")
     coupling_boundary = StraightBoundary()
     remaining_boundary = ExcludeStraightBoundary()
 

partitioned-heat-conduction/fenics/problem_setup.py

@@ -0,0 +1,11 @@
+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
+

partitioned-heat-conduction/fenics/run.sh

@@ -0,0 +1,4 @@
+Neumann-*.vtk
+Dirichlet-*.vtk
+*.log
+precice-*-events.json

partitioned-heat-conduction/nutils/.gitignore

@@ -0,0 +1,9 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+rm -f Dirichlet-*.vtk
+rm -f Neumann-*.vtk
+rm -f *.log
+rm -f *-events.json

partitioned-heat-conduction/nutils/clean.sh

@@ -0,0 +1,170 @@
+#! /usr/bin/env python3
+
+import nutils, treelog
+import numpy as np
+import precice
+
+
+def main(side='Dirichlet'):
+
+  print("Running nutils")
+
+  # domain size
+  y_bottom, y_top = 0, 1
+  x_left, x_right = 0, 2
+  x_coupling = 1  # x coordinate of coupling interface
+
+  n = 10 # number of mesh vertices per dimension 
+
+  if side == 'Dirichlet':
+    x_grid = np.linspace(x_left, x_coupling, n)
+  elif side == 'Neumann':
+    x_grid = np.linspace(x_coupling, x_right, n)
+  else:
+    raise Exception('invalid side {!r}'.format(side))
+
+  y_grid = np.linspace(y_bottom, y_top, n)
+
+  # define the Nutils mesh
+  domain, geom = nutils.mesh.rectilinear([x_grid, y_grid])
+
+  # Nutils namespace
+  ns = nutils.function.Namespace()
+  ns.x = geom
+  
+  degree = 1 # linear finite elements
+  ns.basis = domain.basis('std', degree=degree) 
+  ns.alpha = 3 # parameter of problem
+  ns.beta = 1.3 # parameter of problem
+  ns.u = 'basis_n ?lhs_n' # solution
+  ns.dudt = 'basis_n (?lhs_n - ?lhs0_n) / ?dt' # time derivative
+  ns.flux = 'basis_n ?fluxdofs_n' # heat flux
+  ns.f = 'beta - 2 - 2 alpha' # rhs
+  ns.uexact = '1 + x_0 x_0 + alpha x_1 x_1 + beta ?t' # analytical solution
+
+  # define the weak form
+  res0 = domain.integral('(basis_n dudt - basis_n f + basis_n,i u_,i) d:x' @ ns, degree=degree * 2)
+
+  # set boundary conditions at non-coupling boundaries
+  # top and bottom boundary are non-coupling for both sides
+  sqr0 = domain.boundary['top'].integral('(u - 1 - x_0 x_0 - alpha - beta ?t)^2 d:x' @ ns, degree=degree * 2)
+  sqr0 += domain.boundary['bottom'].integral('(u - 1 - x_0 x_0 - beta ?t)^2 d:x' @ ns, degree=degree * 2)
+  if side == 'Dirichlet':  # left boundary is non-coupling
+    sqr0 += domain.boundary['left'].integral('(u - 1 - alpha x_1 x_1 - beta ?t)^2 d:x' @ ns, degree=degree * 2)
+  elif side == 'Neumann':  # right boundary is non-coupling
+    sqr0 += domain.boundary['right'].integral('(u - 1 - x_0 x_0 - alpha x_1 x_1 - beta ?t)^2 d:x' @ ns, degree=degree * 2)
+
+
+  # preCICE setup
+  interface = precice.Interface(side, "../precice-config.xml", 0, 1)
+  
+  # define coupling mesh
+  mesh_name = side + "-Mesh"
+  mesh_id = interface.get_mesh_id(mesh_name)
+  coupling_boundary = domain.boundary['right' if side == 'Dirichlet' else 'left']
+  coupling_sample = coupling_boundary.sample('gauss', degree=degree * 2)
+  vertices = coupling_sample.eval(ns.x)
+  vertex_ids = interface.set_mesh_vertices(mesh_id, vertices)
+
+  # coupling data
+  write_data = "Temperature" if side == "Neumann" else "Flux"
+  read_data = "Flux" if side == "Neumann" else "Temperature"
+  write_data_id = interface.get_data_id(write_data, mesh_id)
+  read_data_id = interface.get_data_id(read_data, mesh_id)
+
+  # helper functions to project heat flux to coupling boundary
+  projection_matrix = coupling_boundary.integrate(ns.eval_nm('basis_n basis_m d:x'), degree=degree * 2)
+  projection_cons = np.zeros(res0.shape)
+  projection_cons[projection_matrix.rowsupp(1e-15)] = np.nan
+  fluxdofs = lambda v: projection_matrix.solve(v, constrain=projection_cons)
+
+  # helper data structures to integrate heat flux correctly  
+  dx = coupling_sample.eval('d:x' @ ns)
+  dx_function = coupling_sample.asfunction(dx)
+  
+  precice_dt = interface.initialize()
+
+  # write initial data
+  if interface.is_action_required(precice.action_write_initial_data()):
+    write_data = np.zeros(len(vertex_ids))
+    interface.write_block_scalar_data(write_data_id, vertex_ids, write_data)
+    interface.mark_action_fulfilled(precice.action_write_initial_data())
+
+  interface.initialize_data()
+
+  t = 0
+
+  # initial condition
+  sqr = domain.integral('(u - uexact)^2' @ ns, degree=degree*2)
+  lhs0 = nutils.solver.optimize('lhs', sqr, droptol=1e-15, arguments=dict(t=t))
+  bezier = domain.sample('bezier', degree * 2)
+  x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs0, t=t)
+  with treelog.add(treelog.DataLog()):
+    nutils.export.vtk(side + '-0', bezier.tri, x, Temperature=u, reference=uexact)
+                        
+  t += precice_dt
+  timestep = 1
+  dt = 0.1
+
+  while interface.is_coupling_ongoing():
+  
+    # update (time-dependent) boundary condition
+    cons0 = nutils.solver.optimize('lhs', sqr0, droptol=1e-15, arguments=dict(t=t))
+
+    # read data from interface
+    if interface.is_read_data_available():
+      read_data = interface.read_block_scalar_data(read_data_id, vertex_ids)
+      read_function = coupling_sample.asfunction(read_data)
+
+      if side == 'Dirichlet':
+        sqr = coupling_sample.integral((ns.u - read_function)**2)
+        cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons0, arguments=dict(t=t))
+        res = res0
+      else:
+        cons = cons0
+        res = res0 + coupling_sample.integral(ns.basis * read_function * dx_function)
+            
+    # save checkpoint
+    if interface.is_action_required(precice.action_write_iteration_checkpoint()):
+      lhs_checkpoint = lhs0
+      interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
+
+    # potentially adjust non-matching timestep sizes 
+    dt = min(dt, precice_dt)
+
+    # solve nutils timestep
+    lhs = nutils.solver.solve_linear('lhs', res, constrain=cons, arguments=dict(lhs0=lhs0, dt=dt, t=t))
+
+    # write data to interface
+    if interface.is_write_data_required(dt):
+      if side == 'Dirichlet':
+        flux_function = res.eval(lhs0=lhs0, lhs=lhs, dt=dt, t=t)
+        write_data = coupling_sample.eval('flux' @ ns, fluxdofs=fluxdofs(flux_function))
+      else:
+        write_data = coupling_sample.eval('u' @ ns, lhs=lhs)
+
+      interface.write_block_scalar_data(write_data_id, vertex_ids, write_data)
+
+    # do the coupling
+    precice_dt = interface.advance(dt)
+
+    # read checkpoint if required
+    if interface.is_action_required(precice.action_read_iteration_checkpoint()):
+      lhs0 = lhs_checkpoint
+      interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
+    else: # go to next timestep
+      bezier = domain.sample('bezier', degree * 2)
+      x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t)
+
+      with treelog.add(treelog.DataLog()):
+        nutils.export.vtk(side + "-" + str(timestep), bezier.tri, x, Temperature=u, reference=uexact)
+                                
+      t += dt
+      timestep += 1
+      lhs0 = lhs
+
+  interface.finalize()
+
+
+if __name__ == '__main__':
+    nutils.cli.run(main)

partitioned-heat-conduction/nutils/heat.py

@@ -0,0 +1,13 @@
+while getopts ":dn" opt; do
+  case ${opt} in
+    d ) rm -rf Dirichlet-*.vtk
+        python3 heat.py --side=Dirichlet
+      ;;
+    n ) rm -rf Neumann-*.vtk
+	python3 heat.py --side=Neumann
+      ;;
+    \? ) echo "Usage: cmd [-d] [-n]"
+      ;;
+  esac
+done
+