From 9cf67afa360d11fe46f901ee96e1bbf757e29fc5 Mon Sep 17 00:00:00 2001 From: TobiasEppacher Date: Thu, 22 Aug 2024 11:20:18 +0200 Subject: [PATCH 01/13] Add files --- .../dirichlet-pySDC-fenics/clean.sh | 6 + .../precice-adapter-config.json | 9 + .../dirichlet-pySDC-fenics/run.sh | 9 + .../solver-fenics/heat_pySDC.py | 346 ++++++++++++++++++ .../solver-fenics/heat_pySDC_problem_class.py | 130 +++++++ 5 files changed, 500 insertions(+) create mode 100755 partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh create mode 100644 partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json create mode 100755 partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh create mode 100644 partitioned-heat-conduction/solver-fenics/heat_pySDC.py create mode 100644 partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py diff --git a/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh b/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh new file mode 100755 index 000000000..d3f3318b4 --- /dev/null +++ b/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_fenics . diff --git a/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json b/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json new file mode 100644 index 000000000..7989ba340 --- /dev/null +++ b/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.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/dirichlet-pySDC-fenics/run.sh b/partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh new file mode 100755 index 000000000..f06eb4084 --- /dev/null +++ b/partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh @@ -0,0 +1,9 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +python3 ../solver-fenics/heat_pySDC.py Dirichlet + +close_log diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py new file mode 100644 index 000000000..5163fc335 --- /dev/null +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py @@ -0,0 +1,346 @@ +""" +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, grad, inner, dx, 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 sympy as sp + +from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from heat_pySDC_problem_class import fenics_heat_2d + +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") +parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in ProblemType]) +parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-8,) + +args = parser.parse_args() +participant_name = args.participantName + + +# preCICE error tolerance +# Error is bounded by coupling accuracy. In theory we would obtain the analytical solution. +error_tol = args.error_tol + +if participant_name == ProblemType.DIRICHLET.value: + problem = ProblemType.DIRICHLET + domain_part = DomainPart.LEFT +elif participant_name == ProblemType.NEUMANN.value: + problem = ProblemType.NEUMANN + domain_part = DomainPart.RIGHT + exit("Neumann problem not yet supported with pySDC") + +domain_mesh, coupling_boundary, remaining_boundary = get_geometry(domain_part) + +# Define function space using mesh +V = FunctionSpace(domain_mesh, 'P', 2) +V_g = VectorFunctionSpace(domain_mesh, 'P', 1) +W = V_g.sub(0).collapse() + + +################################################################## +# Problem definition and preCICE initialization +################################################################## + +# Time step size +# Should be integer fraction of the used time window size +precice_time_window_size = 0.1 +precice_waveform_deg = 1 +pySDC_dt = precice_time_window_size/precice_waveform_deg + +# Manufactured solution parameters +alpha = 3 # parameter alpha +beta = 1.2 # parameter beta + + +# create sympy expression of manufactured solution +x_sp, y_sp, t_sp = sp.symbols(['x[0]', 'x[1]', 't']) + +temporal_deg = 1 # temporal degree of the manufactured solution +u_D_sp = 1 + x_sp * x_sp + alpha * y_sp * y_sp + beta * (t_sp ** temporal_deg) +u_D = Expression(sp.ccode(u_D_sp), degree=2, alpha=alpha, beta=beta, temporal_deg=temporal_deg, t=0) + +u_D_function = interpolate(u_D, V) +f_sp = u_D_sp.diff(t_sp) - u_D_sp.diff(x_sp).diff(x_sp) - u_D_sp.diff(y_sp).diff(y_sp) +forcing_expr = Expression(sp.ccode(f_sp), degree=2, alpha=alpha, beta=beta, t=0) + + +if problem is ProblemType.DIRICHLET: + # Define flux in x direction + f_N = Expression(sp.ccode(u_D_sp.diff(x_sp)), 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 +precice = Adapter(adapter_config_filename="precice-adapter-config.json") + +if problem is ProblemType.DIRICHLET: + precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function) +elif problem is ProblemType.NEUMANN: + precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D_function) + +precice_dt = precice.get_max_time_step_size() +dt = Constant(0) +dt.assign(np.min([pySDC_dt, precice_dt])) + +coupling_expression = precice.create_coupling_expression() + + +# 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", domain_mesh, domain_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("output/%s.pvd" % precice.get_participant_name()) +ref_out = File("output/ref%s.pvd" % precice.get_participant_name()) +error_out = File("output/error%s.pvd" % precice.get_participant_name()) +ranks = File("output/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)) +ranks << mesh_rank + +error_total, error_pointwise = compute_errors(u_n, u_ref, V) + +# create buffer for output. We need this buffer, because we only want to +# write the converged output at the end of the window, but we also want to +# write the samples that are resulting from substeps inside the window +u_write = [] +ref_write = [] +error_write = [] +# copy data to buffer and rename +uu = u_n.copy() +uu.rename("u", "") +u_write.append((uu, t)) +uu_ref = u_ref.copy() +uu_ref.rename("u_ref", "") +ref_write.append(uu_ref) +err = error_pointwise.copy() +err.rename("err", "") +error_write.append(err) + + +if problem is ProblemType.DIRICHLET: + flux = Function(V_g) + flux.rename("Heat-Flux", "") + + +################################################################## +# Setup pySDC controller +################################################################## +''' +pySDC can be installed via `pip install pySDC`. +If you want to use the developer version, the source code repository can be cloned from "https://github.com/Parallel-in-Time/pySDC". +For more information on pySDC, see also "https://parallel-in-time.org/pySDC/" +''' + +# initialize level parameters +level_params = dict() +level_params['restol'] = 1e-11 +level_params['dt'] = pySDC_dt + +# initialize step parameters +step_params = dict() +step_params['maxiter'] = 40 + +# initialize sweeper parameters +sweeper_params = dict() +sweeper_params['quad_type'] = 'LOBATTO' +sweeper_params['num_nodes'] = 4 + +# initialize problem parameters +problem_params = dict() +problem_params['mesh'] = domain_mesh +problem_params['function_space'] = V +problem_params['coupling_boundary'] = coupling_boundary +problem_params['remaining_boundary'] = remaining_boundary +problem_params['solution_expr'] = u_D +problem_params['forcing_term_expr'] = forcing_expr +problem_params['precice_ref'] = precice +problem_params['coupling_expr'] = coupling_expression +problem_params['participant_name'] = participant_name + +# initialize controller parameters +controller_params = dict() +controller_params['logger_level'] = 30 + +# fill description dictionary for easy step instantiation +description = dict() +description['problem_class'] = fenics_heat_2d +description['problem_params'] = problem_params +description['sweeper_class'] = imex_1st_order_mass +description['sweeper_params'] = sweeper_params +description['level_params'] = level_params +description['step_params'] = step_params + +# Controller for time stepping +controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + +# Reference to problem class for easy access to exact solution +P = controller.MS[0].levels[0].prob + +################################################################## + + + +################################################################# +# Coupling loop +################################################################# + +while precice.is_coupling_ongoing(): + + # write checkpoint + if precice.requires_writing_checkpoint(): + precice.store_checkpoint(u_n, t, n) + + # output solution and reference solution at t_n+1 and substeps (read from buffer) + print('output u^%d and u_ref^%d' % (n, n)) + for sample in u_write: + temperature_out << sample + + for sample in ref_write: + ref_out << sample + + for sample in error_write: + error_out << error_pointwise + + precice_dt = precice.get_max_time_step_size() + dt.assign(np.min([pySDC_dt, precice_dt])) + + # Update start time of the current preCICE time step within + # the problem class. + # This is necessary, that the coupling expression update in the problem class + # is executed correctly. + P.set_t_start(t) + + # Retrieve the result at the end of the timestep from pySDC + # Possible statistics generated by pySDC can be accessed via the 'stats' variable if needed + uend, stats = controller.run(u_n, t0=t, Tend=t + float(dt)) + + # Update the buffer with the new solution values + u_np1 = uend.values + + # 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.advance(dt) + precice_dt = precice.get_max_time_step_size() + + # roll back to checkpoint + if precice.requires_reading_checkpoint(): + u_cp, t_cp, n_cp = precice.retrieve_checkpoint() + u_n.assign(u_cp) + t = t_cp + n = n_cp + # empty buffer if window has not converged + u_write = [] + ref_write = [] + error_write = [] + else: # update solution + u_n.assign(u_np1) + t += float(dt) + n += 1 + # copy data to buffer and rename + uu = u_n.copy() + uu.rename("u", "") + u_write.append((uu, t)) + uu_ref = u_ref.copy() + uu_ref.rename("u_ref", "") + ref_write.append(uu_ref) + err = error_pointwise.copy() + err.rename("err", "") + error_write.append(err) + + 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 and substeps (read from buffer) +print("output u^%d and u_ref^%d" % (n, n)) +for sample in u_write: + temperature_out << sample + +for sample in ref_write: + ref_out << sample + +for sample in error_write: + error_out << error_pointwise + +# Hold plot +precice.finalize() \ No newline at end of file diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py new file mode 100644 index 000000000..860f8ae0e --- /dev/null +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -0,0 +1,130 @@ +from fenics import * +import dolfin as df + +from pySDC.core.Problem import ptype +from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh + + +class fenics_heat_2d(ptype): + dtype_u = fenics_mesh + dtype_f = rhs_fenics_mesh + + def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupling_boundary, remaining_boundary, coupling_expr, precice_ref, participant_name): + # Set precice reference and coupling expression reference to update coupling boundary + # at every step within pySDC + self.precice = precice_ref + self.coupling_expression = coupling_expr + self.t_start = 0.0 + self.participant_name = participant_name + + # set mesh and function space for future reference + self.mesh = mesh + self.V = function_space + + # invoke super init + super(fenics_heat_2d, self).__init__(self.V) + + # Define Trial and Test function + u = df.TrialFunction(self.V) + v = df.TestFunction(self.V) + + # Mass term + a_M = u * v * df.dx + self.M = df.assemble(a_M) + + # Stiffness term (Laplace) + a_K = -1.0 * df.inner(df.nabla_grad(u), df.nabla_grad(v)) * df.dx + self.K = df.assemble(a_K) + + # Forcing term + self.forcing_term_expr = forcing_term_expr + + # Solution expression for error comparison and as boundary condition + # on the non-coupling boundary + self.solution_expr = solution_expr + + # Currently only for Dirichlet boundary, has to be changed for Neumann boundary + if self.participant_name == 'Dirichlet': + self.couplingBC = df.DirichletBC(self.V, coupling_expr, coupling_boundary) + + self.remainingBC = df.DirichletBC(self.V, solution_expr, remaining_boundary) + + # Allow for fixing the boundary conditions for the residual computation + # Necessary if imex-1st-order-mass sweeper is used + self.fix_bc_for_residual = True + + # define the homogeneous Dirichlet boundary for residual correction + def FullBoundary(x, on_boundary): + return on_boundary + self.homogenousBC = df.DirichletBC(self.V, df.Constant(0), FullBoundary) + + + def solve_system(self, rhs, factor, u0, t): + r""" + Solves (M - factor \cdot A) \vec{u} = \vec{rhs}. + """ + u = self.dtype_u(u0) + T = self.M - factor * self.K + b = self.dtype_u(rhs) + + self.solution_expr.t = t + self.remainingBC.apply(T, b.values.vector()) + + if self.participant_name == 'Dirichlet': + # Coupling BC needs to point to correct time + dt = t - self.t_start + read_data = self.precice.read_data(dt) + self.precice.update_coupling_expression(self.coupling_expression, read_data) + + self.couplingBC.apply(T, b.values.vector()) + + df.solve(T, u.values.vector(), b.values.vector()) + + return u + + def eval_f(self, u, t): + """ + Derivative computation. + Returns a dtype_f object (rhs_fenics_mesh in our case), + with an explicit and implicit part. + """ + f = self.dtype_f(self.V) + + # Implicit part: K*u + self.K.mult(u.values.vector(), f.impl.values.vector()) + + # Explicit part: M*g (g = forcing term) + self.forcing_term_expr.t = t + f.expl = self.dtype_u(df.interpolate(self.forcing_term_expr, self.V)) + f.expl = self.apply_mass_matrix(f.expl) + + return f + + def fix_residual(self, res): + """ + Applies homogeneous Dirichlet boundary conditions to the residual + """ + self.homogenousBC.apply(res.values.vector()) + return None + + def apply_mass_matrix(self, u): + r""" + Returns M*u. + """ + + me = self.dtype_u(self.V) + self.M.mult(u.values.vector(), me.values.vector()) + + return me + + def u_exact(self, t): + r""" + Returns the exact solution at time t. + """ + self.solution_expr.t = t + + me = self.dtype_u(interpolate(self.solution_expr, self.V), val=self.V) + return me + + def set_t_start(self, t_start): + self.t_start = t_start \ No newline at end of file From 2785a59dc6ea2d4f1d1d0115fc6c813ab7edfe70 Mon Sep 17 00:00:00 2001 From: TobiasEppacher Date: Sat, 24 Aug 2024 11:29:55 +0200 Subject: [PATCH 02/13] Unify run-script --- .../dirichlet-fenics/run.sh | 15 +++++++++++++-- .../dirichlet-pySDC-fenics/clean.sh | 6 ------ .../precice-adapter-config.json | 9 --------- .../dirichlet-pySDC-fenics/run.sh | 9 --------- 4 files changed, 13 insertions(+), 26 deletions(-) delete mode 100755 partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh delete mode 100644 partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json delete mode 100755 partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index 95cc2d28b..d4149b57b 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -4,6 +4,17 @@ set -e -u . ../../tools/log.sh exec > >(tee --append "$LOGFILE") 2>&1 -python3 ../solver-fenics/heat.py Dirichlet +if [$1 = "Default"]; then + echo "Running simulation with default FEniCS implementation" + python3 ../solver-fenics/heat.py Dirichlet +else if [$1 = "HighOrder"]; then + echo "Running simulation with higher order FEniCS implementation" + python3 ../solver-fenics/heatHigherOrder.py Dirichlet +else if [$1 = "SDC"]; then + echo "Running simulation with pySDC+FEniCS implementation" + python3 ../solver-fenics/heat_pySDC.py Dirichlet +else + echo "Invalid argument. Please provide either 'Default', 'HighOrder', or 'SDC'" +fi -close_log +close_log \ No newline at end of file diff --git a/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh b/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh deleted file mode 100755 index d3f3318b4..000000000 --- a/partitioned-heat-conduction/dirichlet-pySDC-fenics/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/usr/bin/env sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_fenics . diff --git a/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json b/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.json deleted file mode 100644 index 7989ba340..000000000 --- a/partitioned-heat-conduction/dirichlet-pySDC-fenics/precice-adapter-config.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": "Heat-Flux", - "read_data_name": "Temperature" - } -} diff --git a/partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh b/partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh deleted file mode 100755 index f06eb4084..000000000 --- a/partitioned-heat-conduction/dirichlet-pySDC-fenics/run.sh +++ /dev/null @@ -1,9 +0,0 @@ -#!/usr/bin/env bash -set -e -u - -. ../../tools/log.sh -exec > >(tee --append "$LOGFILE") 2>&1 - -python3 ../solver-fenics/heat_pySDC.py Dirichlet - -close_log From dff55b993a7a901404b608e9905ec366f75161c1 Mon Sep 17 00:00:00 2001 From: TobiasEppacher Date: Wed, 28 Aug 2024 16:06:58 +0200 Subject: [PATCH 03/13] Implementation of suggestions --- .../dirichlet-fenics/run.sh | 34 +++- .../solver-fenics/heat_pySDC.py | 178 ++++++++++-------- .../solver-fenics/heat_pySDC_problem_class.py | 114 ++++++----- 3 files changed, 190 insertions(+), 136 deletions(-) diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index d4149b57b..53c2e9006 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -4,17 +4,33 @@ set -e -u . ../../tools/log.sh exec > >(tee --append "$LOGFILE") 2>&1 -if [$1 = "Default"]; then +if [ $# -eq 0 ] +then echo "Running simulation with default FEniCS implementation" python3 ../solver-fenics/heat.py Dirichlet -else if [$1 = "HighOrder"]; then - echo "Running simulation with higher order FEniCS implementation" - python3 ../solver-fenics/heatHigherOrder.py Dirichlet -else if [$1 = "SDC"]; then - echo "Running simulation with pySDC+FEniCS implementation" - python3 ../solver-fenics/heat_pySDC.py Dirichlet else - echo "Invalid argument. Please provide either 'Default', 'HighOrder', or 'SDC'" + case "$1" in + -[hH]|--help) + echo "Usage: $0 [highorder|sdc|*]" + echo "" + echo " highorder: Run simulation with higher order FEniCS implementation" + echo " sdc : Run simulation with pySDC+FEniCS implementation" + echo " * : For every other input run the simulation with default FEniCS implementation" + exit 0 + ;; + highorder) + echo "Running simulation with higher order FEniCS implementation" + python3 ../solver-fenics/heatHigherOrder.py Dirichlet + ;; + sdc) + echo "Running simulation with pySDC+FEniCS implementation" + python3 ../solver-fenics/heat_pySDC.py Dirichlet + ;; + *) + echo "Running simulation with default FEniCS implementation" + python3 ../solver-fenics/heat.py Dirichlet + ;; + esac fi -close_log \ No newline at end of file +close_log diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py index 5163fc335..6c00242f5 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py @@ -39,6 +39,7 @@ from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from heat_pySDC_problem_class import fenics_heat_2d + 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 @@ -55,6 +56,82 @@ def determine_gradient(V_g, u, flux): solve(a == L, flux) +def setup_problem( + participant_name, + domain_mesh, + coupling_boundary, + remaining_boundary, + V, + u_D, + forcing_expr, + precice, + coupling_expression, + dt, + res_tol=1e-11, + maxiter=40, + quad_type='LOBATTO', + num_nodes=4, + logger_level=30): + """ + Setup the problem for pySDC controller + :param participant_name: name of the participant + :param domain_mesh: mesh of the domain + :param coupling_boundary: boundary where coupling happens + :param remaining_boundary: remaining boundary + :param V: function space + :param u_D: initial condition + :param forcing_expr: forcing term + :param precice: preCICE adapter + :param coupling_expression: coupling expression + :param dt: pySDC time step size + :return: controller, problem class + + pySDC can be installed via `pip install pySDC`. + If you want to use the developer version, the source code repository can be cloned from "https://github.com/Parallel-in-Time/pySDC". + For more information on pySDC, see also "https://parallel-in-time.org/pySDC/" + """ + + # initialize controller parameters + controller_params = { + 'logger_level': logger_level + } + + # fill description dictionary for easy instantiation + description = { + 'problem_class': fenics_heat_2d, + 'problem_params': { + 'mesh': domain_mesh, + 'function_space': V, + 'coupling_boundary': coupling_boundary, + 'remaining_boundary': remaining_boundary, + 'solution_expr': u_D, + 'forcing_term_expr': forcing_expr, + 'precice_ref': precice, + 'coupling_expr': coupling_expression, + 'participant_name': participant_name + }, + 'sweeper_class': imex_1st_order_mass, + 'sweeper_params': { + 'quad_type': quad_type, + 'num_nodes': num_nodes + }, + 'level_params': { + 'restol': res_tol, + 'dt': dt + }, + 'step_params': { + 'maxiter': maxiter, + } + } + + # Controller for time stepping + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + # Reference to problem class for easy access to exact solution + P = controller.MS[0].levels[0].prob + + return controller, P + parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in ProblemType]) @@ -84,15 +161,9 @@ def determine_gradient(V_g, u, flux): W = V_g.sub(0).collapse() -################################################################## -# Problem definition and preCICE initialization -################################################################## - # Time step size # Should be integer fraction of the used time window size -precice_time_window_size = 0.1 -precice_waveform_deg = 1 -pySDC_dt = precice_time_window_size/precice_waveform_deg +pySDC_dt = 0.1 # Manufactured solution parameters alpha = 3 # parameter alpha @@ -138,6 +209,25 @@ def determine_gradient(V_g, u, flux): coupling_expression = precice.create_coupling_expression() +controller, P = setup_problem( + participant_name=precice.get_participant_name(), + domain_mesh=domain_mesh, + coupling_boundary=coupling_boundary, + remaining_boundary=remaining_boundary, + V=V, + u_D=u_D, + forcing_expr=forcing_expr, + precice=precice, + coupling_expression=coupling_expression, + dt=dt, + res_tol=error_tol, + maxiter=40, + quad_type='LOBATTO', + num_nodes=4, + logger_level=30 +) + + # Time-stepping u_np1 = Function(V) u_np1.rename("Temperature", "") @@ -185,74 +275,10 @@ def determine_gradient(V_g, u, flux): err.rename("err", "") error_write.append(err) - if problem is ProblemType.DIRICHLET: flux = Function(V_g) flux.rename("Heat-Flux", "") - -################################################################## -# Setup pySDC controller -################################################################## -''' -pySDC can be installed via `pip install pySDC`. -If you want to use the developer version, the source code repository can be cloned from "https://github.com/Parallel-in-Time/pySDC". -For more information on pySDC, see also "https://parallel-in-time.org/pySDC/" -''' - -# initialize level parameters -level_params = dict() -level_params['restol'] = 1e-11 -level_params['dt'] = pySDC_dt - -# initialize step parameters -step_params = dict() -step_params['maxiter'] = 40 - -# initialize sweeper parameters -sweeper_params = dict() -sweeper_params['quad_type'] = 'LOBATTO' -sweeper_params['num_nodes'] = 4 - -# initialize problem parameters -problem_params = dict() -problem_params['mesh'] = domain_mesh -problem_params['function_space'] = V -problem_params['coupling_boundary'] = coupling_boundary -problem_params['remaining_boundary'] = remaining_boundary -problem_params['solution_expr'] = u_D -problem_params['forcing_term_expr'] = forcing_expr -problem_params['precice_ref'] = precice -problem_params['coupling_expr'] = coupling_expression -problem_params['participant_name'] = participant_name - -# initialize controller parameters -controller_params = dict() -controller_params['logger_level'] = 30 - -# fill description dictionary for easy step instantiation -description = dict() -description['problem_class'] = fenics_heat_2d -description['problem_params'] = problem_params -description['sweeper_class'] = imex_1st_order_mass -description['sweeper_params'] = sweeper_params -description['level_params'] = level_params -description['step_params'] = step_params - -# Controller for time stepping -controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) - -# Reference to problem class for easy access to exact solution -P = controller.MS[0].levels[0].prob - -################################################################## - - - -################################################################# -# Coupling loop -################################################################# - while precice.is_coupling_ongoing(): # write checkpoint @@ -277,12 +303,12 @@ def determine_gradient(V_g, u, flux): # the problem class. # This is necessary, that the coupling expression update in the problem class # is executed correctly. - P.set_t_start(t) + P.t_start = t # Retrieve the result at the end of the timestep from pySDC # Possible statistics generated by pySDC can be accessed via the 'stats' variable if needed uend, stats = controller.run(u_n, t0=t, Tend=t + float(dt)) - + # Update the buffer with the new solution values u_np1 = uend.values @@ -329,8 +355,8 @@ def determine_gradient(V_g, u, flux): 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 and substeps (read from buffer) print("output u^%d and u_ref^%d" % (n, n)) for sample in u_write: @@ -343,4 +369,4 @@ def determine_gradient(V_g, u, flux): error_out << error_pointwise # Hold plot -precice.finalize() \ No newline at end of file +precice.finalize() diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py index 860f8ae0e..d20d2ba07 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -1,5 +1,4 @@ -from fenics import * -import dolfin as df +from fenics import TrialFunction, TestFunction, dx, assemble, inner, nabla_grad, DirichletBC, Constant, solve, interpolate from pySDC.core.Problem import ptype from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh @@ -8,60 +7,76 @@ class fenics_heat_2d(ptype): dtype_u = fenics_mesh dtype_f = rhs_fenics_mesh - - def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupling_boundary, remaining_boundary, coupling_expr, precice_ref, participant_name): - # Set precice reference and coupling expression reference to update coupling boundary + + def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupling_boundary, + remaining_boundary, coupling_expr, precice_ref, participant_name): + # Add docstring + """ + Constructor of the 2D heat equation problem class. + + Args: + mesh: FEniCS mesh object + function_space: FEniCS function space object + forcing_term_expr: FEniCS expression for the forcing term + solution_expr: FEniCS expression for the manufactured solution + coupling_boundary: FEniCS SubDomain object for the coupling boundary + remaining_boundary: FEniCS SubDomain object for the remaining boundary + coupling_expr: FEniCS expression for the coupling boundary condition + precice_ref: preCICE-FEniCS adapter object reference + participant_name: Name of the participant (Dirichlet or Neumann) + """ + + # Set precice reference and coupling expression reference to update coupling boundary # at every step within pySDC self.precice = precice_ref self.coupling_expression = coupling_expr self.t_start = 0.0 self.participant_name = participant_name - + # set mesh and function space for future reference self.mesh = mesh self.V = function_space - + # invoke super init super(fenics_heat_2d, self).__init__(self.V) - + # Define Trial and Test function - u = df.TrialFunction(self.V) - v = df.TestFunction(self.V) - + u = TrialFunction(self.V) + v = TestFunction(self.V) + # Mass term - a_M = u * v * df.dx - self.M = df.assemble(a_M) - + a_M = u * v * dx + self.M = assemble(a_M) + # Stiffness term (Laplace) - a_K = -1.0 * df.inner(df.nabla_grad(u), df.nabla_grad(v)) * df.dx - self.K = df.assemble(a_K) - + a_K = -1.0 * inner(nabla_grad(u), nabla_grad(v)) * dx + self.K = assemble(a_K) + # Forcing term self.forcing_term_expr = forcing_term_expr - - # Solution expression for error comparison and as boundary condition + + # Solution expression for error comparison and as boundary condition # on the non-coupling boundary self.solution_expr = solution_expr - + # Currently only for Dirichlet boundary, has to be changed for Neumann boundary if self.participant_name == 'Dirichlet': - self.couplingBC = df.DirichletBC(self.V, coupling_expr, coupling_boundary) - - self.remainingBC = df.DirichletBC(self.V, solution_expr, remaining_boundary) - + self.couplingBC = DirichletBC(self.V, coupling_expr, coupling_boundary) + + self.remainingBC = DirichletBC(self.V, solution_expr, remaining_boundary) + # Allow for fixing the boundary conditions for the residual computation # Necessary if imex-1st-order-mass sweeper is used self.fix_bc_for_residual = True - + # define the homogeneous Dirichlet boundary for residual correction def FullBoundary(x, on_boundary): return on_boundary - self.homogenousBC = df.DirichletBC(self.V, df.Constant(0), FullBoundary) - - + self.homogenousBC = DirichletBC(self.V, Constant(0), FullBoundary) + def solve_system(self, rhs, factor, u0, t): - r""" - Solves (M - factor \cdot A) \vec{u} = \vec{rhs}. + """ + Solves (M - factor * A) u = rhs. """ u = self.dtype_u(u0) T = self.M - factor * self.K @@ -69,16 +84,18 @@ def solve_system(self, rhs, factor, u0, t): self.solution_expr.t = t self.remainingBC.apply(T, b.values.vector()) - - if self.participant_name == 'Dirichlet': - # Coupling BC needs to point to correct time - dt = t - self.t_start - read_data = self.precice.read_data(dt) - self.precice.update_coupling_expression(self.coupling_expression, read_data) - - self.couplingBC.apply(T, b.values.vector()) - df.solve(T, u.values.vector(), b.values.vector()) + # Update the coupling boundary condition for the Dirichlet participant + if self.participant_name == 'Dirichlet': + dt = t - self.t_start # This dt is used to read data from the current time window + read_data = self.precice.read_data(dt) # Read the data to update the coupling expression + self.precice.update_coupling_expression( + self.coupling_expression, + read_data) # Update the coupling expression + + self.couplingBC.apply(T, b.values.vector()) # Apply the coupling boundary condition + + solve(T, u.values.vector(), b.values.vector()) return u @@ -92,23 +109,22 @@ def eval_f(self, u, t): # Implicit part: K*u self.K.mult(u.values.vector(), f.impl.values.vector()) - + # Explicit part: M*g (g = forcing term) self.forcing_term_expr.t = t - f.expl = self.dtype_u(df.interpolate(self.forcing_term_expr, self.V)) + f.expl = self.dtype_u(interpolate(self.forcing_term_expr, self.V)) f.expl = self.apply_mass_matrix(f.expl) return f - + def fix_residual(self, res): """ Applies homogeneous Dirichlet boundary conditions to the residual """ self.homogenousBC.apply(res.values.vector()) - return None - + def apply_mass_matrix(self, u): - r""" + """ Returns M*u. """ @@ -118,13 +134,9 @@ def apply_mass_matrix(self, u): return me def u_exact(self, t): - r""" + """ Returns the exact solution at time t. """ self.solution_expr.t = t - - me = self.dtype_u(interpolate(self.solution_expr, self.V), val=self.V) - return me - def set_t_start(self, t_start): - self.t_start = t_start \ No newline at end of file + return self.dtype_u(interpolate(self.solution_expr, self.V), val=self.V) From 8cea45ea4e25b4c2e09d8d82a82f554b4e55ab60 Mon Sep 17 00:00:00 2001 From: TobiasEppacher Date: Thu, 29 Aug 2024 22:20:55 +0200 Subject: [PATCH 04/13] Corrections and simplifications --- .../solver-fenics/heat_pySDC.py | 111 +++++++++--------- .../solver-fenics/heat_pySDC_problem_class.py | 22 ++-- 2 files changed, 68 insertions(+), 65 deletions(-) diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py index 6c00242f5..b7b22c4da 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py @@ -22,10 +22,16 @@ u = u_0 at t = 0 u = 1 + x^2 + alpha*y^2 + \beta*t f = beta - 2 - 2*alpha + + +This variant of this tutorial example uses the open-source library pySDC for time-stepping. +pySDC can be installed via `pip install pySDC`. +If you want to use the developer version, the source code repository can be cloned from "https://github.com/Parallel-in-Time/pySDC". +For more information on pySDC, see also "https://parallel-in-time.org/pySDC/" """ from __future__ import print_function, division -from fenics import Function, FunctionSpace, Expression, Constant, DirichletBC, TrialFunction, TestFunction, \ +from fenics import Function, FunctionSpace, Expression, Constant, TrialFunction, TestFunction, \ File, solve, grad, inner, dx, interpolate, VectorFunctionSpace, MeshFunction, MPI from fenicsprecice import Adapter from errorcomputation import compute_errors @@ -57,40 +63,43 @@ def determine_gradient(V_g, u, flux): def setup_problem( - participant_name, - domain_mesh, - coupling_boundary, - remaining_boundary, - V, - u_D, - forcing_expr, - precice, - coupling_expression, - dt, - res_tol=1e-11, - maxiter=40, - quad_type='LOBATTO', - num_nodes=4, - logger_level=30): + function_space, + coupling_boundary, + remaining_boundary, + u_D, forcing_expr, + coupling_expression, + precice, + dt, + logger_level=30, + quad_type='LOBATTO', + num_nodes=4, + restol=1e-11, + maxiter=40): + + # Create docstring for this function """ - Setup the problem for pySDC controller - :param participant_name: name of the participant - :param domain_mesh: mesh of the domain - :param coupling_boundary: boundary where coupling happens - :param remaining_boundary: remaining boundary - :param V: function space - :param u_D: initial condition - :param forcing_expr: forcing term - :param precice: preCICE adapter - :param coupling_expression: coupling expression - :param dt: pySDC time step size - :return: controller, problem class - - pySDC can be installed via `pip install pySDC`. - If you want to use the developer version, the source code repository can be cloned from "https://github.com/Parallel-in-Time/pySDC". - For more information on pySDC, see also "https://parallel-in-time.org/pySDC/" + Setup the problem and controller for the heat equation problem. + + Args: + function_space: FEniCS function space object + coupling_boundary: FEniCS SubDomain object for the coupling boundary + remaining_boundary: FEniCS SubDomain object for the remaining boundary + u_D: FEniCS expression for the manufactured solution + forcing_expr: FEniCS expression for the forcing term + coupling_expression: FEniCS expression for the coupling boundary condition + precice: preCICE-FEniCS adapter object reference + dt: time step size + logger_level: logging level + quad_type: quadrature type + num_nodes: number of nodes + restol: residual tolerance + maxiter: maximum number of iterations + + Returns: + controller: pySDC controller object + P: problem object """ - + # initialize controller parameters controller_params = { 'logger_level': logger_level @@ -100,23 +109,21 @@ def setup_problem( description = { 'problem_class': fenics_heat_2d, 'problem_params': { - 'mesh': domain_mesh, - 'function_space': V, + 'function_space': function_space, 'coupling_boundary': coupling_boundary, 'remaining_boundary': remaining_boundary, 'solution_expr': u_D, 'forcing_term_expr': forcing_expr, 'precice_ref': precice, - 'coupling_expr': coupling_expression, - 'participant_name': participant_name + 'coupling_expr': coupling_expression }, 'sweeper_class': imex_1st_order_mass, 'sweeper_params': { 'quad_type': quad_type, - 'num_nodes': num_nodes + 'num_nodes': num_nodes, }, 'level_params': { - 'restol': res_tol, + 'restol': restol, 'dt': dt }, 'step_params': { @@ -129,10 +136,10 @@ def setup_problem( # Reference to problem class for easy access to exact solution P = controller.MS[0].levels[0].prob - return controller, P + parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in ProblemType]) parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-8,) @@ -209,23 +216,13 @@ def setup_problem( coupling_expression = precice.create_coupling_expression() -controller, P = setup_problem( - participant_name=precice.get_participant_name(), - domain_mesh=domain_mesh, - coupling_boundary=coupling_boundary, - remaining_boundary=remaining_boundary, - V=V, - u_D=u_D, - forcing_expr=forcing_expr, - precice=precice, - coupling_expression=coupling_expression, - dt=dt, - res_tol=error_tol, - maxiter=40, - quad_type='LOBATTO', - num_nodes=4, - logger_level=30 -) +controller, P = setup_problem(V, + coupling_boundary, + remaining_boundary, + u_D, forcing_expr, + coupling_expression, + precice, + pySDC_dt) # Time-stepping diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py index d20d2ba07..09034bccb 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -1,4 +1,14 @@ +""" +This file contains the implementation of the pySDC problem class for the 2D heat equation. +This class is used to define the problem pySDC uses in "heat_pySDC.py" for time-stepping. + +This class is based on the pySDC tutorial example for its usage with FEniCS. +More about that tutorial can be found at "https://parallel-in-time.org/pySDC/tutorial/step_7.html#part-a-pysdc-and-fenics". +""" + from fenics import TrialFunction, TestFunction, dx, assemble, inner, nabla_grad, DirichletBC, Constant, solve, interpolate +from my_enums import ProblemType + from pySDC.core.Problem import ptype from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh @@ -8,14 +18,13 @@ class fenics_heat_2d(ptype): dtype_u = fenics_mesh dtype_f = rhs_fenics_mesh - def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupling_boundary, - remaining_boundary, coupling_expr, precice_ref, participant_name): + def __init__(self, function_space, forcing_term_expr, solution_expr, coupling_boundary, + remaining_boundary, coupling_expr, precice_ref): # Add docstring """ Constructor of the 2D heat equation problem class. Args: - mesh: FEniCS mesh object function_space: FEniCS function space object forcing_term_expr: FEniCS expression for the forcing term solution_expr: FEniCS expression for the manufactured solution @@ -23,7 +32,6 @@ def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupl remaining_boundary: FEniCS SubDomain object for the remaining boundary coupling_expr: FEniCS expression for the coupling boundary condition precice_ref: preCICE-FEniCS adapter object reference - participant_name: Name of the participant (Dirichlet or Neumann) """ # Set precice reference and coupling expression reference to update coupling boundary @@ -31,10 +39,8 @@ def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupl self.precice = precice_ref self.coupling_expression = coupling_expr self.t_start = 0.0 - self.participant_name = participant_name # set mesh and function space for future reference - self.mesh = mesh self.V = function_space # invoke super init @@ -60,7 +66,7 @@ def __init__(self, mesh, function_space, forcing_term_expr, solution_expr, coupl self.solution_expr = solution_expr # Currently only for Dirichlet boundary, has to be changed for Neumann boundary - if self.participant_name == 'Dirichlet': + if self.precice.get_participant_name() == ProblemType.DIRICHLET.value: self.couplingBC = DirichletBC(self.V, coupling_expr, coupling_boundary) self.remainingBC = DirichletBC(self.V, solution_expr, remaining_boundary) @@ -86,7 +92,7 @@ def solve_system(self, rhs, factor, u0, t): self.remainingBC.apply(T, b.values.vector()) # Update the coupling boundary condition for the Dirichlet participant - if self.participant_name == 'Dirichlet': + if self.precice.get_participant_name() == ProblemType.DIRICHLET.value: dt = t - self.t_start # This dt is used to read data from the current time window read_data = self.precice.read_data(dt) # Read the data to update the coupling expression self.precice.update_coupling_expression( From a2f2fedde9be1d04197368ec2fbe1d922d4db241 Mon Sep 17 00:00:00 2001 From: TobiasEppacher Date: Thu, 29 Aug 2024 22:46:40 +0200 Subject: [PATCH 05/13] fix comment --- .../solver-fenics/heat_pySDC_problem_class.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py index 09034bccb..be6fdc345 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -40,8 +40,15 @@ def __init__(self, function_space, forcing_term_expr, solution_expr, coupling_bo self.coupling_expression = coupling_expr self.t_start = 0.0 - # set mesh and function space for future reference + # save function space for future reference self.V = function_space + + # Forcing term + self.forcing_term_expr = forcing_term_expr + + # Solution expression for error comparison and as boundary condition + # on the non-coupling boundary + self.solution_expr = solution_expr # invoke super init super(fenics_heat_2d, self).__init__(self.V) @@ -58,14 +65,7 @@ def __init__(self, function_space, forcing_term_expr, solution_expr, coupling_bo a_K = -1.0 * inner(nabla_grad(u), nabla_grad(v)) * dx self.K = assemble(a_K) - # Forcing term - self.forcing_term_expr = forcing_term_expr - - # Solution expression for error comparison and as boundary condition - # on the non-coupling boundary - self.solution_expr = solution_expr - - # Currently only for Dirichlet boundary, has to be changed for Neumann boundary + # Currently only Dirichlet participant is supported if self.precice.get_participant_name() == ProblemType.DIRICHLET.value: self.couplingBC = DirichletBC(self.V, coupling_expr, coupling_boundary) From 3a00c01c7d18af150c90a4c5eee382570cb96135 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 30 Aug 2024 16:18:20 +0200 Subject: [PATCH 06/13] Run case in venv. See #547. --- partitioned-heat-conduction/dirichlet-fenics/run.sh | 4 ++++ partitioned-heat-conduction/neumann-fenics/run.sh | 4 ++++ .../solver-fenics/requirements.txt | 11 +++++++++++ 3 files changed, 19 insertions(+) create mode 100644 partitioned-heat-conduction/solver-fenics/requirements.txt diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index 53c2e9006..ccf3b10f1 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -1,6 +1,10 @@ #!/usr/bin/env bash set -e -u +python3 -m venv --system-site-packages ../solver-fenics/.venv +. ../solver-fenics/.venv/bin/activate +pip install -r ../solver-fenics/requirements.txt + . ../../tools/log.sh exec > >(tee --append "$LOGFILE") 2>&1 diff --git a/partitioned-heat-conduction/neumann-fenics/run.sh b/partitioned-heat-conduction/neumann-fenics/run.sh index 874fdec84..ef1398d56 100755 --- a/partitioned-heat-conduction/neumann-fenics/run.sh +++ b/partitioned-heat-conduction/neumann-fenics/run.sh @@ -1,6 +1,10 @@ #!/usr/bin/env bash set -e -u +python3 -m venv --system-site-packages ../solver-fenics/.venv +. ../solver-fenics/.venv/bin/activate +pip install -r ../solver-fenics/requirements.txt + . ../../tools/log.sh exec > >(tee --append "$LOGFILE") 2>&1 diff --git a/partitioned-heat-conduction/solver-fenics/requirements.txt b/partitioned-heat-conduction/solver-fenics/requirements.txt new file mode 100644 index 000000000..fe3302591 --- /dev/null +++ b/partitioned-heat-conduction/solver-fenics/requirements.txt @@ -0,0 +1,11 @@ +numpy >1, <2 +fenicsprecice~=2.0 +scipy + +# Assuming FEniCS from ppa:fenics-packages/fenics was installed https://fenicsproject.org/download/archive/ +# Use --system-site-packages in venv +fenics-dijitso==2019.2.0.dev0 +fenics-dolfin==2019.2.0.13.dev0 +fenics-ffc==2019.2.0.dev0 +fenics-fiat==2019.2.0.dev0 +fenics-ufl-legacy==2022.3.0 \ No newline at end of file From 15f6d37d6bcd729cc849ad3d1e8e73cad0e3906e Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 30 Aug 2024 16:21:10 +0200 Subject: [PATCH 07/13] Fix autopep8. --- .../solver-fenics/heat_pySDC.py | 43 +++++++++---------- .../solver-fenics/heat_pySDC_problem_class.py | 2 +- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py index b7b22c4da..731989f37 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py @@ -22,7 +22,7 @@ u = u_0 at t = 0 u = 1 + x^2 + alpha*y^2 + \beta*t f = beta - 2 - 2*alpha - + This variant of this tutorial example uses the open-source library pySDC for time-stepping. pySDC can be installed via `pip install pySDC`. @@ -63,23 +63,23 @@ def determine_gradient(V_g, u, flux): def setup_problem( - function_space, - coupling_boundary, - remaining_boundary, - u_D, forcing_expr, - coupling_expression, - precice, - dt, + function_space, + coupling_boundary, + remaining_boundary, + u_D, forcing_expr, + coupling_expression, + precice, + dt, logger_level=30, - quad_type='LOBATTO', - num_nodes=4, - restol=1e-11, + quad_type='LOBATTO', + num_nodes=4, + restol=1e-11, maxiter=40): - + # Create docstring for this function """ Setup the problem and controller for the heat equation problem. - + Args: function_space: FEniCS function space object coupling_boundary: FEniCS SubDomain object for the coupling boundary @@ -94,12 +94,12 @@ def setup_problem( num_nodes: number of nodes restol: residual tolerance maxiter: maximum number of iterations - + Returns: controller: pySDC controller object P: problem object """ - + # initialize controller parameters controller_params = { 'logger_level': logger_level @@ -139,7 +139,6 @@ def setup_problem( return controller, P - parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in ProblemType]) parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-8,) @@ -216,12 +215,12 @@ def setup_problem( coupling_expression = precice.create_coupling_expression() -controller, P = setup_problem(V, - coupling_boundary, - remaining_boundary, - u_D, forcing_expr, - coupling_expression, - precice, +controller, P = setup_problem(V, + coupling_boundary, + remaining_boundary, + u_D, forcing_expr, + coupling_expression, + precice, pySDC_dt) diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py index be6fdc345..ab952c20c 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -42,7 +42,7 @@ def __init__(self, function_space, forcing_term_expr, solution_expr, coupling_bo # save function space for future reference self.V = function_space - + # Forcing term self.forcing_term_expr = forcing_term_expr From f8ca67f46d6895e9c2deb8740ffcd5348f67e09b Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 30 Aug 2024 16:31:56 +0200 Subject: [PATCH 08/13] Rename option --- partitioned-heat-conduction/dirichlet-fenics/run.sh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index ccf3b10f1..628410cab 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -17,13 +17,13 @@ else -[hH]|--help) echo "Usage: $0 [highorder|sdc|*]" echo "" - echo " highorder: Run simulation with higher order FEniCS implementation" - echo " sdc : Run simulation with pySDC+FEniCS implementation" - echo " * : For every other input run the simulation with default FEniCS implementation" + echo " irk : Run simulation with higher-order implicit Runge-Kutta schemes FEniCS implementation" + echo " sdc : Run simulation with pySDC+FEniCS implementation" + echo " * : For every other input run the simulation with default FEniCS implementation" exit 0 ;; - highorder) - echo "Running simulation with higher order FEniCS implementation" + irk) + echo "Running simulation with higher-order implicit Runge-Kutta schemes FEniCS implementation" python3 ../solver-fenics/heatHigherOrder.py Dirichlet ;; sdc) From d6a950f96733f0e908c149946942d384334fb0fc Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 30 Aug 2024 19:03:13 +0200 Subject: [PATCH 09/13] Add checks and instructions for pySDC installation. --- oscillator/mass-right-python/requirements.txt | 2 +- .../dirichlet-fenics/run.sh | 1 + .../solver-fenics/check_pySDC.py | 24 +++++++++++++++++++ .../solver-fenics/requirements.txt | 1 + 4 files changed, 27 insertions(+), 1 deletion(-) create mode 100644 partitioned-heat-conduction/solver-fenics/check_pySDC.py diff --git a/oscillator/mass-right-python/requirements.txt b/oscillator/mass-right-python/requirements.txt index 9dff131c5..e94eeb2d1 100644 --- a/oscillator/mass-right-python/requirements.txt +++ b/oscillator/mass-right-python/requirements.txt @@ -1,3 +1,3 @@ numpy >1, <2 pyprecice~=3.0 -scipy +scipy \ No newline at end of file diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index 628410cab..441737c05 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -28,6 +28,7 @@ else ;; sdc) echo "Running simulation with pySDC+FEniCS implementation" + python3 ../solver-fenics/check_pySDC.py python3 ../solver-fenics/heat_pySDC.py Dirichlet ;; *) diff --git a/partitioned-heat-conduction/solver-fenics/check_pySDC.py b/partitioned-heat-conduction/solver-fenics/check_pySDC.py new file mode 100644 index 000000000..dd35c6749 --- /dev/null +++ b/partitioned-heat-conduction/solver-fenics/check_pySDC.py @@ -0,0 +1,24 @@ +""" +a helper script to test the installation of pySDC and provide advice if an error is detected +""" +import sys +try: + from pySDC.core.Problem import ptype +except ModuleNotFoundError: + print(""" +ERROR: + +Something is wrong with your installation of pySDC. Please make sure that you have cloned the GitHub repository of pySDC and you are using the correct version. Use the following command to clone pySDC version 5.5.0: + +git clone --branch 5.5.0 https://github.com/Parallel-in-Time/pySDC <> + +You then have to add the folder into which you cloned pySDC to your PYTHONPATH: + +export PYTHONPATH=$PYTHONPATH:<> + +If you are running this tutorial in an venv (e.g. using the run script), please make sure that you also correctly set the PYTHONPATH in your venv. + +Refer to https://parallel-in-time.org/pySDC/ for installation instructions. +""") + sys.exit(1) + diff --git a/partitioned-heat-conduction/solver-fenics/requirements.txt b/partitioned-heat-conduction/solver-fenics/requirements.txt index fe3302591..68fcd337b 100644 --- a/partitioned-heat-conduction/solver-fenics/requirements.txt +++ b/partitioned-heat-conduction/solver-fenics/requirements.txt @@ -1,6 +1,7 @@ numpy >1, <2 fenicsprecice~=2.0 scipy +pySDC~=5.5 # Assuming FEniCS from ppa:fenics-packages/fenics was installed https://fenicsproject.org/download/archive/ # Use --system-site-packages in venv From b4bca817fc0758b4540626267b23d07c506fb012 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 2 Sep 2024 10:50:57 +0200 Subject: [PATCH 10/13] Install pySDC in venv using pip install git+...@5.4.3. --- .../dirichlet-fenics/run.sh | 4 +++- .../solver-fenics/check_pySDC.py | 24 ------------------- .../solver-fenics/requirements.txt | 1 - 3 files changed, 3 insertions(+), 26 deletions(-) delete mode 100644 partitioned-heat-conduction/solver-fenics/check_pySDC.py diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index 441737c05..b75e8c84b 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -27,8 +27,10 @@ else python3 ../solver-fenics/heatHigherOrder.py Dirichlet ;; sdc) + # install pySDC + its dependencies only if needed + pip install git+https://github.com/Parallel-in-Time/pySDC@5.4.3 + pip install pySDC~=5.4 echo "Running simulation with pySDC+FEniCS implementation" - python3 ../solver-fenics/check_pySDC.py python3 ../solver-fenics/heat_pySDC.py Dirichlet ;; *) diff --git a/partitioned-heat-conduction/solver-fenics/check_pySDC.py b/partitioned-heat-conduction/solver-fenics/check_pySDC.py deleted file mode 100644 index dd35c6749..000000000 --- a/partitioned-heat-conduction/solver-fenics/check_pySDC.py +++ /dev/null @@ -1,24 +0,0 @@ -""" -a helper script to test the installation of pySDC and provide advice if an error is detected -""" -import sys -try: - from pySDC.core.Problem import ptype -except ModuleNotFoundError: - print(""" -ERROR: - -Something is wrong with your installation of pySDC. Please make sure that you have cloned the GitHub repository of pySDC and you are using the correct version. Use the following command to clone pySDC version 5.5.0: - -git clone --branch 5.5.0 https://github.com/Parallel-in-Time/pySDC <> - -You then have to add the folder into which you cloned pySDC to your PYTHONPATH: - -export PYTHONPATH=$PYTHONPATH:<> - -If you are running this tutorial in an venv (e.g. using the run script), please make sure that you also correctly set the PYTHONPATH in your venv. - -Refer to https://parallel-in-time.org/pySDC/ for installation instructions. -""") - sys.exit(1) - diff --git a/partitioned-heat-conduction/solver-fenics/requirements.txt b/partitioned-heat-conduction/solver-fenics/requirements.txt index 68fcd337b..fe3302591 100644 --- a/partitioned-heat-conduction/solver-fenics/requirements.txt +++ b/partitioned-heat-conduction/solver-fenics/requirements.txt @@ -1,7 +1,6 @@ numpy >1, <2 fenicsprecice~=2.0 scipy -pySDC~=5.5 # Assuming FEniCS from ppa:fenics-packages/fenics was installed https://fenicsproject.org/download/archive/ # Use --system-site-packages in venv From 049ee9fb84cfa3e3a911417234c7de4a821d95c1 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 2 Sep 2024 11:02:00 +0200 Subject: [PATCH 11/13] Required changes for pySDC5.5.0 compatibility. --- partitioned-heat-conduction/dirichlet-fenics/run.sh | 4 ++-- .../solver-fenics/heat_pySDC_problem_class.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/partitioned-heat-conduction/dirichlet-fenics/run.sh b/partitioned-heat-conduction/dirichlet-fenics/run.sh index b75e8c84b..1f7bcb9c0 100755 --- a/partitioned-heat-conduction/dirichlet-fenics/run.sh +++ b/partitioned-heat-conduction/dirichlet-fenics/run.sh @@ -28,8 +28,8 @@ else ;; sdc) # install pySDC + its dependencies only if needed - pip install git+https://github.com/Parallel-in-Time/pySDC@5.4.3 - pip install pySDC~=5.4 + pip install git+https://github.com/Parallel-in-Time/pySDC@5.5.0 + pip install pySDC~=5.5 echo "Running simulation with pySDC+FEniCS implementation" python3 ../solver-fenics/heat_pySDC.py Dirichlet ;; diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py index ab952c20c..7adb896ed 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC_problem_class.py @@ -10,11 +10,11 @@ from my_enums import ProblemType -from pySDC.core.Problem import ptype +from pySDC.core.problem import Problem from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh -class fenics_heat_2d(ptype): +class fenics_heat_2d(Problem): dtype_u = fenics_mesh dtype_f = rhs_fenics_mesh From e889225ca17bee24f04072239004552831f7068d Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 2 Sep 2024 11:02:27 +0200 Subject: [PATCH 12/13] Remove unused. --- partitioned-heat-conduction/solver-fenics/heat_pySDC.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py index 731989f37..e0bad82ac 100644 --- a/partitioned-heat-conduction/solver-fenics/heat_pySDC.py +++ b/partitioned-heat-conduction/solver-fenics/heat_pySDC.py @@ -302,8 +302,7 @@ def setup_problem( P.t_start = t # Retrieve the result at the end of the timestep from pySDC - # Possible statistics generated by pySDC can be accessed via the 'stats' variable if needed - uend, stats = controller.run(u_n, t0=t, Tend=t + float(dt)) + uend, _ = controller.run(u_n, t0=t, Tend=t + float(dt)) # Update the buffer with the new solution values u_np1 = uend.values From e36319d407866a41a12f54d6cea831b8f793ab2c Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Mon, 2 Sep 2024 11:56:32 +0200 Subject: [PATCH 13/13] Remove unrelated change. --- oscillator/mass-right-python/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oscillator/mass-right-python/requirements.txt b/oscillator/mass-right-python/requirements.txt index e94eeb2d1..9dff131c5 100644 --- a/oscillator/mass-right-python/requirements.txt +++ b/oscillator/mass-right-python/requirements.txt @@ -1,3 +1,3 @@ numpy >1, <2 pyprecice~=3.0 -scipy \ No newline at end of file +scipy