Source code for pymordemos.parabolic

#!/usr/bin/env python
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

"""Simple demonstration of solving parabolic equations using pyMOR's builtin discretization toolkit.

Usage:
    parabolic.py [options] heat TOP
    parabolic.py [options] dar SPEED

Arguments:
    TOP          The heat diffusion coefficient for the top bars.
    SPEED        The advection speed.

Options:
    -h, --help   Show this message.
    --fv         Use finite volume discretization instead of finite elements.
    --rect       Use RectGrid instead of TriaGrid.
    --grid=NI    Use grid with NIxNI intervals [default: 100].
    --nt=COUNT   Number of time steps [default: 100].
"""

from docopt import docopt
import numpy as np

from pymor.basic import *


[docs]def parabolic_demo(args): args['--nt'] = int(args['--nt']) args['--grid'] = int(args['--grid']) if args['heat']: args['TOP'] = float(args['TOP']) problem = InstationaryProblem( StationaryProblem( domain=RectDomain(top='dirichlet', bottom='neumann'), diffusion=LincombFunction( [ConstantFunction(1., dim_domain=2), ExpressionFunction('(x[..., 0] > 0.45) * (x[..., 0] < 0.55) * (x[..., 1] < 0.7) * 1.', dim_domain=2), ExpressionFunction('(x[..., 0] > 0.35) * (x[..., 0] < 0.40) * (x[..., 1] > 0.3) * 1. + ' '(x[..., 0] > 0.60) * (x[..., 0] < 0.65) * (x[..., 1] > 0.3) * 1.', dim_domain=2)], [1., 100. - 1., ExpressionParameterFunctional('top - 1.', {'top': 0})] ), rhs=ConstantFunction(value=0., dim_domain=2), dirichlet_data=ConstantFunction(value=0., dim_domain=2), neumann_data=ExpressionFunction('(x[..., 0] > 0.45) * (x[..., 0] < 0.55) * -1000.', dim_domain=2), ), T=1., initial_data=ExpressionFunction('(x[..., 0] > 0.45) * (x[..., 0] < 0.55) * (x[..., 1] < 0.7) * 10.', dim_domain=2), parameter_space=CubicParameterSpace({'top': 0}, minimum=1, maximum=100.) ) else: args['SPEED'] = float(args['SPEED']) problem = InstationaryProblem( StationaryProblem( domain=RectDomain(), diffusion=ConstantFunction(0.01, dim_domain=2), advection=LincombFunction([ConstantFunction(np.array([-1., 0]), dim_domain=2)], [ProjectionParameterFunctional('speed', ())]), reaction=ConstantFunction(0.5, dim_domain=2), rhs=ExpressionFunction('(x[..., 0] > 0.3) * (x[..., 0] < 0.7) * (x[..., 1] > 0.3)*(x[...,1]<0.7) * 0.', dim_domain=2), dirichlet_data=ConstantFunction(value=0., dim_domain=2), ), T=1., initial_data=ExpressionFunction('(x[..., 0] > 0.3) * (x[..., 0] < 0.7) * (x[...,1]>0.3) * (x[..., 1] < 0.7) * 10.', dim_domain=2), ) print('Discretize ...') discretizer = discretize_instationary_fv if args['--fv'] else discretize_instationary_cg m, data = discretizer( analytical_problem=problem, grid_type=RectGrid if args['--rect'] else TriaGrid, diameter=np.sqrt(2) / args['--grid'] if args['--rect'] else 1. / args['--grid'], nt=args['--nt'] ) grid = data['grid'] print(grid) print() print('Solve ...') U = m.solve({'top': args['TOP']} if args['heat'] else {'speed': args['SPEED']}) m.visualize(U, title='Solution') print('')
if __name__ == '__main__': args = docopt(__doc__) parabolic_demo(args)