Source code for pymordemos.fenics_nonlinear

#!/usr/bin/env python
# This file is part of the pyMOR project (
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (

"""Simple example script for reducing a FEniCS-based nonlinear diffusion problem.


    DIM               Spatial dimension of the problem.
    N                 Number of mesh intervals per spatial dimension.
    ORDER             Finite element order.

    -h, --help   Show this message.

from docopt import docopt

[docs]def discretize(DIM, N, ORDER): # ### problem definition import dolfin as df if DIM == 2: mesh = df.UnitSquareMesh(N, N) elif DIM == 3: mesh = df.UnitCubeMesh(N, N, N) else: raise NotImplementedError V = df.FunctionSpace(mesh, "CG", ORDER) g = df.Constant(1.0) c = df.Constant(1.) class DirichletBoundary(df.SubDomain): def inside(self, x, on_boundary): return abs(x[0] - 1.0) < df.DOLFIN_EPS and on_boundary db = DirichletBoundary() bc = df.DirichletBC(V, g, db) u = df.Function(V) v = df.TestFunction(V) f = df.Expression("x[0]*sin(x[1])", degree=2) F = df.inner((1 + c*u**2)*df.grad(u), df.grad(v))*df.dx - f*v*df.dx df.solve(F == 0, u, bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}}) # ### pyMOR wrapping from pymor.bindings.fenics import FenicsVectorSpace, FenicsOperator, FenicsVisualizer from pymor.models.basic import StationaryModel from pymor.operators.constructions import VectorOperator space = FenicsVectorSpace(V) op = FenicsOperator(F, space, space, u, (bc,), parameter_setter=lambda mu: c.assign(mu['c'].item()), parameters={'c': 1}, solver_options={'inverse': {'type': 'newton', 'rtol': 1e-6}}) rhs = VectorOperator(op.range.zeros()) fom = StationaryModel(op, rhs, visualizer=FenicsVisualizer(space)) return fom
[docs]def fenics_nonlinear_demo(args): DIM = int(args['DIM']) N = int(args['N']) ORDER = int(args['ORDER']) from import mpi if mpi.parallel: from pymor.models.mpi import mpi_wrap_model local_models =, discretize, DIM, N, ORDER) fom = mpi_wrap_model(local_models, use_with=True, pickle_local_spaces=False) else: fom = discretize(DIM, N, ORDER) parameter_space =, 1000.)) # ### ROM generation (POD/DEIM) from pymor.algorithms.ei import ei_greedy from pymor.algorithms.newton import newton from pymor.algorithms.pod import pod from pymor.operators.ei import EmpiricalInterpolatedOperator from pymor.reductors.basic import StationaryRBReductor U = fom.solution_space.empty() residuals = fom.solution_space.empty() for mu in parameter_space.sample_uniformly(10): UU, data = newton(fom.operator, fom.rhs.as_vector(), mu=mu, rtol=1e-6, return_residuals=True) U.append(UU) residuals.append(data['residuals']) dofs, cb, _ = ei_greedy(residuals, rtol=1e-7) ei_op = EmpiricalInterpolatedOperator(fom.operator, collateral_basis=cb, interpolation_dofs=dofs, triangular=True) rb, svals = pod(U, rtol=1e-7) fom_ei = fom.with_(operator=ei_op) reductor = StationaryRBReductor(fom_ei, rb) rom = reductor.reduce() # the reductor currently removes all solver_options so we need to add them again rom = rom.with_(operator=rom.operator.with_(solver_options=fom.operator.solver_options)) # ### ROM validation import time import numpy as np # ensure that FFC is not called during runtime measurements rom.solve(1) errs = [] speedups = [] for mu in parameter_space.sample_randomly(10): tic = time.time() U = fom.solve(mu) t_fom = time.time() - tic tic = time.time() u_red = rom.solve(mu) t_rom = time.time() - tic U_red = reductor.reconstruct(u_red) errs.append(((U - U_red).l2_norm() / U.l2_norm())[0]) speedups.append(t_fom / t_rom) print(f'Maximum relative ROM error: {max(errs)}') print(f'Median of ROM speedup: {np.median(speedups)}')
if __name__ == '__main__': args = docopt(__doc__) fenics_nonlinear_demo(args)