Source code for pymordemos.parametric_heat

#!/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 (

import numpy as np
import matplotlib.pyplot as plt
from typer import Argument, run

from pymor.analyticalproblems.domaindescriptions import LineDomain
from pymor.analyticalproblems.elliptic import StationaryProblem
from pymor.analyticalproblems.functions import ConstantFunction, ExpressionFunction, LincombFunction
from pymor.analyticalproblems.instationary import InstationaryProblem
from pymor.core.config import config
from pymor.core.logger import set_log_levels
from pymor.discretizers.builtin import discretize_instationary_cg
from pymor.parameters.functionals import ProjectionParameterFunctional
from import BTReductor, LQGBTReductor, BRBTReductor
from pymor.reductors.h2 import IRKAReductor, TSIAReductor, OneSidedIRKAReductor

[docs]def run_mor_method_param(fom, r, w, mus, reductor_cls, reductor_short_name, **reductor_kwargs): """Plot reductor errors for different parameter values. Parameters ---------- fom The full-order |LTIModel|. r The order of the reduced-order model. w Array of frequencies. mus An array of parameter values. reductor_cls The reductor class. reductor_short_name A short name for the reductor. reductor_kwargs Optional keyword arguments for the reductor class. """ # Reduction roms = [] for mu in mus: rom = reductor_cls(fom, mu=mu, **reductor_kwargs).reduce(r) roms.append(rom) # Poles fig, ax = plt.subplots() for rom in roms: poles_rom = rom.poles() ax.plot(poles_rom.real, poles_rom.imag, '.', label=fr'$\mu = {mu}$') ax.set_title(f"{reductor_short_name} reduced model's poles") # Magnitude plots fig, ax = plt.subplots() for mu, rom in zip(mus, roms): rom.mag_plot(w, ax=ax, label=fr'$\mu = {mu}$') ax.set_title('Magnitude plot of {reductor_short_name} reduced models') ax.legend() fig, ax = plt.subplots() for mu, rom in zip(mus, roms): (fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$') ax.set_title('Magnitude plot of the {reductor_short_name} error system') ax.legend() # Errors for mu, rom in zip(mus, roms): err = fom - rom print(f'mu = {mu}') print(f' {reductor_short_name} relative H_2-error:' f' {err.h2_norm(mu=mu) / fom.h2_norm(mu=mu):e}') if config.HAVE_SLYCOT: print(f' {reductor_short_name} relative H_inf-error:' f' {err.hinf_norm(mu=mu) / fom.hinf_norm(mu=mu):e}') print(f' {reductor_short_name} relative Hankel-error:' f' {err.hankel_norm(mu=mu) / fom.hankel_norm(mu=mu):e}')
[docs]def main( diameter: float = Argument(0.01, help='Diameter option for the domain discretizer.'), r: int = Argument(5, help='Order of the ROMs.'), ): """Parametric 1D heat equation example.""" set_log_levels({'pymor.algorithms.gram_schmidt.gram_schmidt': 'WARNING'}) # Model p = InstationaryProblem( StationaryProblem( domain=LineDomain([0., 1.], left='robin', right='robin'), diffusion=LincombFunction([ExpressionFunction('(x[...,0] <= 0.5) * 1.', 1), ExpressionFunction('(0.5 < x[...,0]) * 1.', 1)], [1, ProjectionParameterFunctional('diffusion')]), robin_data=(ConstantFunction(1., 1), ExpressionFunction('(x[...,0] < 1e-10) * 1.', 1)), outputs=(('l2_boundary', ExpressionFunction('(x[...,0] > (1 - 1e-10)) * 1.', 1)),), ), ConstantFunction(0., 1), T=3. ) fom, _ = discretize_instationary_cg(p, diameter=diameter, nt=100) fom.visualize(fom.solve(mu=0.1)) fom.visualize(fom.solve(mu=1)) fom.visualize(fom.solve(mu=10)) lti = fom.to_lti() print(f'order of the model = {lti.order}') print(f'number of inputs = {lti.dim_input}') print(f'number of outputs = {lti.dim_output}') mu_list = [0.1, 1, 10] w_list = np.logspace(-1, 3, 100) # System poles fig, ax = plt.subplots() for mu in mu_list: poles = lti.poles(mu=mu) ax.plot(poles.real, poles.imag, '.', label=fr'$\mu = {mu}$') ax.set_title('System poles') ax.legend() # Magnitude plots fig, ax = plt.subplots() for mu in mu_list: lti.mag_plot(w_list, ax=ax, mu=mu, label=fr'$\mu = {mu}$') ax.set_title('Magnitude plot of the full model') ax.legend() # Hankel singular values fig, ax = plt.subplots() for mu in mu_list: hsv = lti.hsv(mu=mu) ax.semilogy(range(1, len(hsv) + 1), hsv, label=fr'$\mu = {mu}$') ax.set_title('Hankel singular values') ax.legend() # System norms for mu in mu_list: print(f'mu = {mu}:') print(f' H_2-norm of the full model: {lti.h2_norm(mu=mu):e}') if config.HAVE_SLYCOT: print(f' H_inf-norm of the full model: {lti.hinf_norm(mu=mu):e}') print(f' Hankel-norm of the full model: {lti.hankel_norm(mu=mu):e}') # Model order reduction run_mor_method_param(lti, r, w_list, mu_list, BTReductor, 'BT') run_mor_method_param(lti, r, w_list, mu_list, LQGBTReductor, 'LQGBT') run_mor_method_param(lti, r, w_list, mu_list, BRBTReductor, 'BRBT') run_mor_method_param(lti, r, w_list, mu_list, IRKAReductor, 'IRKA') run_mor_method_param(lti, r, w_list, mu_list, TSIAReductor, 'TSIA') run_mor_method_param(lti, r, w_list, mu_list, OneSidedIRKAReductor, 'OS-IRKA', version='V')
if __name__ == "__main__": run(main)