Source code for pymordemos.hapod

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

from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor
import time

import numpy as np
from typer import Argument, Option, run

from pymor.analyticalproblems.burgers import burgers_problem_2d
from pymor.discretizers.builtin import discretize_instationary_fv, RectGrid
from pymor.algorithms.hapod import dist_vectorarray_hapod, inc_vectorarray_hapod
from pymor.algorithms.pod import pod
from pymor.tools.table import format_table


[docs]def main( tol: float = Argument(..., help='Prescribed mean l2 approximation error.'), dist: int = Argument(..., help='Number of slices for distributed HAPOD.'), inc: int = Argument(..., help='Number of steps for incremental HAPOD.'), grid: int = Option(60, help='Use grid with (2*NI)*NI elements.'), nt: int = Option(100, help='Number of time steps.'), omega: float = Option(0.9, help='Parameter omega from HAPOD algorithm.'), procs: int = Option(0, help='Number of processes to use for parallelization.'), snap: int = Option(20, help='Number of snapshot trajectories to compute.'), threads: int = Option(0, help='Number of threads to use for parallelization.'), ): """Compression of snapshot data with the HAPOD algorithm from [HLR18].""" assert procs == 0 or threads == 0 executor = ProcessPoolExecutor(procs) if procs > 0 else \ ThreadPoolExecutor(threads) if threads > 0 else \ None p = burgers_problem_2d() m, data = discretize_instationary_fv(p, grid_type=RectGrid, diameter=np.sqrt(2)/grid, nt=nt) U = m.solution_space.empty() for mu in p.parameter_space.sample_randomly(snap): U.append(m.solve(mu)) tic = time.perf_counter() pod_modes = pod(U, l2_err=tol * np.sqrt(len(U)), product=m.l2_product)[0] pod_time = time.perf_counter() - tic tic = time.perf_counter() dist_modes = dist_vectorarray_hapod(dist, U, tol, omega, product=m.l2_product, executor=executor)[0] dist_time = time.perf_counter() - tic tic = time.perf_counter() inc_modes = inc_vectorarray_hapod(inc, U, tol, omega, product=m.l2_product)[0] inc_time = time.perf_counter() - tic print(f'Snapshot matrix: {U.dim} x {len(U)}') print(format_table([ ['Method', 'Error', 'Modes', 'Time'], ['POD', np.linalg.norm(m.l2_norm(U-pod_modes.lincomb(m.l2_product.apply2(U, pod_modes)))/np.sqrt(len(U))), len(pod_modes), pod_time], ['DIST HAPOD', np.linalg.norm(m.l2_norm(U-dist_modes.lincomb(m.l2_product.apply2(U, dist_modes)))/np.sqrt(len(U))), len(dist_modes), dist_time], ['INC HAPOD', np.linalg.norm(m.l2_norm(U-inc_modes.lincomb(m.l2_product.apply2(U, inc_modes)))/np.sqrt(len(U))), len(inc_modes), inc_time]] ))
if __name__ == '__main__': run(main)