From 72c8cdf781dcb4716a70a720275ada12c5f2480b Mon Sep 17 00:00:00 2001 From: Henry Yu Date: Tue, 5 Sep 2023 09:45:38 -0700 Subject: [PATCH] dmd/euler example (#12) * bindings for AdaptiveDMD and NonuniformDMD. * added ci test for DMD class. * parallel test for DMD. * fix on parallel dmd test. * initial commit * initial commit * dg_euler.py in progress... * finished dg_euler testing example * minor changes and added mesh file * fom timer corrected * minor corrections to dg_euler.py * corrections to dg_euler.py and pyAdaptiveDMD.cpp --------- Co-authored-by: Kevin Chung Co-authored-by: Henry Yu Co-authored-by: Henry Yu Co-authored-by: Henry Yu --- .github/workflows/ci.yml | 9 + CMakeLists.txt | 25 +- bindings/pylibROM/algo/pyAdaptiveDMD.cpp | 32 ++ bindings/pylibROM/algo/pyDMD.cpp | 2 +- bindings/pylibROM/algo/pyNonuniformDMD.cpp | 40 ++ bindings/pylibROM/pylibROM.cpp | 4 + examples/data/periodic-square.mesh | 85 ++++ examples/dmd/dg_euler.py | 499 +++++++++++++++++++++ examples/dmd/dg_euler_common.py | 460 +++++++++++++++++++ tests/test_DMD.py | 289 ------------ tests/test_pyDMD.py | 49 ++ 11 files changed, 1192 insertions(+), 302 deletions(-) create mode 100644 bindings/pylibROM/algo/pyAdaptiveDMD.cpp create mode 100644 bindings/pylibROM/algo/pyNonuniformDMD.cpp create mode 100644 examples/data/periodic-square.mesh create mode 100644 examples/dmd/dg_euler.py create mode 100644 examples/dmd/dg_euler_common.py delete mode 100644 tests/test_DMD.py create mode 100644 tests/test_pyDMD.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 63a9c19..d782f7c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -63,6 +63,9 @@ jobs: pytest test_pyBasisReader.py --verbose echo run pyBasisWriter unit test pytest test_pyBasisWriter.py --verbose + echo run pyDMD unit test + pytest test_pyDMD.py --verbose + mpirun -n 2 pytest test_pyDMD.py --verbose echo run pyDEIM unit test pytest test_pyDEIM.py --verbose echo run pyGNAT unit test @@ -127,6 +130,9 @@ jobs: pytest test_pyBasisReader.py --verbose echo run pyBasisWriter unit test pytest test_pyBasisWriter.py --verbose + echo run pyDMD unit test + pytest test_pyDMD.py --verbose + mpirun -n 2 pytest test_pyDMD.py --verbose echo run pyDEIM unit test pytest test_pyDEIM.py --verbose echo run pyGNAT unit test @@ -188,6 +194,9 @@ jobs: pytest test_pyBasisReader.py --verbose echo run pyBasisWriter unit test pytest test_pyBasisWriter.py --verbose + echo run pyDMD unit test + pytest test_pyDMD.py --verbose + mpirun -n 2 pytest test_pyDMD.py --verbose echo run pyDEIM unit test pytest test_pyDEIM.py --verbose echo run pyGNAT unit test diff --git a/CMakeLists.txt b/CMakeLists.txt index 4485f80..125cb3d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -128,22 +128,23 @@ link_libraries( add_subdirectory("extern/pybind11") pybind11_add_module(_pylibROM - bindings/pylibROM/pylibROM.cpp + bindings/pylibROM/pylibROM.cpp - bindings/pylibROM/linalg/pyMatrix.cpp - bindings/pylibROM/linalg/pyVector.cpp - bindings/pylibROM/linalg/pyBasisGenerator.cpp - bindings/pylibROM/linalg/pyBasisReader.cpp + bindings/pylibROM/linalg/pyMatrix.cpp + bindings/pylibROM/linalg/pyVector.cpp + bindings/pylibROM/linalg/pyBasisGenerator.cpp + bindings/pylibROM/linalg/pyBasisReader.cpp bindings/pylibROM/linalg/pyBasisWriter.cpp - bindings/pylibROM/linalg/pyOptions.cpp - bindings/pylibROM/linalg/pyNNLS.cpp - bindings/pylibROM/linalg/svd/pySVD.cpp - bindings/pylibROM/linalg/svd/pyStaticSVD.cpp - bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp + bindings/pylibROM/linalg/pyOptions.cpp + bindings/pylibROM/linalg/pyNNLS.cpp + bindings/pylibROM/linalg/svd/pySVD.cpp + bindings/pylibROM/linalg/svd/pyStaticSVD.cpp + bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp - bindings/pylibROM/algo/pyDMD.cpp + bindings/pylibROM/algo/pyDMD.cpp bindings/pylibROM/algo/pyParametricDMD.cpp - + bindings/pylibROM/algo/pyNonuniformDMD.cpp + bindings/pylibROM/algo/pyAdaptiveDMD.cpp bindings/pylibROM/hyperreduction/pyDEIM.cpp bindings/pylibROM/hyperreduction/pyGNAT.cpp bindings/pylibROM/hyperreduction/pyQDEIM.cpp diff --git a/bindings/pylibROM/algo/pyAdaptiveDMD.cpp b/bindings/pylibROM/algo/pyAdaptiveDMD.cpp new file mode 100644 index 0000000..9428ceb --- /dev/null +++ b/bindings/pylibROM/algo/pyAdaptiveDMD.cpp @@ -0,0 +1,32 @@ +#include +#include +#include "algo/AdaptiveDMD.h" +#include "linalg/Vector.h" +#include "linalg/Matrix.h" + +namespace py = pybind11; +using namespace CAROM; + +void init_AdaptiveDMD(pybind11::module_ &m){ + + py::class_(m, "AdaptiveDMD") + .def(py::init(), + py::arg("dim"), + py::arg("desired_dt") = -1.0, + py::arg("rbf") = "G", + py::arg("interp_method") = "LS", + py::arg("closest_rbf_val") = 0.9, + py::arg("alt_output_basis") = false, + py::arg("state_offset") = nullptr) + .def("train", (void (AdaptiveDMD::*)(double, const Matrix*, double)) &AdaptiveDMD::train, + py::arg("energy_fraction").noconvert(), + py::arg("W0") = nullptr, + py::arg("linearity_tol") = 0.0) + .def("train", (void (AdaptiveDMD::*)(int, const Matrix*, double)) &AdaptiveDMD::train, + py::arg("k").noconvert(), + py::arg("W0") = nullptr, + py::arg("linearity_tol") = 0.0) + .def("getTrueDt", &AdaptiveDMD::getTrueDt) + .def("getInterpolatedSnapshots", &AdaptiveDMD::getInterpolatedSnapshots); +} + diff --git a/bindings/pylibROM/algo/pyDMD.cpp b/bindings/pylibROM/algo/pyDMD.cpp index 0331833..7ba67a9 100644 --- a/bindings/pylibROM/algo/pyDMD.cpp +++ b/bindings/pylibROM/algo/pyDMD.cpp @@ -26,7 +26,7 @@ void init_DMD(pybind11::module_ &m) { return new DMD(dim, dt, alt_output_basis, vec); }), py::arg("dim"), py::arg("dt"), py::arg("alt_output_basis") = false, py::arg("vec") = nullptr) - // .def("setOffset", &PyDMD::setOffset, py::arg("offset_vector"), py::arg("order")) //problem if we want to name the wrapper as DMD. Could get rid of the using namespace directive? + .def("setOffset", &DMD::setOffset, py::arg("offset_vector"), py::arg("order")) .def("takeSample", [](DMD &self, py::array_t &u_in, double t) { self.takeSample(getVectorPointer(u_in), t); }) diff --git a/bindings/pylibROM/algo/pyNonuniformDMD.cpp b/bindings/pylibROM/algo/pyNonuniformDMD.cpp new file mode 100644 index 0000000..1ff752f --- /dev/null +++ b/bindings/pylibROM/algo/pyNonuniformDMD.cpp @@ -0,0 +1,40 @@ +// +// Created by barrow9 on 6/4/23. +// +#include +#include +#include +#include +#include "librom.h" +#include "python_utils/cpp_utils.hpp" + +namespace py = pybind11; +using namespace CAROM; + +void init_NonuniformDMD(pybind11::module_ &m) { + + py::class_(m, "NonuniformDMD") + + //constructor, default. + .def(py::init()) //constructor a + + .def(py::init([](int dim, + bool alt_output_basis, + Vector* state_offset, + Vector* derivative_offset) { + return new NonuniformDMD(dim, alt_output_basis, state_offset, derivative_offset); + }), + py::arg("dim"), + py::arg("alt_output_basis") = false, + py::arg("state_offset") = nullptr, + py::arg("derivative_offset") = nullptr) + + //TODO: needed explicitly? + .def("__del__", [](NonuniformDMD& self) { self.~NonuniformDMD(); }) // Destructor + + .def("setOffset", &NonuniformDMD::setOffset, py::arg("offset_vector"), py::arg("order")) + + .def("load", &NonuniformDMD::load, py::arg("base_file_name")) + .def("save", &NonuniformDMD::save, py::arg("base_file_name")); + +} diff --git a/bindings/pylibROM/pylibROM.cpp b/bindings/pylibROM/pylibROM.cpp index b54f3c5..a919eab 100644 --- a/bindings/pylibROM/pylibROM.cpp +++ b/bindings/pylibROM/pylibROM.cpp @@ -20,6 +20,8 @@ void init_IncrementalSVD(pybind11::module_ &m); //algo void init_DMD(pybind11::module_ &); void init_ParametricDMD(pybind11::module_ &m); +void init_NonuniformDMD(pybind11::module_ &m); +void init_AdaptiveDMD(pybind11::module_ &m); //hyperreduction void init_DEIM(pybind11::module_ &m); @@ -64,6 +66,8 @@ PYBIND11_MODULE(_pylibROM, m) { py::module algo = m.def_submodule("algo"); init_DMD(algo); init_ParametricDMD(algo); + init_AdaptiveDMD(algo); + init_NonuniformDMD(algo); py::module hyperreduction = m.def_submodule("hyperreduction"); init_DEIM(hyperreduction); diff --git a/examples/data/periodic-square.mesh b/examples/data/periodic-square.mesh new file mode 100644 index 0000000..7d0bf5a --- /dev/null +++ b/examples/data/periodic-square.mesh @@ -0,0 +1,85 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see mesh/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# + +dimension +2 + +# format: ... +elements +9 +1 3 0 1 4 3 +2 3 1 2 5 4 +3 3 2 0 3 5 +4 3 3 4 7 6 +5 3 4 5 8 7 +6 3 5 3 6 8 +7 3 6 7 1 0 +8 3 7 8 2 1 +9 3 8 6 0 2 + +boundary +0 + +vertices +9 + +nodes +FiniteElementSpace +FiniteElementCollection: L2_T1_2D_P1 +VDim: 2 +Ordering: 1 + +-1 -1 +-0.333333333 -1 +-1 -0.333333333 +-0.333333333 -0.333333333 + +-0.333333333 -1 ++0.333333333 -1 +-0.333333333 -0.333333333 ++0.333333333 -0.333333333 + ++0.333333333 -1 ++1 -1 ++0.333333333 -0.333333333 ++1 -0.333333333 + +-1 -0.333333333 +-0.333333333 -0.333333333 +-1 +0.333333333 +-0.333333333 +0.333333333 + +-0.333333333 -0.333333333 ++0.333333333 -0.333333333 +-0.333333333 +0.333333333 ++0.333333333 +0.333333333 + ++0.333333333 -0.333333333 ++1 -0.333333333 ++0.333333333 +0.333333333 ++1 +0.333333333 + +-1 +0.333333333 +-0.333333333 +0.333333333 +-1 +1 +-0.333333333 +1 + +-0.333333333 +0.333333333 ++0.333333333 +0.333333333 +-0.333333333 +1 ++0.333333333 +1 + ++0.333333333 +0.333333333 ++1 +0.333333333 ++0.333333333 +1 ++1 +1 diff --git a/examples/dmd/dg_euler.py b/examples/dmd/dg_euler.py new file mode 100644 index 0000000..2f63e1f --- /dev/null +++ b/examples/dmd/dg_euler.py @@ -0,0 +1,499 @@ +''' +****************************************************************************** +* +* Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC +* and other libROM project developers. See the top-level COPYRIGHT +* file for details. +* +* SPDX-License-Identifier: (Apache-2.0 OR MIT) +* +***************************************************************************** + +This code is adapted from PyMFEM/examples/ex18p.py, and is a mirror of +its cpp version located at libROM/examples/dmd/dg_euler.cpp + +Below is the description from libROM/examples/dmd/dg_euler.cpp: + +// Compile with: make dg_euler + +// ================================================================================= +// +// Sample runs and results for adaptive DMD: +// +// Command 1: +// mpirun -np 8 python dg_euler.py -p 1 -rs 1 -rp 1 -o 5 -s 6 -tf 0.1 -visit +// +// Output 1: +// Relative error of DMD density (dens) at t_final: 0.1 is 0.00015272589 +// Relative error of DMD x-momentum (x_mom) at t_final: 0.1 is 2.8719908e-05 +// Relative error of DMD y-momentum (y_mom) at t_final: 0.1 is 8.9435003e-05 +// Relative error of DMD energy (e) at t_final: 0.1 is 6.85403e-05 +// +// Command 2: +// mpirun -np 8 python dg_euler.py -p 2 -rs 2 -rp 1 -o 1 -s 3 -tf 0.1 -visit +// +// Output 2: +// Relative error of DMD density (dens) at t_final: 0.1 is 1.573349e-06 +// Relative error of DMD x-momentum (x_mom) at t_final: 0.1 is 4.3846865e-05 +// Relative error of DMD y-momentum (y_mom) at t_final: 0.1 is 0.0026493438 +// Relative error of DMD energy (e) at t_final: 0.1 is 1.7326842e-06 +// +// Command 3: +// mpirun -np 8 python dg_euler.py -p 2 -rs 2 -rp 1 -o 1 -s 3 -visit +// +// Output 3: +// Relative error of DMD density (dens) at t_final: 2 is 0.00022777614 +// Relative error of DMD x-momentum (x_mom) at t_final: 2 is 0.00022107792 +// Relative error of DMD y-momentum (y_mom) at t_final: 2 is 0.00030374609 +// Relative error of DMD energy (e) at t_final: 2 is 0.0002277899 +// +// ================================================================================= +// +// Sample runs and results for nonuniform DMD: +// +// Command 1: +// mpirun -np 8 python dg_euler.py -p 1 -rs 1 -rp 1 -o 5 -s 6 -tf 0.1 -nonunif -visit +// +// Output 1: +// Relative error of DMD density (dens) at t_final: 0.1 is 0.00015499558 +// Relative error of DMD x-momentum (x_mom) at t_final: 0.1 is 4.5300074e-05 +// Relative error of DMD y-momentum (y_mom) at t_final: 0.1 is 0.0034796374 +// Relative error of DMD energy (e) at t_final: 0.1 is 7.0110651e-05 +// +// Command 2: +// mpirun -np 8 python dg_euler.py -p 2 -rs 2 -rp 1 -o 1 -s 3 -tf 0.1 -nonunif -visit +// +// Output 2: +// Relative error of DMD density (dens) at t_final: 0.1 is 4.1676355e-07 +// Relative error of DMD x-momentum (x_mom) at t_final: 0.1 is 4.4263729e-05 +// Relative error of DMD y-momentum (y_mom) at t_final: 0.1 is 0.0017438412 +// Relative error of DMD energy (e) at t_final: 0.1 is 8.3869658e-07 +// +// Command 3: +// mpirun -np 8 python dg_euler.py -p 2 -rs 2 -rp 1 -o 1 -s 3 -nonunif -visit +// +// Output 3: +// Relative error of DMD density (dens) at t_final: 0.1 is 7.9616991e-07 +// Relative error of DMD x-momentum (x_mom) at t_final: 0.1 is 0.00011741735 +// Relative error of DMD y-momentum (y_mom) at t_final: 0.1 is 0.016937741 +// Relative error of DMD energy (e) at t_final: 0.1 is 2.6258626e-06 +// +// ================================================================================= +// +// Description: This example code solves the compressible Euler system of +// equations, a model nonlinear hyperbolic PDE, with a +// discontinuous Galerkin (DG) formulation. +// +// Specifically, it solves for an exact solution of the equations +// whereby a vortex is transported by a uniform flow. Since all +// boundaries are periodic here, the method's accuracy can be +// assessed by measuring the difference between the solution and +// the initial condition at a later time when the vortex returns +// to its initial location. +// +// Note that as the order of the spatial discretization increases, +// the timestep must become smaller. This example currently uses a +// simple estimate derived by Cockburn and Shu for the 1D RKDG +// method. An additional factor can be tuned by passing the --cfl +// (or -c shorter) flag. +// +// The example demonstrates user-defined bilinear and nonlinear +// form integrators for systems of equations that are defined with +// block vectors, and how these are used with an operator for +// explicit time integrators. In this case the system also +// involves an external approximate Riemann solver for the DG +// interface flux. It also demonstrates how to use GLVis for +// in-situ visualization of vector grid functions. +// +// We recommend viewing examples 9, 14 and 17 before viewing this +// example. + +''' +import mfem.par as mfem + +from dg_euler_common import FE_Evolution, InitialCondition, \ + RiemannSolver, DomainIntegrator, FaceIntegrator +from mfem.common.arg_parser import ArgParser + +from os.path import expanduser, join, dirname +import numpy as np +from numpy import sqrt, pi, cos, sin, hypot, arctan2 +from scipy.special import erfc +from ctypes import * + +# Equation constant parameters.(using globals to share them with dg_euler_common) +import dg_euler_common +from pylibROM.python_utils.StopWatch import StopWatch +from pylibROM.algo import NonuniformDMD, AdaptiveDMD + +# 1. Initialize MPI.from mpi4py import MPI + +from mpi4py import MPI +num_procs = MPI.COMM_WORLD.size +myid = MPI.COMM_WORLD.rank + +parser = ArgParser(description='dg_euler') +parser.add_argument('-m', '--mesh', + default='periodic-square.mesh', + action='store', type=str, + help='Mesh file to use.') +parser.add_argument('-p', '--problem', + action='store', default=1, type=int, + help='Problem setup to use. See options in velocity_function().') +parser.add_argument('-rs', '--refine_serial', + action='store', default=0, type=int, + help="Number of times to refine the mesh uniformly before parallel.") +parser.add_argument('-rp', '--refine_parallel', + action='store', default=1, type=int, + help="Number of times to refine the mesh uniformly after parallel.") +parser.add_argument('-o', '--order', + action='store', default=3, type=int, + help="Finite element order (polynomial degree)") +parser.add_argument('-s', '--ode_solver', + action='store', default=4, type=int, + help="ODE solver: 1 - Forward Euler,\n\t" + + " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.") +parser.add_argument('-tf', '--t_final', + action='store', default=2.0, type=float, + help="Final time; start time is 0.") +parser.add_argument("-dt", "--time_step", + action='store', default=-0.01, type=float, + help="Time step.") +parser.add_argument('-c', '--cfl_number', + action='store', default=0.3, type=float, + help="CFL number for timestep calculation.") +parser.add_argument('-vis', '--visualization', + action='store_true', + help='Enable GLVis visualization') +parser.add_argument('-visit', '--visit-datafiles', + action='store_true', default=False, + help='Save data files for VisIt (visit.llnl.gov) visualize.') +parser.add_argument('-vs', '--visualization_steps', + action='store', default=50, type=float, + help="Visualize every n-th timestep.") + +# additional args for DMD +parser.add_argument('-ef','--energy_fraction', + action='store',default=0.9999, type=float, + help='Energy fraction for DMD') +parser.add_argument('-rdim','--rdim', + action='store',default=-1,type=int, + help='Reduced dimension for DMD') +parser.add_argument('-crbf','--crbf', + action='store',default=0.9,type=float, + help='Closest RBF value') +parser.add_argument('-nonunif','--nonunif',dest='nonunif', + action='store_true',help='Use NonuniformDMD') + + +# assign args +args = parser.parse_args() +mesh = args.mesh +ser_ref_levels = args.refine_serial +par_ref_levels = args.refine_parallel +order = args.order +ode_solver_type = args.ode_solver +t_final = args.t_final +dt = args.time_step +cfl = args.cfl_number +visualization = args.visualization +visit = args.visit_datafiles +vis_steps = args.visualization_steps +ef = args.energy_fraction +rdim = args.rdim +crbf = args.crbf +nonunif = args.nonunif + +if myid == 0: + parser.print_options(args) + +device = mfem.Device('cpu') +if myid == 0: + device.Print() + +dg_euler_common.num_equation = 4 +dg_euler_common.specific_heat_ratio = 1.4 +dg_euler_common.gas_constant = 1.0 +dg_euler_common.problem = args.problem +num_equation = dg_euler_common.num_equation + + +# 3. Read the mesh from the given mesh file. This example requires a 2D +# periodic mesh, such as ../data/periodic-square.mesh. +meshfile = expanduser(join(dirname(__file__), '..', 'data', mesh)) +mesh = mfem.Mesh(meshfile, 1, 1) +dim = mesh.Dimension() + +# 4. Define the ODE solver used for time integration. Several explicit +# Runge-Kutta methods are available. +ode_solver = None +if ode_solver_type == 1: + ode_solver = mfem.ForwardEulerSolver() +elif ode_solver_type == 2: + ode_solver = mfem.RK2Solver(1.0) +elif ode_solver_type == 3: + ode_solver = mfem.RK3SSolver() +elif ode_solver_type == 4: + ode_solver = mfem.RK4Solver() +elif ode_solver_type == 6: + ode_solver = mfem.RK6Solver() +else: + print("Unknown ODE solver type: " + str(ode_solver_type)) + exit + +# 5. Refine the mesh in serial to increase the resolution. In this example +# we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is +# a command-line parameter. +for lev in range(ser_ref_levels): + mesh.UniformRefinement() + +# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine +# this mesh further in parallel to increase the resolution. Once the +# parallel mesh is defined, the serial mesh can be deleted. + +pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh) +del mesh +for lev in range(par_ref_levels): + pmesh.UniformRefinement() + +# 7. Define the discontinuous DG finite element space of the given +# polynomial order on the refined mesh. +fec = mfem.DG_FECollection(order, dim) +# Finite element space for a scalar (thermodynamic quantity) +fes = mfem.ParFiniteElementSpace(pmesh, fec) +# Finite element space for a mesh-dim vector quantity (momentum) +dfes = mfem.ParFiniteElementSpace(pmesh, fec, dim, mfem.Ordering.byNODES) +# Finite element space for all variables together (total thermodynamic state) +vfes = mfem.ParFiniteElementSpace( + pmesh, fec, num_equation, mfem.Ordering.byNODES) + +assert fes.GetOrdering() == mfem.Ordering.byNODES, "Ordering must be byNODES" +glob_size = vfes.GlobalTrueVSize() +if myid == 0: + print("Number of unknowns: " + str(glob_size)) + +# 8. Define the initial conditions, save the corresponding mesh and grid +# functions to a file. This can be opened with GLVis with the -gc option. +# The solution u has components {density, x-momentum, y-momentum, energy}. +# These are stored contiguously in the BlockVector u_block. + +offsets = [k*vfes.GetNDofs() for k in range(num_equation+1)] +offsets = mfem.intArray(offsets) +u_block = mfem.BlockVector(offsets) + +# Momentum grid function on dfes for visualization. +mom = mfem.ParGridFunction(dfes, u_block, offsets[1]) + +# Initialize the state. +u0 = InitialCondition(num_equation) +sol = mfem.ParGridFunction(vfes, u_block.GetData()) +sol.ProjectCoefficient(u0) + +smyid = '{:0>6d}'.format(myid) +pmesh.Print("vortex-mesh."+smyid, 8) +for k in range(num_equation): + uk = mfem.ParGridFunction(fes, u_block.GetBlock(k).GetData()) + sol_name = "vortex-" + str(k) + "-init."+smyid + uk.Save(sol_name, 8) + +# 9. Set up the nonlinear form corresponding to the DG discretization of the +# flux divergence, and assemble the corresponding mass matrix. +Aflux = mfem.MixedBilinearForm(dfes, fes) +Aflux.AddDomainIntegrator(DomainIntegrator(dim)) +Aflux.Assemble() + +A = mfem.ParNonlinearForm(vfes) +rsolver = RiemannSolver() +ii = FaceIntegrator(rsolver, dim) +A.AddInteriorFaceIntegrator(ii) + +# 10. Define the time-dependent evolution operator describing the ODE +# right-hand side, and perform time-integration (looping over the time +# iterations, ti, with a time-step dt). +euler = FE_Evolution(vfes, A, Aflux.SpMat()) + +if (visualization): + MPI.COMM_WORLD.Barrier() + sout = mfem.socketstream("localhost", 19916) + sout.send_text("parallel " + str(num_procs) + " " + str(myid)) + sout.precision(8) + sout.send_solution(pmesh, mom) + sout.send_text("pause") + sout.flush() + if myid == 0: + print("GLVis visualization paused.") + print(" Press space (in the GLVis window) to resume it.") + +# Determine the minimum element size. +my_hmin = 0 +if (cfl > 0): + my_hmin = min([pmesh.GetElementSize(i, 1) for i in range(pmesh.GetNE())]) +hmin = MPI.COMM_WORLD.allreduce(my_hmin, op=MPI.MIN) + +# initialize timers +fom_timer, dmd_training_timer, dmd_prediction_timer = \ + StopWatch(), StopWatch(), StopWatch() + +fom_timer.Start() +t = 0.0 +ts = [] +euler.SetTime(t) +ode_solver.Init(euler) +fom_timer.Stop() + +if (cfl > 0): + # Find a safe dt, using a temporary vector. Calling Mult() computes the + # maximum char speed at all quadrature points on all faces. + z = mfem.Vector(A.Width()) + A.Mult(sol, z) + max_char_speed = MPI.COMM_WORLD.allreduce(dg_euler_common.max_char_speed, op=MPI.MAX) + dg_euler_common.max_char_speed = max_char_speed + dt = cfl * hmin / dg_euler_common.max_char_speed / (2*order+1) + +#- DMD setup +# Initialize dmd_vars = [dmd_dens, dmd_x_mom, dmd_y_mom, dmd_e] +dmd_training_timer.Start() +if nonunif: + dmd_vars = [NonuniformDMD(u_block.GetBlock(i).Size()) + for i in range(4)] +else: + dmd_vars = [AdaptiveDMD(u_block.GetBlock(i).Size(),dt,'G','LS',crbf) \ + for i in range(4)] +for i in range(4): + dmd_vars[i].takeSample(u_block.GetBlock(i).GetDataArray(),t) +ts += [t] +dmd_training_timer.Stop() + +# Integrate in time. +done = False +ti = 0 +while not done: + + fom_timer.Start() + + dt_real = min(dt, t_final - t) + t, dt_real = ode_solver.Step(sol, t, dt_real) + + if (cfl > 0): + max_char_speed = MPI.COMM_WORLD.allreduce( + dg_euler_common.max_char_speed, op=MPI.MAX) + dg_euler_common.max_char_speed = max_char_speed + dt = cfl * hmin / dg_euler_common.max_char_speed / (2*order+1) + + ti = ti+1 + done = (t >= t_final - 1e-8*dt) + + fom_timer.Stop() + + #- DMD take sample + dmd_training_timer.Start() + for i in range(4): + dmd_vars[i].takeSample(u_block.GetBlock(i).GetDataArray(),t) + ts.append(t) + dmd_training_timer.Stop() + + + if (done or ti % vis_steps == 0): + if myid == 0: + print("time step: " + str(ti) + ", time: " + "{:g}".format(t)) + if (visualization): + sout.send_text("parallel " + str(num_procs) + " " + str(myid)) + sout.send_solution(pmesh, mom) + sout.flush() + if visit: + visit_dc.SetCycle(ti) + visit_dc.SetTime(t) + visit_dc.Save() + + +if myid == 0: + print("done") + +# 11. Save the final solution. This output can be viewed later using GLVis: +# "glvis -np 4 -m vortex-mesh -g vortex-1-final". + +for k in range(num_equation): + uk = mfem.ParGridFunction(fes, u_block.GetBlock(k).GetData()) + sol_name = "vortex-" + str(k) + "-final."+smyid + uk.Save(sol_name, 8) + +# 12. Compute the L2 solution error summed for all components. +if (t_final == 2.0): + error = sol.ComputeLpError(2., u0) + if myid == 0: + print("Solution error: " + "{:g}".format(error)) + +# 13. Calculate the DMD modes +if myid==0 and rdim != -1 and ef != -1: + print('Both rdim and ef are set. ef will be ignored') + +dmd_training_timer.Start() + +if rdim != -1: + if myid==0: + print(f'Creating DMD with rdim: {rdim}') + for dmd_var in dmd_vars: + dmd_var.train(rdim) +elif ef != -1: + if myid == 0: + print(f'Creating DMD with energy fraction: {ef}') + for dmd_var in dmd_vars: + dmd_var.train(ef) + +dmd_training_timer.Stop() + +true_solution_vars = [u_block.GetBlock(i) for i in range(4)] + +# 14. Predict the state at t_final using DMD +dmd_prediction_timer.Start() +if myid == 0: + print('Predicting density, momentum, and energy using DMD') + +result_vars = [dmd_var.predict(ts[0]) for dmd_var in dmd_vars] +initial_dmd_sols = [mfem.Vector(result_var.getData(), result_var.dim()) for result_var in result_vars] +for i in range(4): + block = u_block.GetBlock(i) + block = initial_dmd_sols[i] + #u_block.Update(initial_dmd_sols[i],offsets[i]) + +dmd_visit_dc = mfem.VisItDataCollection('DMD_DG_Euler', pmesh) +dmd_visit_dc.RegisterField('solution',mom) +if (visit): + dmd_visit_dc.SetCycle(0) + dmd_visit_dc.SetTime(0.0) + dmd_visit_dc.Save() + +if visit: + for i in range(len(ts.size)): + if (i==(len(ts)-1)) or (i%vis_steps==0): + result_vars = [var.predict(ts[i]) for var in dmd_vars] + dmd_sols = [mfem.Vector(result_var.getData(),result_var.dim()) for result_var in result_vars] + for k in range(4): + u_block.Update(dmd_sols[k],offsets[k]) + + dmd_visit_dc.SetCycle(i) + dmd_visit_dc.SetTime(ts[i]) + dmd_visit_dc.Save() +dmd_prediction_timer.Stop() +result_vars = [dmd_var.predict(t_final) for dmd_var in dmd_vars] + +# 15. Calculate the relative error between the DMD final solution and the true solution. +dmd_solution_vars = [mfem.Vector(result_var.getData(),result_var.dim()) for result_var in result_vars] +diff_vars = [mfem.Vector(result_var.dim()) for result_var in result_vars] +for i in range(4): + mfem.subtract_vector(dmd_solution_vars[i],true_solution_vars[i],\ + diff_vars[i]) +tot_diff_norm_vars = [sqrt(mfem.InnerProduct(MPI.COMM_WORLD,diff_var,diff_var)) for diff_var in diff_vars] +tot_true_solution_norm_vars = [sqrt(mfem.InnerProduct(MPI.COMM_WORLD,true_solution_var,true_solution_var))\ + for true_solution_var in true_solution_vars] + +if myid==0: + var_names = ['dens', 'x_mom', 'y_mom', 'e'] + for i in range(len(var_names)): + rel_error = tot_diff_norm_vars[i]/tot_true_solution_norm_vars[i] + print(f'Relative error of DMD {var_names[i]} at t_final: {t_final} is {rel_error:.10f}') + + print(f'Elapsed time for solving FOM: {fom_timer.duration:.6e}') + print(f'Elapsed time for training DMD: {dmd_training_timer.duration:.6e}') + print(f'Elapsed time for predicting DMD: {dmd_prediction_timer.duration:.6e}') + diff --git a/examples/dmd/dg_euler_common.py b/examples/dmd/dg_euler_common.py new file mode 100644 index 0000000..770d17d --- /dev/null +++ b/examples/dmd/dg_euler_common.py @@ -0,0 +1,460 @@ +''' +****************************************************************************** +* +* Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC +* and other libROM project developers. See the top-level COPYRIGHT +* file for details. +* +* SPDX-License-Identifier: (Apache-2.0 OR MIT) +* +***************************************************************************** + +This code is adapted from PyMFEM/examples/ex18_common.py, and is a mirror of +its cpp version located at libROM/examples/dmd/dg_euler.hpp. + +Original description from ex18_common.py is listed below: + +>>> + MFEM Example 18 - Serial/Parallel Shared Code + + This is a python translation of ex18.hpp + + note: following variabls are set from ex18 or ex18p + problem + num_equation + max_char_speed + specific_heat_ratio; + gas_constant; +<<< +''' + +import numpy as np +from mfem import mfem_mode + +if mfem_mode is None or mfem_mode == 'serial': + import mfem.ser as mfem +else: + import mfem.par as mfem + +num_equation = 0 +specific_heat_ratio = 0 +gas_constant = 0 +problem = 0 +max_char_speed = 0 + + +class FE_Evolution(mfem.TimeDependentOperator): + def __init__(self, vfes, A, A_flux): + self.dim = vfes.GetFE(0).GetDim() + self.vfes = vfes + self.A = A + self.Aflux = A_flux + self.Me_inv = mfem.DenseTensor(vfes.GetFE(0).GetDof(), + vfes.GetFE(0).GetDof(), + vfes.GetNE()) + + self.state = mfem.Vector(num_equation) + self.f = mfem.DenseMatrix(num_equation, self.dim) + self.flux = mfem.DenseTensor(vfes.GetNDofs(), self.dim, num_equation) + self.z = mfem.Vector(A.Height()) + + dof = vfes.GetFE(0).GetDof() + Me = mfem.DenseMatrix(dof) + inv = mfem.DenseMatrixInverse(Me) + mi = mfem.MassIntegrator() + for i in range(vfes.GetNE()): + mi.AssembleElementMatrix(vfes.GetFE( + i), vfes.GetElementTransformation(i), Me) + inv.Factor() + inv.GetInverseMatrix(self.Me_inv(i)) + super(FE_Evolution, self).__init__(A.Height()) + + def GetFlux(self, x, flux): + state = self.state + dof = self.flux.SizeI() + dim = self.flux.SizeJ() + + flux_data = [] + for i in range(dof): + for k in range(num_equation): + self.state[k] = x[i, k] + ComputeFlux(state, dim, self.f) + + flux_data.append(self.f.GetDataArray().transpose().copy()) + # flux[i].Print() + # print(self.f.GetDataArray()) + # for d in range(dim): + # for k in range(num_equation): + # flux[i, d, k] = self.f[k, d] + + mcs = ComputeMaxCharSpeed(state, dim) + if (mcs > globals()['max_char_speed']): + globals()['max_char_speed'] = mcs + + flux.Assign(np.stack(flux_data)) + #print("max char speed", globals()['max_char_speed']) + + def Mult(self, x, y): + globals()['max_char_speed'] = 0. + num_equation = globals()['num_equation'] + # 1. Create the vector z with the face terms -. + self.A.Mult(x, self.z) + + # 2. Add the element terms. + # i. computing the flux approximately as a grid function by interpolating + # at the solution nodes. + # ii. multiplying this grid function by a (constant) mixed bilinear form for + # each of the num_equation, computing (F(u), grad(w)) for each equation. + + xmat = mfem.DenseMatrix( + x.GetData(), self.vfes.GetNDofs(), num_equation) + self.GetFlux(xmat, self.flux) + + for k in range(num_equation): + fk = mfem.Vector(self.flux[k].GetData(), + self.dim * self.vfes.GetNDofs()) + o = k * self.vfes.GetNDofs() + zk = self.z[o: o+self.vfes.GetNDofs()] + self.Aflux.AddMult(fk, zk) + + # 3. Multiply element-wise by the inverse mass matrices. + zval = mfem.Vector() + vdofs = mfem.intArray() + dof = self.vfes.GetFE(0).GetDof() + zmat = mfem.DenseMatrix() + ymat = mfem.DenseMatrix(dof, num_equation) + + for i in range(self.vfes.GetNE()): + # Return the vdofs ordered byNODES + vdofs = mfem.intArray(self.vfes.GetElementVDofs(i)) + self.z.GetSubVector(vdofs, zval) + zmat.UseExternalData(zval.GetData(), dof, num_equation) + mfem.Mult(self.Me_inv[i], zmat, ymat) + y.SetSubVector(vdofs, ymat.GetData()) + + +class DomainIntegrator(mfem.BilinearFormIntegrator): + def __init__(self, dim): + num_equation = globals()['num_equation'] + self.flux = mfem.DenseMatrix(num_equation, dim) + self.shape = mfem.Vector() + self.dshapedr = mfem.DenseMatrix() + self.dshapedx = mfem.DenseMatrix() + super(DomainIntegrator, self).__init__() + + def AssembleElementMatrix2(self, trial_fe, test_fe, Tr, elmat): + # Assemble the form (vec(v), grad(w)) + + # Trial space = vector L2 space (mesh dim) + # Test space = scalar L2 space + + dof_trial = trial_fe.GetDof() + dof_test = test_fe.GetDof() + dim = trial_fe.GetDim() + + self.shape.SetSize(dof_trial) + self.dshapedr.SetSize(dof_test, dim) + self.dshapedx.SetSize(dof_test, dim) + + elmat.SetSize(dof_test, dof_trial * dim) + elmat.Assign(0.0) + + maxorder = max(trial_fe.GetOrder(), test_fe.GetOrder()) + intorder = 2 * maxorder + ir = mfem.IntRules.Get(trial_fe.GetGeomType(), intorder) + + for i in range(ir.GetNPoints()): + ip = ir.IntPoint(i) + + # Calculate the shape functions + trial_fe.CalcShape(ip, self.shape) + self.shape *= ip.weight + + # Compute the physical gradients of the test functions + Tr.SetIntPoint(ip) + test_fe.CalcDShape(ip, self.dshapedr) + mfem.Mult(self.dshapedr, Tr.AdjugateJacobian(), self.dshapedx) + + for d in range(dim): + for j in range(dof_test): + for k in range(dof_trial): + elmat[j, k + d * dof_trial] += self.shape[k] * \ + self.dshapedx[j, d] + + +class FaceIntegrator(mfem.NonlinearFormIntegrator): + def __init__(self, rsolver, dim): + self.rsolver = rsolver + self.shape1 = mfem.Vector() + self.shape2 = mfem.Vector() + self.funval1 = mfem.Vector(num_equation) + self.funval2 = mfem.Vector(num_equation) + self.nor = mfem.Vector(dim) + self.fluxN = mfem.Vector(num_equation) + self.eip1 = mfem.IntegrationPoint() + self.eip2 = mfem.IntegrationPoint() + super(FaceIntegrator, self).__init__() + + self.fluxNA = np.atleast_2d(self.fluxN.GetDataArray()) + + def AssembleFaceVector(self, el1, el2, Tr, elfun, elvect): + num_equation = globals()['num_equation'] + # Compute the term on the interior faces. + dof1 = el1.GetDof() + dof2 = el2.GetDof() + + self.shape1.SetSize(dof1) + self.shape2.SetSize(dof2) + + elvect.SetSize((dof1 + dof2) * num_equation) + elvect.Assign(0.0) + + elfun1_mat = mfem.DenseMatrix(elfun.GetData(), dof1, num_equation) + elfun2_mat = mfem.DenseMatrix( + elfun[dof1*num_equation:].GetData(), dof2, num_equation) + + elvect1_mat = mfem.DenseMatrix(elvect.GetData(), dof1, num_equation) + elvect2_mat = mfem.DenseMatrix( + elvect[dof1*num_equation:].GetData(), dof2, num_equation) + + # Integration order calculation from DGTraceIntegrator + if (Tr.Elem2No >= 0): + intorder = (min(Tr.Elem1.OrderW(), Tr.Elem2.OrderW()) + + 2*max(el1.GetOrder(), el2.GetOrder())) + else: + intorder = Tr.Elem1.OrderW() + 2*el1.GetOrder() + + if (el1.Space() == mfem.FunctionSpace().Pk): + intorder += 1 + + ir = mfem.IntRules.Get(Tr.GetGeometryType(), int(intorder)) + + mat1A = elvect1_mat.GetDataArray() + mat2A = elvect2_mat.GetDataArray() + shape1A = np.atleast_2d(self.shape1.GetDataArray()) + shape2A = np.atleast_2d(self.shape2.GetDataArray()) + + for i in range(ir.GetNPoints()): + ip = ir.IntPoint(i) + Tr.Loc1.Transform(ip, self.eip1) + Tr.Loc2.Transform(ip, self.eip2) + + # Calculate basis functions on both elements at the face + el1.CalcShape(self.eip1, self.shape1) + el2.CalcShape(self.eip2, self.shape2) + + # Interpolate elfun at the point + elfun1_mat.MultTranspose(self.shape1, self.funval1) + elfun2_mat.MultTranspose(self.shape2, self.funval2) + Tr.Face.SetIntPoint(ip) + + # Get the normal vector and the flux on the face + + mfem.CalcOrtho(Tr.Face.Jacobian(), self.nor) + + mcs = self.rsolver.Eval( + self.funval1, self.funval2, self.nor, self.fluxN) + + # Update max char speed + if mcs > globals()['max_char_speed']: + globals()['max_char_speed'] = mcs + + self.fluxN *= ip.weight + + # + mat1A -= shape1A.transpose().dot(self.fluxNA) + mat2A += shape2A.transpose().dot(self.fluxNA) + ''' + for k in range(num_equation): + for s in range(dof1): + elvect1_mat[s, k] -= self.fluxN[k] * self.shape1[s] + for s in range(dof2): + elvect2_mat[s, k] += self.fluxN[k] * self.shape2[s] + ''' + + +class RiemannSolver(object): + def __init__(self): + num_equation = globals()['num_equation'] + self.flux1 = mfem.Vector(num_equation) + self.flux2 = mfem.Vector(num_equation) + + def Eval(self, state1, state2, nor, flux): + + # NOTE: nor in general is not a unit normal + dim = nor.Size() + + assert StateIsPhysical(state1, dim), "" + assert StateIsPhysical(state2, dim), "" + + maxE1 = ComputeMaxCharSpeed(state1, dim) + maxE2 = ComputeMaxCharSpeed(state2, dim) + maxE = max(maxE1, maxE2) + + ComputeFluxDotN(state1, nor, self.flux1) + ComputeFluxDotN(state2, nor, self.flux2) + + #normag = np.sqrt(np.sum(nor.GetDataArray()**2)) + normag = nor.Norml2() + + ''' + for i in range(num_equation): + flux[i] = (0.5 * (self.flux1[i] + self.flux2[i]) + - 0.5 * maxE * (state2[i] - state1[i]) * normag) + ''' + f = (0.5 * (self.flux1.GetDataArray() + self.flux2.GetDataArray()) + - 0.5 * maxE * (state2.GetDataArray() - state1.GetDataArray()) * normag) + flux.Assign(f) + + return maxE + + +def StateIsPhysical(state, dim): + specific_heat_ratio = globals()["specific_heat_ratio"] + + den = state[0] + #den_vel = state.GetDataArray()[1:1+dim] + den_energy = state[1 + dim] + + if (den < 0): + print("Negative density: " + str(state.GetDataArray())) + return False + if (den_energy <= 0): + print("Negative energy: " + str(state.GetDataArray())) + return False + + #den_vel2 = np.sum(den_vel**2)/den + den_vel2 = (state[1:1+dim].Norml2())**2/den + pres = (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2) + if (pres <= 0): + print("Negative pressure: " + str(state.GetDataArray())) + return False + return True + + +class InitialCondition(mfem.VectorPyCoefficient): + def __init__(self, dim): + mfem.VectorPyCoefficient.__init__(self, dim) + + def EvalValue(self, x): + dim = x.shape[0] + assert dim == 2, "" + problem = globals()['problem'] + if (problem == 1): + # "Fast vortex" + radius = 0.2 + Minf = 0.5 + beta = 1. / 5. + elif (problem == 2): + # "Slow vortex" + radius = 0.2 + Minf = 0.05 + beta = 1. / 50. + else: + assert False, "Cannot recognize problem. Options are: 1 - fast vortex, 2 - slow vortex" + + xc = 0.0 + yc = 0.0 + # Nice units + vel_inf = 1. + den_inf = 1. + + specific_heat_ratio = globals()["specific_heat_ratio"] + gas_constant = globals()["gas_constant"] + + pres_inf = (den_inf / specific_heat_ratio) * \ + (vel_inf / Minf) * (vel_inf / Minf) + temp_inf = pres_inf / (den_inf * gas_constant) + + r2rad = 0.0 + r2rad += (x[0] - xc) * (x[0] - xc) + r2rad += (x[1] - yc) * (x[1] - yc) + r2rad /= (radius * radius) + + shrinv1 = 1.0 / (specific_heat_ratio - 1.) + + velX = vel_inf * \ + (1 - beta * (x[1] - yc) / radius * np.exp(-0.5 * r2rad)) + velY = vel_inf * beta * (x[0] - xc) / radius * np.exp(-0.5 * r2rad) + vel2 = velX * velX + velY * velY + + specific_heat = gas_constant * specific_heat_ratio * shrinv1 + + temp = temp_inf - (0.5 * (vel_inf * beta) * + (vel_inf * beta) / specific_heat * np.exp(-r2rad)) + + den = den_inf * (temp/temp_inf)**shrinv1 + pres = den * gas_constant * temp + energy = shrinv1 * pres / den + 0.5 * vel2 + + y = np.array([den, den * velX, den * velY, den * energy]) + return y + + +def ComputePressure(state, dim): + den = state[0] + #den_vel = state.GetDataArray()[1:1+dim] + den_energy = state[1 + dim] + + specific_heat_ratio = globals()["specific_heat_ratio"] + #den_vel2 = np.sum(den_vel**2)/den + den_vel2 = (state[1:1+dim].Norml2())**2/den + pres = (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2) + + return pres + + +def ComputeFlux(state, dim, flux): + den = state[0] + den_vel = state.GetDataArray()[1:1+dim] + den_energy = state[1 + dim] + + assert StateIsPhysical(state, dim), "" + + pres = ComputePressure(state, dim) + + den_vel2 = np.atleast_2d(den_vel) + fluxA = flux.GetDataArray() + fluxA[0, :] = den_vel + fluxA[1:1+dim, :] = den_vel2.transpose().dot(den_vel2) / den + for d in range(dim): + fluxA[1+d, d] += pres + + H = (den_energy + pres) / den + flux.GetDataArray()[1+dim, :] = den_vel * H + + +def ComputeFluxDotN(state, nor, fluxN): + # NOTE: nor in general is not a unit normal + dim = nor.Size() + nor = nor.GetDataArray() + fluxN = fluxN.GetDataArray() + + den = state[0] + den_vel = state.GetDataArray()[1:1+dim] + den_energy = state[1 + dim] + + assert StateIsPhysical(state, dim), "" + + pres = ComputePressure(state, dim) + + den_velN = den_vel.dot(nor) + + fluxN[0] = den_velN + fluxN[1:1+dim] = den_velN * den_vel / den + pres * nor + + H = (den_energy + pres) / den + fluxN[1+dim] = den_velN * H + + +def ComputeMaxCharSpeed(state, dim): + specific_heat_ratio = globals()["specific_heat_ratio"] + + den = state[0] + den_vel2 = (state[1:1+dim].Norml2())**2/den + pres = ComputePressure(state, dim) + + sound = np.sqrt(specific_heat_ratio * pres / den) + vel = np.sqrt(den_vel2 / den) + + return vel + sound diff --git a/tests/test_DMD.py b/tests/test_DMD.py deleted file mode 100644 index ca5bcd3..0000000 --- a/tests/test_DMD.py +++ /dev/null @@ -1,289 +0,0 @@ -''' - adapted from pyMFEM example 16 - - How to run: - mpirun -np 4 python - - Example of arguments: - ex16p.py -m inline-tri.mesh - ex16p.py -m disc-nurbs.mesh -tf 2 - ex16p.py -s 1 -a 0.0 -k 1.0 - ex16p.py -s 2 -a 1.0 -k 0.0 - ex16p.py -s 3 -a 0.5 -k 0.5 -o 4 - ex16p.py -s 14 -dt 1.0e-4 -tf 4.0e-2 -vs 40 - ex16p.py -m fichera-q2.mesh - ex16p.py -m escher.mesh - ex16p.py -m beam-tet.mesh -tf 10 -dt 0.1 - ex16p.py -m amr-quad.mesh -o 4 -rs 0 -rp 0 - ex16p.py -m amr-hex.mesh -o 2 -rs 0 -rp 0 -''' -import sys -from mfem.common.arg_parser import ArgParser -from os.path import expanduser, join, dirname -import numpy as np -from numpy import sqrt, pi, cos, sin, hypot, arctan2 -from scipy.special import erfc - -from mfem.par import intArray, add_vector -import mfem.par as mfem -from mpi4py import MPI - -num_procs = MPI.COMM_WORLD.size -myid = MPI.COMM_WORLD.rank - - -class ConductionOperator(mfem.PyTimeDependentOperator): - def __init__(self, fespace, alpha, kappa, u): - mfem.PyTimeDependentOperator.__init__( - self, fespace.GetTrueVSize(), 0.0) - rel_tol = 1e-8 - self.alpha = alpha - self.kappa = kappa - self.T = None - self.K = None - self.M = None - self.fespace = fespace - - self.ess_tdof_list = intArray() - self.Mmat = mfem.HypreParMatrix() - self.Kmat = mfem.HypreParMatrix() - self.M_solver = mfem.CGSolver(fespace.GetComm()) - self.M_prec = mfem.HypreSmoother() - self.T_solver = mfem.CGSolver(fespace.GetComm()) - self.T_prec = mfem.HypreSmoother() - self.z = mfem.Vector(self.Height()) - - self.M = mfem.ParBilinearForm(fespace) - self.M.AddDomainIntegrator(mfem.MassIntegrator()) - self.M.Assemble() - self.M.FormSystemMatrix(self.ess_tdof_list, self.Mmat) - - self.M_solver.iterative_mode = False - self.M_solver.SetRelTol(rel_tol) - self.M_solver.SetAbsTol(0.0) - self.M_solver.SetMaxIter(100) - self.M_solver.SetPrintLevel(0) - self.M_prec.SetType(mfem.HypreSmoother.Jacobi) - self.M_solver.SetPreconditioner(self.M_prec) - self.M_solver.SetOperator(self.Mmat) - - self.T_solver.iterative_mode = False - self.T_solver.SetRelTol(rel_tol) - self.T_solver.SetAbsTol(0.0) - self.T_solver.SetMaxIter(100) - self.T_solver.SetPrintLevel(0) - self.T_solver.SetPreconditioner(self.T_prec) - - self.SetParameters(u) - - def Mult(self, u, u_dt): - # Compute: - # du_dt = M^{-1}*-K(u) for du_dt - self.Kmat.Mult(u, self.z) - self.z.Neg() # z = -z - self.M_solver.Mult(self.z, du_dt) - - def ImplicitSolve(self, dt, u, du_dt): - # Solve the equation: - # du_dt = M^{-1}*[-K(u + dt*du_dt)] - # for du_dt - if self.T is None: - self.T = mfem.Add(1.0, self.Mmat, dt, self.Kmat) - current_dt = dt - self.T_solver.SetOperator(self.T) - self.Kmat.Mult(u, self.z) - self.z.Neg() - self.T_solver.Mult(self.z, du_dt) - - def SetParameters(self, u): - u_alpha_gf = mfem.ParGridFunction(self.fespace) - u_alpha_gf.SetFromTrueDofs(u) - for i in range(u_alpha_gf.Size()): - u_alpha_gf[i] = self.kappa + self.alpha * u_alpha_gf[i] - - self.K = mfem.ParBilinearForm(self.fespace) - u_coeff = mfem.GridFunctionCoefficient(u_alpha_gf) - self.K.AddDomainIntegrator(mfem.DiffusionIntegrator(u_coeff)) - self.K.Assemble(0) - self.K.FormSystemMatrix(self.ess_tdof_list, self.Kmat) - self.T = None - - -class InitialTemperature(mfem.PyCoefficient): - def EvalValue(self, x): - xx = np.array(x) - norm2 = np.sqrt(float(np.sum(xx**2))) - if norm2 < 0.5: - return 2.0 - return 1.0 - - -parser = ArgParser(description='Ex16p') -parser.add_argument('-m', '--mesh', - default='star.mesh', - action='store', type=str, - help='Mesh file to use.') -parser.add_argument('-rs', '--refine-serial', - action='store', default=2, type=int, - help="Number of times to refine the mesh uniformly in serial") -parser.add_argument('-rp', '--refine-parallel', - action='store', default=1, type=int, - help="Number of times to refine the mesh uniformly in parallel") - -parser.add_argument('-o', '--order', - action='store', default=2, type=int, - help="Finite element order (polynomial degree)") -parser.add_argument('-s', '--ode-solver', - action='store', default=3, type=int, - help='\n'.join(["ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3", - "\t\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4."])) -parser.add_argument('-t', '--t-final', - action='store', default=0.5, type=float, - help="Final time; start time is 0.") -parser.add_argument("-dt", "--time-step", - action='store', default=0.01, type=float, - help="Time step.") -parser.add_argument('-a', '--alpha', - action='store', default=0.01, type=float, - help='Alpha coefficient') -parser.add_argument('-k', '--kappa', - action='store', default=0.5, type=float, - help='Kappa coefficient') -parser.add_argument('-vis', '--visualization', - action='store_true', - help='Enable GLVis visualization') -parser.add_argument('-vs', '--visualization-steps', - action='store', default=5, type=int, - help="Visualize every n-th timestep.") - -args = parser.parse_args() -ser_ref_levels = args.refine_serial -par_ref_levels = args.refine_parallel -order = args.order -dt = args.time_step -t_final = args.t_final -alpha = args.alpha -kappa = args.kappa -visualization = args.visualization -vis_steps = args.visualization_steps -ode_solver_type = args.ode_solver -if myid == 0: - parser.print_options(args) - -device = mfem.Device('cpu') -if myid == 0: - device.Print() - -# 3. Read the serial mesh from the given mesh file on all processors. We can -# handle triangular, quadrilateral, tetrahedral and hexahedral meshes -# with the same code. -meshfile = expanduser(join(dirname(__file__), '..', 'data', args.mesh)) -mesh = mfem.Mesh(meshfile, 1, 1) -dim = mesh.Dimension() - -# 4. Define the ODE solver used for time integration. Several implicit -# singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as -# explicit Runge-Kutta methods are available. -if ode_solver_type == 1: - ode_solver = BackwardEulerSolver() -elif ode_solver_type == 2: - ode_solver = mfem.SDIRK23Solver(2) -elif ode_solver_type == 3: - ode_solver = mfem.SDIRK33Solver() -elif ode_solver_type == 11: - ode_solver = ForwardEulerSolver() -elif ode_solver_type == 12: - ode_solver = mfem.RK2Solver(0.5) -elif ode_solver_type == 13: - ode_solver = mfem.RK3SSPSolver() -elif ode_solver_type == 14: - ode_solver = mfem.RK4Solver() -elif ode_solver_type == 22: - ode_solver = mfem.ImplicitMidpointSolver() -elif ode_solver_type == 23: - ode_solver = mfem.SDIRK23Solver() -elif ode_solver_type == 24: - ode_solver = mfem.SDIRK34Solver() -else: - print("Unknown ODE solver type: " + str(ode_solver_type)) - exit -# 5. Refine the mesh in serial to increase the resolution. In this example -# we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is -# a command-line parameter. -for lev in range(ser_ref_levels): - mesh.UniformRefinement() - -# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine -# this mesh further in parallel to increase the resolution. Once the -# parallel mesh is defined, the serial mesh can be deleted. -pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh) -del mesh -for x in range(par_ref_levels): - pmesh.UniformRefinement() - -# 7. Define the vector finite element space representing the current and the -# initial temperature, u_ref. -fe_coll = mfem.H1_FECollection(order, dim) -fespace = mfem.ParFiniteElementSpace(pmesh, fe_coll) - -fe_size = fespace.GlobalTrueVSize() -if myid == 0: - print("Number of temperature unknowns: " + str(fe_size)) -u_gf = mfem.ParGridFunction(fespace) - -# 8. Set the initial conditions for u. All boundaries are considered -# natural. -u_0 = InitialTemperature() -u_gf.ProjectCoefficient(u_0) -u = mfem.Vector() -u_gf.GetTrueDofs(u) - -# 9. Initialize the conduction operator and the visualization. -oper = ConductionOperator(fespace, alpha, kappa, u) -u_gf.SetFromTrueDofs(u) - -smyid = '{:0>6d}'.format(myid) -mesh_name = "ex16-mesh."+smyid -sol_name = "ex16-init."+smyid -pmesh.Print(mesh_name, 8) -u_gf.Save(sol_name, 8) - -if visualization: - sout = mfem.socketstream("localhost", 19916) - sout.send_text("parallel " + str(num_procs) + " " + str(myid)) - sout.precision(8) - sout.send_solution(pmesh, u_gf) - sout.send_text("pause") - sout.flush() - if myid == 0: - print("GLVis visualization paused.") - print(" Press space (in the GLVis window) to resume it.") - -# 8. Perform time-integration (looping over the time iterations, ti, with a -# time-step dt). -ode_solver.Init(oper) -t = 0.0 -ti = 1 -last_step = False -while not last_step: - - if (t + dt >= t_final - dt/2): - last_step = True - - t, dt = ode_solver.Step(u, t, dt) - - if (last_step or (ti % vis_steps) == 0): - if myid == 0: - print("step " + str(ti) + ", t = " + "{:g}".format(t)) - u_gf.SetFromTrueDofs(u) - if (visualization): - sout.send_text("parallel " + str(num_procs) + " " + str(myid)) - sout.send_solution(pmesh, u_gf) - sout.flush() - ti = ti + 1 - oper.SetParameters(u) - -# 11. Save the final solution in parallel. This output can be viewed later -# using GLVis: "glvis -np -m ex16-mesh -g ex16-final". -sol_name = "ex16-final."+smyid -u_gf.Save(sol_name, 8) \ No newline at end of file diff --git a/tests/test_pyDMD.py b/tests/test_pyDMD.py new file mode 100644 index 0000000..7a93168 --- /dev/null +++ b/tests/test_pyDMD.py @@ -0,0 +1,49 @@ +import pytest +import numpy as np +import sys +try: + # import pip-installed package + import pylibROM + import pylibROM.linalg as linalg + import pylibROM.algo as algo + import pylibROM.utils as utils +except ModuleNotFoundError: + # If pip-installed package is not found, import cmake-built package + sys.path.append("../build") + import _pylibROM as pylibROM + import _pylibROM.linalg as linalg + import _pylibROM.algo as algo + import _pylibROM.utils as utils + +def test_DMD(): + from mpi4py import MPI + d_rank = MPI.COMM_WORLD.Get_rank() + d_num_procs = MPI.COMM_WORLD.Get_size() + + num_total_rows = 5 + d_num_rows = utils.split_dimension(num_total_rows, MPI.COMM_WORLD) + num_total_rows_check, row_offset = utils.get_global_offsets(d_num_rows, MPI.COMM_WORLD) + assert(num_total_rows == num_total_rows_check) + + samples = [[0.5377, 1.8339, -2.2588, 0.8622, 0.3188], + [-1.3077, -0.4336, 0.3426, 3.5784, 2.7694], + [-1.3499, 3.0349, 0.7254, -0.0631, 0.7147]] + prediction_baseline = [-0.4344, -0.0974, 0.0542, 1.2544, 0.9610] + + dmd = algo.DMD(d_num_rows, 1.0) + for k, sample in enumerate(samples): + dmd.takeSample(sample[row_offset[d_rank]:row_offset[d_rank]+d_num_rows], k * 1.0) + + dmd.train(2) + result = dmd.predict(3.0) + # print("rank: %d, " % d_rank, result.getData()) + assert np.allclose(result.getData(), prediction_baseline[row_offset[d_rank]:row_offset[d_rank]+d_num_rows], atol=1e-3) + + dmd.save("test_DMD") + dmd_load = algo.DMD("test_DMD") + result_load = dmd_load.predict(3.0) + + assert np.allclose(result_load.getData(), prediction_baseline[row_offset[d_rank]:row_offset[d_rank]+d_num_rows], atol=1e-3) + +if __name__ == '__main__': + pytest.main()