diff --git a/README.md b/README.md index db98075..6bdb9a2 100644 --- a/README.md +++ b/README.md @@ -46,6 +46,29 @@ python3 setup.py install --with-parallel --with-gslib --user ``` Make sure [`swig`](https://pypi.org/project/swig) is installed first. Also, the binary file must be located in `PATH` environment variable. +### pylibROM-Jupyter Docker Container + +This Docker container provides an environment with Jupyter Notebook for the pylibROM library. Follow the steps below to build the Docker image and run a Jupyter Notebook server. + +#### Build the Docker Image + +Navigate to the directory containing the Dockerfile: + ``` +cd /path/to/folder/pylibROM/docker/jupyter + ``` + +Now, run the following command to build the Docker image: + + ``` +docker build -t pylibrom-jupyter:latest . + ``` + +Once the image is built, you can run a container and start a Jupyter Notebook server. Replace /path/to/host/folder with the absolute path to the local directory you want to mount inside the container for Jupyter notebooks: + + ``` +docker run -p 8888:8888 -v /path/to/host/folder:/notebooks -w /notebooks pylibrom-jupyter:latest + ``` + ## License pylibROM is distributed under the terms of the MIT license. All new contributions must be made under the MIT. See diff --git a/docker/jupyter/Dockerfile b/docker/jupyter/Dockerfile new file mode 100644 index 0000000..3aaf0d3 --- /dev/null +++ b/docker/jupyter/Dockerfile @@ -0,0 +1,31 @@ +# start from pylibrom docker container +FROM --platform=linux/amd64 ghcr.io/llnl/pylibrom/pylibrom_env:latest +ENV ENVDIR=env +ENV LIB_DIR=/$ENVDIR/dependencies +WORKDIR $LIB_DIR + +# WORKDIR /env/dependencies +RUN sudo git clone --recursive https://github.com/LLNL/pylibROM.git +WORKDIR pylibROM +RUN sudo -E pip install ./ --global-option="--librom_dir=/env/dependencies/libROM" + +# Install Jupyter Notebook +RUN sudo apt-get install -yq python3-pip +RUN sudo pip3 install jupyter + +# Create a directory for Jupyter notebooks +RUN mkdir /home/$USERNAME/notebooks +WORKDIR /home/$USERNAME/notebooks + +# Configure Jupyter Notebook +RUN jupyter notebook --generate-config +RUN echo "c.NotebookApp.ip = '*'" >> /home/$USERNAME/.jupyter/jupyter_notebook_config.py + +# Expose the Jupyter Notebook port +EXPOSE 8888 + +# Run Jupyter Notebook +CMD ["jupyter", "notebook", "--ip=0.0.0.0", "--port=8888", "--no-browser", "--allow-root"] + +# create and switch to a user +WORKDIR /home/$USERNAME diff --git a/examples/dmd/dg_euler.ipynb b/examples/dmd/dg_euler.ipynb new file mode 100644 index 0000000..7252acc --- /dev/null +++ b/examples/dmd/dg_euler.ipynb @@ -0,0 +1,564 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "7ddc2f06-881c-4814-9752-063116769403", + "metadata": {}, + "outputs": [], + "source": [ + "import mfem.par as mfem\n", + "\n", + "from dg_euler_common import FE_Evolution, InitialCondition, \\\n", + " RiemannSolver, DomainIntegrator, FaceIntegrator\n", + "from mfem.common.arg_parser import ArgParser\n", + "\n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sqrt, pi, cos, sin, hypot, arctan2\n", + "from scipy.special import erfc\n", + "from ctypes import *\n", + "\n", + "# Equation constant parameters.(using globals to share them with dg_euler_common)\n", + "import dg_euler_common\n", + "from pylibROM.python_utils.StopWatch import StopWatch\n", + "from pylibROM.algo import NonuniformDMD, AdaptiveDMD\n", + "from mpi4py import MPI\n", + "import os\n", + "num_procs = MPI.COMM_WORLD.size\n", + "myid = MPI.COMM_WORLD.rank" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb74c4b1-9055-4b1e-af71-6317b423c0f2", + "metadata": {}, + "outputs": [], + "source": [ + "parser = ArgParser(description='dg_euler')\n", + "parser.add_argument('-m', '--mesh',\n", + " default='periodic-square.mesh',\n", + " action='store', type=str,\n", + " help='Mesh file to use.')\n", + "parser.add_argument('-p', '--problem',\n", + " action='store', default=1, type=int,\n", + " help='Problem setup to use. See options in velocity_function().')\n", + "parser.add_argument('-rs', '--refine_serial',\n", + " action='store', default=0, type=int,\n", + " help=\"Number of times to refine the mesh uniformly before parallel.\")\n", + "parser.add_argument('-rp', '--refine_parallel',\n", + " action='store', default=1, type=int,\n", + " help=\"Number of times to refine the mesh uniformly after parallel.\")\n", + "parser.add_argument('-o', '--order',\n", + " action='store', default=3, type=int,\n", + " help=\"Finite element order (polynomial degree)\")\n", + "parser.add_argument('-s', '--ode_solver',\n", + " action='store', default=4, type=int,\n", + " help=\"ODE solver: 1 - Forward Euler,\\n\\t\" +\n", + " \" 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.\")\n", + "parser.add_argument('-tf', '--t_final',\n", + " action='store', default=2.0, type=float,\n", + " help=\"Final time; start time is 0.\")\n", + "parser.add_argument(\"-dt\", \"--time_step\",\n", + " action='store', default=-0.01, type=float,\n", + " help=\"Time step.\")\n", + "parser.add_argument('-c', '--cfl_number',\n", + " action='store', default=0.3, type=float,\n", + " help=\"CFL number for timestep calculation.\")\n", + "parser.add_argument('-vis', '--visualization',\n", + " action='store_true',\n", + " help='Enable GLVis visualization')\n", + "parser.add_argument('-visit', '--visit-datafiles',\n", + " action='store_true', default=False,\n", + " help='Save data files for VisIt (visit.llnl.gov) visualize.')\n", + "parser.add_argument('-vs', '--visualization_steps',\n", + " action='store', default=50, type=float,\n", + " help=\"Visualize every n-th timestep.\")\n", + "\n", + "# additional args for DMD\n", + "parser.add_argument('-ef','--energy_fraction',\n", + " action='store',default=0.9999, type=float,\n", + " help='Energy fraction for DMD')\n", + "parser.add_argument('-rdim','--rdim',\n", + " action='store',default=-1,type=int,\n", + " help='Reduced dimension for DMD')\n", + "parser.add_argument('-crbf','--crbf',\n", + " action='store',default=0.9,type=float,\n", + " help='Closest RBF value')\n", + "parser.add_argument('-nonunif','--nonunif',dest='nonunif',\n", + " action='store_true',help='Use NonuniformDMD')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01f91b50-b46c-47d0-92d0-4936064b3d4e", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run for adaptive DMD:\n", + "\n", + "# args = parser.parse_args(\"-p 1 -rs 1 -rp 1 -o 5 -s 6 -tf 0.1 -visit\".split())\n", + "# args = parser.parse_args(\"-p 2 -rs 2 -rp 1 -o 1 -s 3 -tf 0.1 -visit\".split())\n", + "# args = parser.parse_args(\"-p 2 -rs 2 -rp 1 -o 1 -s 3 -visit\".split())\n", + "\n", + "# Sample run for nonuniform DMD:\n", + "\n", + "# args = parser.parse_args(\"-p 1 -rs 1 -rp 1 -o 5 -s 6 -tf 0.1 -nonunif -visit\".split()) \n", + "# args = parser.parse_args(\"-p 2 -rs 2 -rp 1 -o 1 -s 3 -tf 0.1 -nonunif -visit\".split()) \n", + "# args = parser.parse_args(\"-p 2 -rs 2 -rp 1 -o 1 -s 3 -nonunif -visit\".split()) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05c0a673-0043-4696-8ec2-7c88f456fc64", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "# The rest of your script\n", + "mesh = os.path.abspath(os.path.join('../data', args.mesh))\n", + "ser_ref_levels = args.refine_serial\n", + "par_ref_levels = args.refine_parallel\n", + "order = args.order\n", + "ode_solver_type = args.ode_solver\n", + "t_final = args.t_final\n", + "dt = args.time_step\n", + "cfl = args.cfl_number\n", + "visualization = args.visualization\n", + "visit = args.visit_datafiles\n", + "vis_steps = args.visualization_steps\n", + "ef = args.energy_fraction\n", + "rdim = args.rdim\n", + "crbf = args.crbf\n", + "nonunif = args.nonunif" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63aa9975-f4f7-4154-b088-72506a5f7162", + "metadata": {}, + "outputs": [], + "source": [ + "if myid == 0:\n", + " parser.print_options(args)\n", + "\n", + "device = mfem.Device('cpu')\n", + "if myid == 0:\n", + " device.Print()\n", + "\n", + "dg_euler_common.num_equation = 4\n", + "dg_euler_common.specific_heat_ratio = 1.4\n", + "dg_euler_common.gas_constant = 1.0\n", + "dg_euler_common.problem = args.problem\n", + "num_equation = dg_euler_common.num_equation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1be6c7af-bebc-4594-a846-ad04697251b0", + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Read the mesh from the given mesh file. This example requires a 2D\n", + "# periodic mesh, such as ../data/periodic-square.mesh.\n", + "meshfile = expanduser(join(dirname(mesh), '..', 'data', mesh))\n", + "mesh = mfem.Mesh(meshfile, 1, 1)\n", + "dim = mesh.Dimension()\n", + "\n", + "# 4. Define the ODE solver used for time integration. Several explicit\n", + "# Runge-Kutta methods are available.\n", + "ode_solver = None\n", + "if ode_solver_type == 1:\n", + " ode_solver = mfem.ForwardEulerSolver()\n", + "elif ode_solver_type == 2:\n", + " ode_solver = mfem.RK2Solver(1.0)\n", + "elif ode_solver_type == 3:\n", + " ode_solver = mfem.RK3SSolver()\n", + "elif ode_solver_type == 4:\n", + " ode_solver = mfem.RK4Solver()\n", + "elif ode_solver_type == 6:\n", + " ode_solver = mfem.RK6Solver()\n", + "else:\n", + " print(\"Unknown ODE solver type: \" + str(ode_solver_type))\n", + " exit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6734fecb-26ca-47ec-8bc3-2c6ca8d2a153", + "metadata": {}, + "outputs": [], + "source": [ + "# 5. Refine the mesh in serial to increase the resolution. In this example\n", + "# we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is\n", + "# a command-line parameter.\n", + "for lev in range(ser_ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "701e36df-2b80-4bc2-81e9-297230768ea8", + "metadata": {}, + "outputs": [], + "source": [ + "# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine\n", + "# this mesh further in parallel to increase the resolution. Once the\n", + "# parallel mesh is defined, the serial mesh can be deleted.\n", + "\n", + "pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh)\n", + "del mesh\n", + "for lev in range(par_ref_levels):\n", + " pmesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b75b4260-72f1-4ead-9310-d10701dce84b", + "metadata": {}, + "outputs": [], + "source": [ + "# 7. Define the discontinuous DG finite element space of the given\n", + "# polynomial order on the refined mesh.\n", + "fec = mfem.DG_FECollection(order, dim)\n", + "# Finite element space for a scalar (thermodynamic quantity)\n", + "fes = mfem.ParFiniteElementSpace(pmesh, fec)\n", + "# Finite element space for a mesh-dim vector quantity (momentum)\n", + "dfes = mfem.ParFiniteElementSpace(pmesh, fec, dim, mfem.Ordering.byNODES)\n", + "# Finite element space for all variables together (total thermodynamic state)\n", + "vfes = mfem.ParFiniteElementSpace(\n", + " pmesh, fec, num_equation, mfem.Ordering.byNODES)\n", + "\n", + "assert fes.GetOrdering() == mfem.Ordering.byNODES, \"Ordering must be byNODES\"\n", + "glob_size = vfes.GlobalTrueVSize()\n", + "if myid == 0:\n", + " print(\"Number of unknowns: \" + str(glob_size))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe303417-3440-4bf2-a3eb-c1f606e1a3d0", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Define the initial conditions, save the corresponding mesh and grid\n", + "# functions to a file. This can be opened with GLVis with the -gc option.\n", + "# The solution u has components {density, x-momentum, y-momentum, energy}.\n", + "# These are stored contiguously in the BlockVector u_block.\n", + "\n", + "offsets = [k*vfes.GetNDofs() for k in range(num_equation+1)]\n", + "offsets = mfem.intArray(offsets)\n", + "u_block = mfem.BlockVector(offsets)\n", + "\n", + "# Momentum grid function on dfes for visualization.\n", + "mom = mfem.ParGridFunction(dfes, u_block, offsets[1])\n", + "\n", + "# Initialize the state.\n", + "u0 = InitialCondition(num_equation)\n", + "sol = mfem.ParGridFunction(vfes, u_block.GetData())\n", + "sol.ProjectCoefficient(u0)\n", + "\n", + "smyid = '{:0>6d}'.format(myid)\n", + "pmesh.Print(\"vortex-mesh.\"+smyid, 8)\n", + "for k in range(num_equation):\n", + " uk = mfem.ParGridFunction(fes, u_block.GetBlock(k).GetData())\n", + " sol_name = \"vortex-\" + str(k) + \"-init.\"+smyid\n", + " uk.Save(sol_name, 8)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5513d593-574e-4af6-8ed6-f07f34588cce", + "metadata": {}, + "outputs": [], + "source": [ + "# 9. Set up the nonlinear form corresponding to the DG discretization of the\n", + "# flux divergence, and assemble the corresponding mass matrix.\n", + "Aflux = mfem.MixedBilinearForm(dfes, fes)\n", + "Aflux.AddDomainIntegrator(DomainIntegrator(dim))\n", + "Aflux.Assemble()\n", + "\n", + "A = mfem.ParNonlinearForm(vfes)\n", + "rsolver = RiemannSolver()\n", + "ii = FaceIntegrator(rsolver, dim)\n", + "A.AddInteriorFaceIntegrator(ii)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05c371cc-375d-4d44-84f6-da698f45b444", + "metadata": {}, + "outputs": [], + "source": [ + "# 10. Define the time-dependent evolution operator describing the ODE\n", + "# right-hand side, and perform time-integration (looping over the time\n", + "# iterations, ti, with a time-step dt).\n", + "euler = FE_Evolution(vfes, A, Aflux.SpMat())\n", + "\n", + "if (visualization):\n", + " MPI.COMM_WORLD.Barrier()\n", + " sout = mfem.socketstream(\"localhost\", 19916)\n", + " sout.send_text(\"parallel \" + str(num_procs) + \" \" + str(myid))\n", + " sout.precision(8)\n", + " sout.send_solution(pmesh, mom)\n", + " sout.send_text(\"pause\")\n", + " sout.flush()\n", + " if myid == 0:\n", + " print(\"GLVis visualization paused.\")\n", + " print(\" Press space (in the GLVis window) to resume it.\")\n", + "visit_dc = mfem.VisItDataCollection('DG_Euler',pmesh)\n", + "visit_dc.RegisterField('solution',mom)\n", + "if visit:\n", + " visit_dc.SetCycle(0)\n", + " visit_dc.SetTime(0.0)\n", + " visit_dc.Save()\n", + "\n", + "# Determine the minimum element size.\n", + "my_hmin = 0\n", + "if (cfl > 0):\n", + " my_hmin = min([pmesh.GetElementSize(i, 1) for i in range(pmesh.GetNE())])\n", + "hmin = MPI.COMM_WORLD.allreduce(my_hmin, op=MPI.MIN)\n", + "\n", + "# initialize timers\n", + "fom_timer, dmd_training_timer, dmd_prediction_timer = \\\n", + " StopWatch(), StopWatch(), StopWatch()\n", + "\n", + "fom_timer.Start()\n", + "t = 0.0\n", + "ts = []\n", + "euler.SetTime(t)\n", + "ode_solver.Init(euler)\n", + "fom_timer.Stop()\n", + "\n", + "if (cfl > 0):\n", + " # Find a safe dt, using a temporary vector. Calling Mult() computes the\n", + " # maximum char speed at all quadrature points on all faces.\n", + " z = mfem.Vector(A.Width())\n", + " A.Mult(sol, z)\n", + " max_char_speed = MPI.COMM_WORLD.allreduce(dg_euler_common.max_char_speed, op=MPI.MAX) \n", + " dg_euler_common.max_char_speed = max_char_speed \n", + " dt = cfl * hmin / dg_euler_common.max_char_speed / (2*order+1)\n", + "\n", + "#- DMD setup\n", + "# Initialize dmd_vars = [dmd_dens, dmd_x_mom, dmd_y_mom, dmd_e]\n", + "dmd_training_timer.Start()\n", + "if nonunif:\n", + " dmd_vars = [NonuniformDMD(u_block.GetBlock(i).Size()) \n", + " for i in range(4)]\n", + "else:\n", + " dmd_vars = [AdaptiveDMD(u_block.GetBlock(i).Size(),dt,'G','LS',crbf) \\\n", + " for i in range(4)]\n", + "for i in range(4):\n", + " dmd_vars[i].takeSample(u_block.GetBlock(i).GetDataArray(),t)\n", + "ts += [t]\n", + "dmd_training_timer.Stop()\n", + "\n", + "# Integrate in time.\n", + "done = False\n", + "ti = 0\n", + "while not done:\n", + "\n", + " fom_timer.Start()\n", + "\n", + " dt_real = min(dt, t_final - t)\n", + " t, dt_real = ode_solver.Step(sol, t, dt_real)\n", + "\n", + " if (cfl > 0):\n", + " max_char_speed = MPI.COMM_WORLD.allreduce(\n", + " dg_euler_common.max_char_speed, op=MPI.MAX)\n", + " dg_euler_common.max_char_speed = max_char_speed\n", + " dt = cfl * hmin / dg_euler_common.max_char_speed / (2*order+1)\n", + "\n", + " ti = ti+1\n", + " done = (t >= t_final - 1e-8*dt)\n", + "\n", + " fom_timer.Stop()\n", + " \n", + " #- DMD take sample\n", + " dmd_training_timer.Start()\n", + " for i in range(4):\n", + " dmd_vars[i].takeSample(u_block.GetBlock(i).GetDataArray(),t)\n", + " ts.append(t)\n", + " dmd_training_timer.Stop()\n", + "\n", + "\n", + " if (done or ti % vis_steps == 0):\n", + " if myid == 0:\n", + " print(\"time step: \" + str(ti) + \", time: \" + \"{:g}\".format(t))\n", + " if (visualization):\n", + " sout.send_text(\"parallel \" + str(num_procs) + \" \" + str(myid))\n", + " sout.send_solution(pmesh, mom)\n", + " sout.flush()\n", + " if visit:\n", + " visit_dc.SetCycle(ti)\n", + " visit_dc.SetTime(t)\n", + " visit_dc.Save()\n", + "\n", + "\n", + "if myid == 0:\n", + " print(\"done\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "726039c6-f8ce-425e-8849-c90144441b0a", + "metadata": {}, + "outputs": [], + "source": [ + "# 11. Save the final solution. This output can be viewed later using GLVis:\n", + "# \"glvis -np 4 -m vortex-mesh -g vortex-1-final\".\n", + "\n", + "for k in range(num_equation):\n", + " uk = mfem.ParGridFunction(fes, u_block.GetBlock(k).GetData())\n", + " sol_name = \"vortex-\" + str(k) + \"-final.\"+smyid\n", + " uk.Save(sol_name, 8)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9c61b98-49f2-4248-ba7b-a54a508140c1", + "metadata": {}, + "outputs": [], + "source": [ + "# 12. Compute the L2 solution error summed for all components.\n", + "if (t_final == 2.0):\n", + " error = sol.ComputeLpError(2., u0)\n", + " if myid == 0:\n", + " print(\"Solution error: \" + \"{:g}\".format(error))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17ee664a-0c3e-4698-9e0d-eca03d8099f2", + "metadata": {}, + "outputs": [], + "source": [ + "# 13. Calculate the DMD modes\n", + "if myid==0 and rdim != -1 and ef != -1:\n", + " print('Both rdim and ef are set. ef will be ignored')\n", + "\n", + "dmd_training_timer.Start()\n", + "\n", + "if rdim != -1:\n", + " if myid==0:\n", + " print(f'Creating DMD with rdim: {rdim}')\n", + " for dmd_var in dmd_vars:\n", + " dmd_var.train(rdim)\n", + "elif ef != -1:\n", + " if myid == 0:\n", + " print(f'Creating DMD with energy fraction: {ef}')\n", + " for dmd_var in dmd_vars:\n", + " dmd_var.train(ef)\n", + "\n", + "dmd_training_timer.Stop()\n", + "\n", + "true_solution_vars = [u_block.GetBlock(i) for i in range(4)]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c55cd0d-e0e7-45c1-a2f6-8cec295c2b25", + "metadata": {}, + "outputs": [], + "source": [ + "# 14. Predict the state at t_final using DMD\n", + "dmd_prediction_timer.Start()\n", + "if myid == 0:\n", + " print('Predicting density, momentum, and energy using DMD')\n", + "\n", + "result_vars = [dmd_var.predict(ts[0]) for dmd_var in dmd_vars]\n", + "initial_dmd_sols = [mfem.Vector(result_var.getData(), result_var.dim()) for result_var in result_vars]\n", + "for i in range(4):\n", + " block = u_block.GetBlock(i)\n", + " block = initial_dmd_sols[i]\n", + " #u_block.Update(initial_dmd_sols[i],offsets[i])\n", + "\n", + "dmd_visit_dc = mfem.VisItDataCollection('DMD_DG_Euler', pmesh)\n", + "dmd_visit_dc.RegisterField('solution',mom)\n", + "if (visit):\n", + " dmd_visit_dc.SetCycle(0)\n", + " dmd_visit_dc.SetTime(0.0)\n", + " dmd_visit_dc.Save()\n", + "\n", + "if visit:\n", + " for i in range(len(ts)):\n", + " if (i==(len(ts)-2)) or (i%vis_steps==0):\n", + " result_vars = [var.predict(ts[i]) for var in dmd_vars]\n", + " dmd_sols = [mfem.Vector(result_var.getData(),result_var.dim()) for result_var in result_vars]\n", + " for k in range(4):\n", + " block = u_block.GetBlock(i)\n", + " block = dmd_sols[k]\n", + "\n", + " dmd_visit_dc.SetCycle(i)\n", + " dmd_visit_dc.SetTime(ts[i])\n", + " dmd_visit_dc.Save()\n", + "dmd_prediction_timer.Stop()\n", + "result_vars = [dmd_var.predict(t_final) for dmd_var in dmd_vars]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a90ddf1b-c15a-4053-9705-4142cf747dcd", + "metadata": {}, + "outputs": [], + "source": [ + "# 15. Calculate the relative error between the DMD final solution and the true solution.\n", + "dmd_solution_vars = [mfem.Vector(result_var.getData(),result_var.dim()) for result_var in result_vars]\n", + "diff_vars = [mfem.Vector(result_var.dim()) for result_var in result_vars]\n", + "for i in range(4):\n", + " mfem.subtract_vector(dmd_solution_vars[i],true_solution_vars[i],\\\n", + " diff_vars[i])\n", + "tot_diff_norm_vars = [sqrt(mfem.InnerProduct(MPI.COMM_WORLD,diff_var,diff_var)) for diff_var in diff_vars]\n", + "tot_true_solution_norm_vars = [sqrt(mfem.InnerProduct(MPI.COMM_WORLD,true_solution_var,true_solution_var))\\\n", + " for true_solution_var in true_solution_vars]\n", + "\n", + "if myid==0:\n", + " var_names = ['dens', 'x_mom', 'y_mom', 'e']\n", + " for i in range(len(var_names)):\n", + " rel_error = tot_diff_norm_vars[i]/tot_true_solution_norm_vars[i]\n", + " print(f'Relative error of DMD {var_names[i]} at t_final: {t_final} is {rel_error:.10f}')\n", + "\n", + " print(f'Elapsed time for solving FOM: {fom_timer.duration:.6e}')\n", + " print(f'Elapsed time for training DMD: {dmd_training_timer.duration:.6e}')\n", + " print(f'Elapsed time for predicting DMD: {dmd_prediction_timer.duration:.6e}')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/dmd/nonlinear_elasticity.ipynb b/examples/dmd/nonlinear_elasticity.ipynb new file mode 100644 index 0000000..1ad6f52 --- /dev/null +++ b/examples/dmd/nonlinear_elasticity.ipynb @@ -0,0 +1,799 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "84c2db4d-dbde-4732-b251-51a39e21331b", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import mfem.par as mfem\n", + "\n", + "from mfem.common.arg_parser import ArgParser\n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sqrt, pi, cos, sin, hypot, arctan2\n", + "from scipy.special import erfc\n", + "\n", + "from pylibROM.python_utils.StopWatch import StopWatch\n", + "from pylibROM.algo import DMD\n", + "\n", + "from mfem.par import intArray, add_vector\n", + "from mpi4py import MPI\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db14c3f4-10c6-4c2e-a223-c7a2e8c8e088", + "metadata": {}, + "outputs": [], + "source": [ + "num_procs = MPI.COMM_WORLD.size\n", + "myid = MPI.COMM_WORLD.rank" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d949c86b-3e42-4270-8666-2b59acc0c23d", + "metadata": {}, + "outputs": [], + "source": [ + "parser = ArgParser(description='nonlinear_elasticity')\n", + "parser.add_argument('-m', '--mesh',\n", + " default='beam-quad.mesh',\n", + " action='store', type=str,\n", + " help='Mesh file to use.')\n", + "parser.add_argument('-rs', '--refine-serial',\n", + " action='store', default=2, type=int,\n", + " help=\"Number of times to refine the mesh uniformly before parallel\")\n", + "parser.add_argument('-rp', '--refine-parallel',\n", + " action='store', default=0, type=int,\n", + " help=\"Number of times to refine the mesh uniformly after parallel\")\n", + "parser.add_argument('-o', '--order',\n", + " action='store', default=2, type=int,\n", + " help=\"Finite element order (polynomial degree)\")\n", + "help_ode = \"\\n\".join([\"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\",\n", + " \"\\t11 - Forward Euler, 12 - RK2\",\n", + " \"\\t13 - RK3 SSP, 14 - RK4.\"])\n", + "parser.add_argument('-s', '--ode-solver',\n", + " action='store', default=3, type=int,\n", + " help=help_ode)\n", + "parser.add_argument('-tf', '--t-final',\n", + " action='store', default=300.0, type=float,\n", + " help=\"Final time; start time is 0.\")\n", + "parser.add_argument('-dt', '--time-step',\n", + " action='store', default=3.0, type=float,\n", + " help=\"Time step\")\n", + "parser.add_argument(\"-v\", \"--viscosity\",\n", + " action='store', default=1e-2, type=float,\n", + " help=\"Viscosity coefficient.\")\n", + "parser.add_argument(\"-mu\", \"--shear-modulus\",\n", + " action='store', default=0.25, type=float,\n", + " help=\"Shear modulus in the Neo-Hookean hyperelastic model.\")\n", + "parser.add_argument(\"-K\", \"--bulk-modulus\",\n", + " action='store', default=5.0, type=float,\n", + " help=\"Bulk modulus in the Neo-Hookean hyperelastic model.\")\n", + "parser.add_argument('-vis', '--visualization',\n", + " action='store_true', default=True,\n", + " help='Enable GLVis visualization')\n", + "parser.add_argument('-visit', '--visit-datafiles',\n", + " action='store_true', default=False,\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-vs\", \"--visualization-steps\",\n", + " action='store', default=1, type=int,\n", + " help=\"Visualize every n-th timestep.\")\n", + "parser.add_argument('-ef','--energy_fraction',\n", + " action='store', default=0.9999, type=float,\n", + " help='Energy fraction for DMD')\n", + "parser.add_argument('-rdim','--rdim',\n", + " action='store', default=-1, type=int,\n", + " help='Reduced dimension for DMD')\n", + "parser.add_argument('-nwinsamp','--numwindowsamples',\n", + " action='store', default=sys.maxsize, type=int,\n", + " help='Number of samples in DMD windows.');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a567de4-af6a-42cc-a63d-a8fd27767121", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run for serial case:\n", + "# args = parser.parse_args(\"-s 2 -rs 1 -dt 0.01 -tf 5 -visit\".split())\n", + "\n", + "# Sample run for time windowing case:\n", + "# args = parser.parse_args(\"-s 2 -rs 1 -dt 0.01 -tf 5 -nwinsamp 10 -visit\".split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "734d3c6b-1e8c-4c84-bddc-ba2838ba27de", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "ser_ref_levels = args.refine_serial\n", + "par_ref_levels = args.refine_parallel\n", + "order = args.order\n", + "ode_solver_type = args.ode_solver\n", + "t_final = args.t_final\n", + "dt = args.time_step\n", + "visc = args.viscosity\n", + "mu = args.shear_modulus\n", + "K = args.bulk_modulus\n", + "visualization = args.visualization\n", + "visit = args.visit_datafiles\n", + "vis_steps = args.visualization_steps\n", + "ef = args.energy_fraction\n", + "rdim = args.rdim\n", + "windowNumSamples = args.numwindowsamples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96364bd5-8ddd-4e01-9dcd-45905004b9ee", + "metadata": {}, + "outputs": [], + "source": [ + "if (myid == 0):\n", + " parser.print_options(args)\n", + "\n", + "device = mfem.Device('cpu')\n", + "if myid == 0:\n", + " device.Print()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dfb6f4de-0374-4314-93b9-c0d458f7dc0d", + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Read the serial mesh from the given mesh file on all processors. We can\n", + "# handle triangular, quadrilateral, tetrahedral and hexahedral meshes\n", + "# with the same code.\n", + "meshfile = expanduser(os.path.join('../data', args.mesh))\n", + "mesh = mfem.Mesh(meshfile, 1, 1)\n", + "dim = mesh.Dimension()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2097e25-4a8e-4537-a8f1-aa9788991d27", + "metadata": {}, + "outputs": [], + "source": [ + "# 4. Define the ODE solver used for time integration. Several implicit\n", + "# singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as\n", + "# explicit Runge-Kutta methods are available.\n", + "if ode_solver_type == 1:\n", + " ode_solver = BackwardEulerSolver()\n", + "elif ode_solver_type == 2:\n", + " ode_solver = mfem.SDIRK23Solver(2)\n", + "elif ode_solver_type == 3:\n", + " ode_solver = mfem.SDIRK33Solver()\n", + "elif ode_solver_type == 11:\n", + " ode_solver = ForwardEulerSolver()\n", + "elif ode_solver_type == 12:\n", + " ode_solver = mfem.RK2Solver(0.5)\n", + "elif ode_solver_type == 13:\n", + " ode_solver = mfem.RK3SSPSolver()\n", + "elif ode_solver_type == 14:\n", + " ode_solver = mfem.RK4Solver()\n", + "elif ode_solver_type == 22:\n", + " ode_solver = mfem.ImplicitMidpointSolver()\n", + "elif ode_solver_type == 23:\n", + " ode_solver = mfem.SDIRK23Solver()\n", + "elif ode_solver_type == 24:\n", + " ode_solver = mfem.SDIRK34Solver()\n", + "else:\n", + " if myid == 0:\n", + " print(\"Unknown ODE solver type: \" + str(ode_solver_type))\n", + " sys.exit()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "83608d9d-2c8f-4b4b-94bd-0e1b0f662916", + "metadata": {}, + "outputs": [], + "source": [ + "# 5. Refine the mesh in serial to increase the resolution. In this example\n", + "# we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is\n", + "# a command-line parameter.\n", + "for lev in range(ser_ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dba97b4e-9e95-4248-ba81-9ded94e8fb6b", + "metadata": {}, + "outputs": [], + "source": [ + "# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine\n", + "# this mesh further in parallel to increase the resolution. Once the\n", + "# parallel mesh is defined, the serial mesh can be deleted.\n", + "pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh)\n", + "del mesh\n", + "for x in range(par_ref_levels):\n", + " pmesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41caa721-b6e8-4c73-8401-d8ef1bca6d10", + "metadata": {}, + "outputs": [], + "source": [ + "# 7. Define the parallel vector finite element spaces representing the mesh\n", + "# deformation x_gf, the velocity v_gf, and the initial configuration,\n", + "# x_ref. Define also the elastic energy density, w_gf, which is in a\n", + "# discontinuous higher-order space. Since x and v are integrated in time\n", + "# as a system, we group them together in block vector vx, on the unique\n", + "# parallel degrees of freedom, with offsets given by array true_offset.\n", + "\n", + "fec = mfem.H1_FECollection(order, dim)\n", + "fespace = mfem.ParFiniteElementSpace(pmesh, fec, dim)\n", + "glob_size = fespace.GlobalTrueVSize()\n", + "if (myid == 0):\n", + " print('Number of velocity/deformation unknowns: ' + str(glob_size))\n", + "\n", + "true_size = fespace.TrueVSize()\n", + "true_offset = mfem.intArray(3)\n", + "true_offset[0] = 0\n", + "true_offset[1] = true_size\n", + "true_offset[2] = 2*true_size\n", + "\n", + "vx = mfem.BlockVector(true_offset)\n", + "\n", + "v_gf = mfem.ParGridFunction(fespace)\n", + "x_gf = mfem.ParGridFunction(fespace)\n", + "\n", + "x_ref = mfem.ParGridFunction(fespace)\n", + "pmesh.GetNodes(x_ref)\n", + "\n", + "w_fec = mfem.L2_FECollection(order + 1, dim)\n", + "w_fespace = mfem.ParFiniteElementSpace(pmesh, w_fec)\n", + "w_gf = mfem.ParGridFunction(w_fespace)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1188f2c6-50c2-47b6-a982-66d15cc0614a", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Set the initial conditions for v_gf, x_gf and vx, and define the\n", + "# boundary conditions on a beam-like mesh (see description above).\n", + "\n", + "\n", + "class InitialVelocity(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " dim = len(x)\n", + " s = 0.1/64.\n", + "\n", + " v = np.zeros(len(x))\n", + " v[-1] = s*x[0]**2*(8.0-x[0])\n", + " v[0] = -s*x[0]**2\n", + " return v\n", + "\n", + "\n", + "class InitialDeformation(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " return x.copy()\n", + "\n", + "\n", + "velo = InitialVelocity(dim)\n", + "v_gf.ProjectCoefficient(velo)\n", + "deform = InitialDeformation(dim)\n", + "x_gf.ProjectCoefficient(deform)\n", + "\n", + "v_gf.GetTrueDofs(vx.GetBlock(0))\n", + "x_gf.GetTrueDofs(vx.GetBlock(1))\n", + "\n", + "ess_bdr = mfem.intArray(fespace.GetMesh().bdr_attributes.Max())\n", + "ess_bdr.Assign(0)\n", + "ess_bdr[0] = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e64b9de5-9906-446c-9888-f09b78c22f06", + "metadata": {}, + "outputs": [], + "source": [ + "# 9. Initialize the hyperelastic operator, the GLVis visualization and print\n", + "# the initial energies.\n", + "class ElasticEnergyCoefficient(mfem.PyCoefficient):\n", + " def __init__(self, model, x):\n", + " self.x = x\n", + " self.model = model\n", + " self.J = mfem.DenseMatrix()\n", + " mfem.PyCoefficient.__init__(self)\n", + "\n", + " def Eval(self, T, ip):\n", + " self.model.SetTransformation(T)\n", + " self.x.GetVectorGradient(T, self.J)\n", + " return self.model.EvalW(self.J)/(self.J.Det())\n", + "\n", + "\n", + "class ReducedSystemOperator(mfem.PyOperator):\n", + " def __init__(self, M, S, H, ess_tdof_list):\n", + " mfem.PyOperator.__init__(self, M.ParFESpace().TrueVSize())\n", + " self.M = M\n", + " self.S = S\n", + " self.H = H\n", + " self.Jacobian = None\n", + " h = M.ParFESpace().TrueVSize()\n", + " self.w = mfem.Vector(h)\n", + " self.z = mfem.Vector(h)\n", + " self.dt = 0.0\n", + " self.ess_tdof_list = ess_tdof_list\n", + "\n", + " def SetParameters(self, dt, v, x):\n", + " self.dt = dt\n", + " self.v = v\n", + " self.x = x\n", + "\n", + " def Mult(self, k, y):\n", + " add_vector(self.v, self.dt, k, self.w)\n", + " add_vector(self.x, self.dt, self.w, self.z)\n", + " self.H.Mult(self.z, y)\n", + " self.M.TrueAddMult(k, y)\n", + " self.S.TrueAddMult(self.w, y)\n", + " y.SetSubVector(self.ess_tdof_list, 0.0)\n", + "\n", + " def GetGradient(self, k):\n", + " localJ = mfem.Add(1.0, self.M.SpMat(), self.dt, self.S.SpMat())\n", + " add_vector(self.v, self.dt, k, self.w)\n", + " add_vector(self.x, self.dt, self.w, self.z)\n", + " localJ.Add(self.dt * self.dt, self.H.GetLocalGradient(self.z))\n", + " Jacobian = self.M.ParallelAssemble(localJ)\n", + " Jacobian.EliminateRowsCols(self.ess_tdof_list)\n", + " return Jacobian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9106bae9-a3f2-40e5-a1d8-115c7634ae9f", + "metadata": {}, + "outputs": [], + "source": [ + "class HyperelasticOperator(mfem.PyTimeDependentOperator):\n", + " def __init__(self, fespace, ess_bdr, visc, mu, K):\n", + " mfem.PyTimeDependentOperator.__init__(self, 2*fespace.TrueVSize(), 0.0)\n", + "\n", + " rel_tol = 1e-8\n", + " skip_zero_entries = 0\n", + " ref_density = 1.0\n", + "\n", + " self.ess_tdof_list = intArray()\n", + " self.z = mfem.Vector(self.Height()//2)\n", + " self.fespace = fespace\n", + " self.viscosity = visc\n", + " self.newton_solver = mfem.NewtonSolver(fespace.GetComm())\n", + "\n", + " M = mfem.ParBilinearForm(fespace)\n", + " S = mfem.ParBilinearForm(fespace)\n", + " H = mfem.ParNonlinearForm(fespace)\n", + " self.M = M\n", + " self.H = H\n", + " self.S = S\n", + "\n", + " rho = mfem.ConstantCoefficient(ref_density)\n", + " M.AddDomainIntegrator(mfem.VectorMassIntegrator(rho))\n", + " M.Assemble(skip_zero_entries)\n", + " M.Finalize(skip_zero_entries)\n", + " self.Mmat = M.ParallelAssemble()\n", + "\n", + " fespace.GetEssentialTrueDofs(ess_bdr, self.ess_tdof_list)\n", + " self.Mmat.EliminateRowsCols(self.ess_tdof_list)\n", + "\n", + " M_solver = mfem.CGSolver(fespace.GetComm())\n", + " M_prec = mfem.HypreSmoother()\n", + " M_solver.iterative_mode = False\n", + " M_solver.SetRelTol(rel_tol)\n", + " M_solver.SetAbsTol(0.0)\n", + " M_solver.SetMaxIter(30)\n", + " M_solver.SetPrintLevel(0)\n", + " M_prec.SetType(mfem.HypreSmoother.Jacobi)\n", + " M_solver.SetPreconditioner(M_prec)\n", + " M_solver.SetOperator(self.Mmat)\n", + "\n", + " self.M_solver = M_solver\n", + " self.M_prec = M_prec\n", + "\n", + " model = mfem.NeoHookeanModel(mu, K)\n", + " H.AddDomainIntegrator(mfem.HyperelasticNLFIntegrator(model))\n", + " H.SetEssentialTrueDofs(self.ess_tdof_list)\n", + " self.model = model\n", + "\n", + " visc_coeff = mfem.ConstantCoefficient(visc)\n", + " S.AddDomainIntegrator(mfem.VectorDiffusionIntegrator(visc_coeff))\n", + " S.Assemble(skip_zero_entries)\n", + " S.Finalize(skip_zero_entries)\n", + "\n", + " self.reduced_oper = ReducedSystemOperator(M, S, H, self.ess_tdof_list)\n", + "\n", + " J_hypreSmoother = mfem.HypreSmoother()\n", + " J_hypreSmoother.SetType(mfem.HypreSmoother.l1Jacobi)\n", + " J_hypreSmoother.SetPositiveDiagonal(True)\n", + " J_prec = J_hypreSmoother\n", + "\n", + " J_minres = mfem.MINRESSolver(fespace.GetComm())\n", + " J_minres.SetRelTol(rel_tol)\n", + " J_minres.SetAbsTol(0.0)\n", + " J_minres.SetMaxIter(300)\n", + " J_minres.SetPrintLevel(-1)\n", + " J_minres.SetPreconditioner(J_prec)\n", + "\n", + " self.J_solver = J_minres\n", + " self.J_prec = J_prec\n", + "\n", + " newton_solver = mfem.NewtonSolver(fespace.GetComm())\n", + " newton_solver.iterative_mode = False\n", + " newton_solver.SetSolver(self.J_solver)\n", + " newton_solver.SetOperator(self.reduced_oper)\n", + " newton_solver.SetPrintLevel(1) # print Newton iterations\n", + " newton_solver.SetRelTol(rel_tol)\n", + " newton_solver.SetAbsTol(0.0)\n", + " newton_solver.SetAdaptiveLinRtol(2, 0.5, 0.9)\n", + " newton_solver.SetMaxIter(10)\n", + " self.newton_solver = newton_solver\n", + "\n", + " def Mult(self, vx, vx_dt):\n", + " sc = self.Height()//2\n", + " v = mfem.Vector(vx, 0, sc)\n", + " x = mfem.Vector(vx, sc, sc)\n", + " dv_dt = mfem.Vector(dvx_dt, 0, sc)\n", + " dx_dt = mfem.Vector(dvx_dt, sc, sc)\n", + "\n", + " self.H.Mult(x, z)\n", + " if (self.viscosity != 0.0):\n", + " S.TrueAddMult(v, z)\n", + " z.SetSubVector(self.ess_tdof_list, 0.0)\n", + " z.Neg()\n", + " self.M_solver.Mult(z, dv_dt)\n", + " dx_dt = v\n", + "\n", + " def ImplicitSolve(self, dt, vx, dvx_dt):\n", + " sc = self.Height()//2\n", + " v = mfem.Vector(vx, 0, sc)\n", + " x = mfem.Vector(vx, sc, sc)\n", + " dv_dt = mfem.Vector(dvx_dt, 0, sc)\n", + " dx_dt = mfem.Vector(dvx_dt, sc, sc)\n", + "\n", + " # By eliminating kx from the coupled system:\n", + " # kv = -M^{-1}*[H(x + dt*kx) + S*(v + dt*kv)]\n", + " # kx = v + dt*kv\n", + " # we reduce it to a nonlinear equation for kv, represented by the\n", + " # backward_euler_oper. This equation is solved with the newton_solver\n", + " # object (using J_solver and J_prec internally).\n", + " self.reduced_oper.SetParameters(dt, v, x)\n", + " zero = mfem.Vector() # empty vector is interpreted as\n", + " # zero r.h.s. by NewtonSolver\n", + " self.newton_solver.Mult(zero, dv_dt)\n", + " add_vector(v, dt, dv_dt, dx_dt)\n", + "\n", + " def ElasticEnergy(self, x):\n", + " return self.H.GetEnergy(x)\n", + "\n", + " def KineticEnergy(self, v):\n", + " local_energy = 0.5*self.M.InnerProduct(v, v)\n", + " energy = MPI.COMM_WORLD.allreduce(local_energy, op=MPI.SUM)\n", + " return energy\n", + "\n", + " def GetElasticEnergyDensity(self, x, w):\n", + " w_coeff = ElasticEnergyCoefficient(self.model, x)\n", + " w.ProjectCoefficient(w_coeff)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d49e8bd7-46ea-4fb8-b4ff-85ba18bc0fbd", + "metadata": {}, + "outputs": [], + "source": [ + "def visualize(out, pmesh, deformed_nodes, field,\n", + " field_name='', init_vis=False):\n", + " nodes = deformed_nodes\n", + " owns_nodes = 0\n", + "\n", + " nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " out.send_text(\"parallel \" + str(num_procs) + \" \" + str(myid))\n", + " out.send_solution(pmesh, field)\n", + "\n", + " nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " if (init_vis):\n", + " out.send_text(\"window_size 400 400\")\n", + " out.send_text(\"window_title '\" + field_name)\n", + " if (pmesh.SpaceDimension() == 2):\n", + " out.send_text(\"view 0 0\")\n", + " out.send_text(\"keys jl\")\n", + " out.send_text(\"keys cm\") # show colorbar and mesh\n", + " # update value-range; keep mesh-extents fixed\n", + " out.send_text(\"autoscale value\")\n", + " out.send_text(\"pause\")\n", + " out.flush()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eeba1602-237c-4a93-9ef3-47963bc0410d", + "metadata": {}, + "outputs": [], + "source": [ + "oper = HyperelasticOperator(fespace, ess_bdr, visc, mu, K)\n", + "if (visualization):\n", + " vis_v = mfem.socketstream(\"localhost\", 19916)\n", + " vis_v.precision(8)\n", + " visualize(vis_v, pmesh, x_gf, v_gf, \"Velocity\", True)\n", + "\n", + " MPI.COMM_WORLD.Barrier()\n", + " vis_w = mfem.socketstream(\"localhost\", 19916)\n", + " oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + " vis_w.precision(8)\n", + " visualize(vis_w, pmesh, x_gf, w_gf, \"Elastic energy density\", True)\n", + "\n", + "ee0 = oper.ElasticEnergy(x_gf)\n", + "ke0 = oper.KineticEnergy(v_gf)\n", + "\n", + "if myid == 0:\n", + " print(\"initial elastic energy (EE) = \" + str(ee0))\n", + " print(\"initial kinetic energy (KE) = \" + str(ke0))\n", + " print(\"initial total energy (TE) = \" + str(ee0 + ke0))\n", + "\n", + "# initialize timers\n", + "fom_timer, dmd_training_timer, dmd_prediction_timer = \\\n", + " StopWatch(), StopWatch(), StopWatch()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04673fcb-cdd3-4f93-b294-982e6a87ac30", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Perform time-integration (looping over the time iterations, ti, with a\n", + "# time-step dt).\n", + "t = 0.\n", + "ti = 1\n", + "fom_timer.Start()\n", + "oper.SetTime(t)\n", + "ode_solver.Init(oper)\n", + "last_step = False\n", + "fom_timer.Stop()\n", + "\n", + "curr_window = 0\n", + "dmd_v = []\n", + "dmd_x = []\n", + "ts = []\n", + "dmd_training_timer.Start()\n", + "dmd_v.append(DMD(vx.GetBlock(0).Size(),dt))\n", + "dmd_x.append(DMD(vx.GetBlock(1).Size(),dt))\n", + "dmd_v[curr_window].takeSample(vx.GetBlock(0).GetDataArray(), t)\n", + "dmd_x[curr_window].takeSample(vx.GetBlock(1).GetDataArray(), t)\n", + "ts.append(t)\n", + "dmd_training_timer.Stop()\n", + "\n", + "while not last_step:\n", + " fom_timer.Start()\n", + " dt_real = min(dt, t_final - t)\n", + " t, dt = ode_solver.Step(vx, t, dt_real)\n", + " if (t >= t_final - 1e-8*dt):\n", + " last_step = True\n", + " fom_timer.Stop()\n", + "\n", + " dmd_training_timer.Start()\n", + " dmd_v[curr_window].takeSample(vx.GetBlock(0).GetDataArray(), t)\n", + " dmd_x[curr_window].takeSample(vx.GetBlock(1).GetDataArray(), t)\n", + " if (last_step or (ti % windowNumSamples) == 0):\n", + " if myid == 0 and rdim != -1 and ef != -1:\n", + " print('Both rdim and ef are set. ef will be ignored')\n", + " if rdim != -1:\n", + " if myid == 0:\n", + " print(f'Creating DMD with rdim: {rdim}')\n", + " dmd_v[curr_window].train(rdim)\n", + " dmd_x[curr_window].train(rdim)\n", + " elif ef != 1:\n", + " if myid == 0:\n", + " print(f'Creating DMD with energy fraction: {ef}')\n", + " dmd_v[curr_window].train(ef)\n", + " dmd_x[curr_window].train(ef)\n", + " if not last_step:\n", + " curr_window += 1\n", + " dmd_v.append(DMD(vx.GetBlock(0).Size(),dt))\n", + " dmd_x.append(DMD(vx.GetBlock(1).Size(),dt))\n", + " dmd_v[curr_window].takeSample(vx.GetBlock(0).GetDataArray(), t)\n", + " dmd_x[curr_window].takeSample(vx.GetBlock(1).GetDataArray(), t)\n", + " ts.append(t)\n", + " dmd_training_timer.Stop()\n", + "\n", + " if (last_step or (ti % vis_steps) == 0):\n", + " v_gf.Distribute(vx.GetBlock(0))\n", + " x_gf.Distribute(vx.GetBlock(1))\n", + "\n", + " ee = oper.ElasticEnergy(x_gf)\n", + " ke = oper.KineticEnergy(v_gf)\n", + "\n", + " text = (\"step \" + str(ti) + \", t = \" + str(t) +\n", + " \", EE = \" + \"{:g}\".format(ee) +\n", + " \", KE = \" + \"{:g}\".format(ke) +\n", + " \", dTE = \" + \"{:g}\".format((ee+ke)-(ee0+ke0)))\n", + "\n", + " if myid == 0:\n", + " print(text)\n", + " if visualization:\n", + " visualize(vis_v, pmesh, x_gf, v_gf)\n", + " oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + " visualize(vis_w, pmesh, x_gf, w_gf)\n", + "\n", + " ti = ti + 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a059fafc-b2cc-4463-a91c-599e6185c27a", + "metadata": {}, + "outputs": [], + "source": [ + "true_solution_v = vx.GetBlock(0)\n", + "true_solution_x = vx.GetBlock(1)\n", + "\n", + "# Predict using DMD.\n", + "dmd_prediction_timer.Start()\n", + "if myid == 0:\n", + " print('Predicting position and velocity using DMD')\n", + "\n", + "curr_window = 0\n", + "result_v = dmd_v[curr_window].predict(ts[0]) \n", + "result_x = dmd_x[curr_window].predict(ts[0]) \n", + "initial_dmd_v = mfem.Vector(result_v.getData(), result_v.dim())\n", + "initial_dmd_x = mfem.Vector(result_x.getData(), result_x.dim())\n", + "v_gf.SetFromTrueDofs(initial_dmd_v)\n", + "x_gf.SetFromTrueDofs(initial_dmd_x)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b725f51-f7f9-43ed-92d3-0e3a44aac0ea", + "metadata": {}, + "outputs": [], + "source": [ + "if visit:\n", + " dmd_dc = mfem.VisItDataCollection('DMD_nonlinear_elasticity', pmesh)\n", + " dmd_dc.RegisterField('v',v_gf)\n", + " dmd_dc.RegisterField('x',x_gf)\n", + " dmd_dc.SetCycle(0)\n", + " dmd_dc.SetTime(0.0)\n", + " dmd_dc.Save()\n", + "\n", + " for i in range(1,len(ts)):\n", + " result_v = dmd_v[curr_window].predict(ts[i]) \n", + " result_x = dmd_x[curr_window].predict(ts[i]) \n", + " dmd_solution_v = mfem.Vector(result_v.getData(), result_v.dim())\n", + " dmd_solution_x = mfem.Vector(result_x.getData(), result_x.dim())\n", + " v_gf.SetFromTrueDofs(dmd_solution_v)\n", + " x_gf.SetFromTrueDofs(dmd_solution_x)\n", + "\n", + " nodes = x_gf\n", + " owns_nodes = 0\n", + " nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)\n", + " dmd_dc.SetCycle(i)\n", + " dmd_dc.SetTime(ts[i])\n", + " dmd_dc.Save()\n", + " nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " if i % windowNumSamples == 0 and i < len(ts)-1:\n", + " curr_window += 1\n", + "else:\n", + " curr_window = len(dmd_v) - 1\n", + "\n", + "result_v = dmd_v[curr_window].predict(t_final) \n", + "result_x = dmd_x[curr_window].predict(t_final) \n", + "\n", + "dmd_prediction_timer.Stop()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f7d31d8-3e6e-4734-b253-29746b1c0765", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate the relative error between the DMD final solution and the true solution.\n", + "dmd_solution_v = mfem.Vector(result_v.getData(), result_v.dim())\n", + "diff_v = mfem.Vector(true_solution_v.Size())\n", + "mfem.subtract_vector(dmd_solution_v, true_solution_v, diff_v)\n", + "tot_diff_norm_v = np.sqrt(mfem.InnerProduct(diff_v, diff_v))\n", + "tot_true_solution_v_norm = np.sqrt(mfem.InnerProduct(true_solution_v, true_solution_v))\n", + "\n", + "dmd_solution_x = mfem.Vector(result_x.getData(), result_x.dim())\n", + "diff_x = mfem.Vector(true_solution_x.Size())\n", + "mfem.subtract_vector(dmd_solution_x, true_solution_x, diff_x)\n", + "tot_diff_norm_x = np.sqrt(mfem.InnerProduct(diff_x, diff_x))\n", + "tot_true_solution_x_norm = np.sqrt(mfem.InnerProduct(true_solution_x, true_solution_x))\n", + "\n", + "if myid == 0:\n", + " print(\"Relative error of DMD velocity (v) at t_final: %f is %.3E\" % (t_final, tot_diff_norm_v / tot_true_solution_v_norm))\n", + " print(\"Relative error of DMD position (x) at t_final: %f is %.3E\" % (t_final, tot_diff_norm_x / tot_true_solution_x_norm))\n", + " print(\"Elapsed time for solving FOM: %e second\\n\" % fom_timer.duration)\n", + " print(\"Elapsed time for training DMD: %e second\\n\" % dmd_training_timer.duration)\n", + " print(\"Elapsed time for predicting DMD: %e second\\n\" % dmd_prediction_timer.duration)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d19e06ac-67e5-491b-b79a-19765e09bb80", + "metadata": {}, + "outputs": [], + "source": [ + "#\n", + "# if i translate c++ line-by-line, ti seems the second swap does not work...\n", + "#\n", + "\n", + "smyid = '{:0>6d}'.format(myid)\n", + "mesh_name = \"deformed.\"+smyid\n", + "velo_name = \"velocity.\"+smyid\n", + "ee_name = \"elastic_energy.\"+smyid\n", + "\n", + "nodes = x_gf\n", + "owns_nodes = 0\n", + "nodes, owns_nodes = pmesh.SwapNodes(nodes, owns_nodes)\n", + "pmesh.Print(mesh_name, 8)\n", + "pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + "v_gf.Save(velo_name, 8)\n", + "oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + "w_gf.Save(ee_name, 8)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/dmd/parametric_heat_conduction.ipynb b/examples/dmd/parametric_heat_conduction.ipynb new file mode 100644 index 0000000..e9d2382 --- /dev/null +++ b/examples/dmd/parametric_heat_conduction.ipynb @@ -0,0 +1,918 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "b7677341-3c79-4d04-b46b-9ab2398c95dc", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import io\n", + "import pathlib\n", + "import sys\n", + "try:\n", + " import mfem.par as mfem\n", + "except ModuleNotFoundError:\n", + " msg = \"PyMFEM is not installed yet. Install PyMFEM:\\n\"\n", + " msg += \"\\tgit clone https://github.com/mfem/PyMFEM.git\\n\"\n", + " msg += \"\\tcd PyMFEM\\n\"\n", + " msg += \"\\tpython3 setup.py install --with-parallel\\n\"\n", + " raise ModuleNotFoundError(msg)\n", + "\n", + "from mfem.par import intArray\n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum\n", + "\n", + "sys.path.append(\"../../build\")\n", + "import pylibROM.algo as algo\n", + "import pylibROM.linalg as linalg\n", + "import pylibROM.utils as utils\n", + "from pylibROM.python_utils import StopWatch\n", + "from pylibROM.mfem import PointwiseSnapshot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc785e8f-b758-4459-bdae-f957f7ff3cad", + "metadata": {}, + "outputs": [], + "source": [ + "class ConductionOperator(mfem.PyTimeDependentOperator):\n", + " def __init__(self, fespace, alpha, kappa, u):\n", + " mfem.PyTimeDependentOperator.__init__(\n", + " self, fespace.GetTrueVSize(), 0.0)\n", + " rel_tol = 1e-8\n", + " self.alpha = alpha\n", + " self.kappa = kappa\n", + " self.T = None\n", + " self.K = None\n", + " self.M = None\n", + " self.fespace = fespace\n", + "\n", + " self.ess_tdof_list = intArray()\n", + " self.Mmat = mfem.HypreParMatrix()\n", + " self.Kmat = mfem.HypreParMatrix()\n", + " self.M_solver = mfem.CGSolver(fespace.GetComm())\n", + " self.M_prec = mfem.HypreSmoother()\n", + " self.T_solver = mfem.CGSolver(fespace.GetComm())\n", + " self.T_prec = mfem.HypreSmoother()\n", + " self.z = mfem.Vector(self.Height())\n", + "\n", + " self.M = mfem.ParBilinearForm(fespace)\n", + " self.M.AddDomainIntegrator(mfem.MassIntegrator())\n", + " self.M.Assemble()\n", + " self.M.FormSystemMatrix(self.ess_tdof_list, self.Mmat)\n", + "\n", + " self.M_solver.iterative_mode = False\n", + " self.M_solver.SetRelTol(rel_tol)\n", + " self.M_solver.SetAbsTol(0.0)\n", + " self.M_solver.SetMaxIter(100)\n", + " self.M_solver.SetPrintLevel(0)\n", + " self.M_prec.SetType(mfem.HypreSmoother.Jacobi)\n", + " self.M_solver.SetPreconditioner(self.M_prec)\n", + " self.M_solver.SetOperator(self.Mmat)\n", + "\n", + " self.T_solver.iterative_mode = False\n", + " self.T_solver.SetRelTol(rel_tol)\n", + " self.T_solver.SetAbsTol(0.0)\n", + " self.T_solver.SetMaxIter(100)\n", + " self.T_solver.SetPrintLevel(0)\n", + " self.T_solver.SetPreconditioner(self.T_prec)\n", + "\n", + " self.SetParameters(u)\n", + "\n", + " def Mult(self, u, u_dt):\n", + " # Compute:\n", + " # du_dt = M^{-1}*-K(u) for du_dt\n", + " self.Kmat.Mult(u, self.z)\n", + " self.z.Neg() # z = -z\n", + " self.M_solver.Mult(self.z, du_dt)\n", + "\n", + " def ImplicitSolve(self, dt, u, du_dt):\n", + " # Solve the equation:\n", + " # du_dt = M^{-1}*[-K(u + dt*du_dt)]\n", + " # for du_dt\n", + " if self.T is None:\n", + " self.T = mfem.Add(1.0, self.Mmat, dt, self.Kmat)\n", + " current_dt = dt\n", + " self.T_solver.SetOperator(self.T)\n", + " self.Kmat.Mult(u, self.z)\n", + " self.z.Neg()\n", + " self.T_solver.Mult(self.z, du_dt)\n", + "\n", + " def SetParameters(self, u):\n", + " u_alpha_gf = mfem.ParGridFunction(self.fespace)\n", + " u_alpha_gf.SetFromTrueDofs(u)\n", + " for i in range(u_alpha_gf.Size()):\n", + " u_alpha_gf[i] = self.kappa + self.alpha * u_alpha_gf[i]\n", + "\n", + " self.K = mfem.ParBilinearForm(self.fespace)\n", + " u_coeff = mfem.GridFunctionCoefficient(u_alpha_gf)\n", + " self.K.AddDomainIntegrator(mfem.DiffusionIntegrator(u_coeff))\n", + " self.K.Assemble(0)\n", + " self.K.FormSystemMatrix(self.ess_tdof_list, self.Kmat)\n", + " self.T = None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6cf4052e-604a-48c1-8f11-a6e9ccfb5661", + "metadata": {}, + "outputs": [], + "source": [ + "class InitialTemperature(mfem.PyCoefficient):\n", + " def __init__(self, radius_, cx_, cy_):\n", + " self.radius = radius_\n", + " self.cx = cx_\n", + " self.cy = cy_\n", + "\n", + " mfem.PyCoefficient.__init__(self)\n", + " return\n", + " \n", + " def EvalValue(self, x):\n", + " xx = np.array(x) - np.array([self.cx, self.cy])\n", + " norm2 = np.sqrt(float(np.sum(xx**2)))\n", + " if norm2 < self.radius:\n", + " return 2.0\n", + " return 1.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5a59512-6d31-4a66-a500-a3668c1cff88", + "metadata": {}, + "outputs": [], + "source": [ + "from mpi4py import MPI\n", + "comm = MPI.COMM_WORLD\n", + "myid = comm.Get_rank()\n", + "num_procs = comm.Get_size()\n", + "\n", + "from mfem.common.arg_parser import ArgParser\n", + "parser = ArgParser(description=\"DMD - MFEM wave equation (ex23) example.\")\n", + "parser.add_argument('-m', '--mesh',\n", + " default='../data/star.mesh',\n", + " action='store', type=str,\n", + " help='Mesh file to use.')\n", + "parser.add_argument('-rs', '--refine-serial',\n", + " action='store', default=2, type=int,\n", + " help=\"Number of times to refine the mesh uniformly in serial\")\n", + "parser.add_argument('-rp', '--refine-parallel',\n", + " action='store', default=1, type=int,\n", + " help=\"Number of times to refine the mesh uniformly in parallel\")\n", + "parser.add_argument('-o', '--order',\n", + " action='store', default=2, type=int,\n", + " help=\"Finite element order (polynomial degree)\")\n", + "parser.add_argument('-s', '--ode-solver',\n", + " action='store', default=3, type=int,\n", + " help='\\n'.join([\"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3\",\n", + " \"\\t\\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.\"]))\n", + "parser.add_argument('-t', '--t-final',\n", + " action='store', default=0.5, type=float,\n", + " help=\"Final time; start time is 0.\")\n", + "parser.add_argument(\"-dt\", \"--time-step\",\n", + " action='store', default=0.01, type=float,\n", + " help=\"Time step.\")\n", + "parser.add_argument('-a', '--alpha',\n", + " action='store', default=0.01, type=float,\n", + " help='Alpha coefficient')\n", + "parser.add_argument('-k', '--kappa',\n", + " action='store', default=0.5, type=float,\n", + " help='Kappa coefficient')\n", + "parser.add_argument(\"-r\", \"--radius\",\n", + " action='store', default=0.5, type=float,\n", + " help=\"Radius of the interface of initial temperature.\")\n", + "parser.add_argument(\"-cx\", \"--center_x\",\n", + " action='store', default=0.0, type=float,\n", + " help=\"Center offset in the x direction.\")\n", + "parser.add_argument(\"-cy\", \"--center_y\",\n", + " action='store', default=0.0, type=float,\n", + " help=\"Center offset in the y direction.\")\n", + "parser.add_argument(\"-crv\", \"--crv\",\n", + " action='store', default=0.9, type=float,\n", + " help=\"DMD Closest RBF Value.\")\n", + "parser.add_argument('-vis', '--visualization',\n", + " action='store_true', default=True,\n", + " help='Enable GLVis visualization')\n", + "parser.add_argument('-no-vis', '--no-visualization',\n", + " action='store_false', dest='visualization',\n", + " help='Enable GLVis visualization')\n", + "parser.add_argument('-visit', '--visit-datafiles',\n", + " action='store_true', default=False,\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-vs\", \"--visualization-steps\",\n", + " action='store', default=5, type=int,\n", + " help=\"Visualize every n-th timestep.\")\n", + "parser.add_argument(\"-rdim\", \"--rdim\",\n", + " action='store', default=-1, type=int,\n", + " help=\"Reduced dimension for DMD.\")\n", + "parser.add_argument(\"-offline\", \"--offline\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the offline phase.\")\n", + "parser.add_argument(\"-online\", \"--online\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the online phase.\")\n", + "parser.add_argument(\"-predict\", \"--predict\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable DMD prediction.\")\n", + "parser.add_argument(\"-adios2\", \"--adios2-streams\",\n", + " action='store_true', default=False,\n", + " help=\"Save data using adios2 streams.\")\n", + "parser.add_argument(\"-save\", \"--save\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable MFEM DOF solution snapshot files).\")\n", + "parser.add_argument(\"-csv\", \"--csv\",\n", + " action='store_true', default=False, dest='csvFormat',\n", + " help=\"Use CSV or HDF format for files output by -save option.\")\n", + "parser.add_argument(\"-hdf\", \"--hdf\",\n", + " action='store_false', dest='csvFormat',\n", + " help=\"Use CSV or HDF format for files output by -save option.\")\n", + "parser.add_argument(\"-out\", \"--outputfile-name\",\n", + " action='store', default=\"\", type=str,\n", + " help=\"Name of the sub-folder to dump files within the run directory.\")\n", + "parser.add_argument(\"-pwsnap\", \"--pw-snap\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable writing pointwise snapshots.\")\n", + "parser.add_argument(\"-pwx\", \"--pwx\",\n", + " action='store', default=0, type=int,\n", + " help=\"Number of snapshot points in x\")\n", + "parser.add_argument(\"-pwy\", \"--pwy\",\n", + " action='store', default=0, type=int,\n", + " help=\"Number of snapshot points in y\")\n", + "parser.add_argument(\"-pwz\", \"--pwz\",\n", + " action='store', default=0, type=int,\n", + " help=\"Number of snapshot points in z\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1d9b7c2-72c1-4da8-b696-d1542834be48", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run for Parametric DMD:\n", + "# args = parser.parse_args(\"-r 0.1 -cx 0.1 -cy 0.1 -o 4 -visit -offline -rdim 16\".split())\n", + "# args = parser.parse_args(\"-r 0.1 -cx 0.1 -cy 0.5 -o 4 -visit -offline -rdim 16\".split())\n", + "# args = parser.parse_args(\"-r 0.1 -cx 0.5 -cy 0.5 -o 4 -visit -offline -rdim 16\".split())\n", + "# args = parser.parse_args(\"-r 0.5 -cx 0.1 -cy 0.1 -o 4 -visit -offline -rdim 16\".split())\n", + "# args = parser.parse_args(\"-r 0.25 -cx 0.2 -cy 0.4 -o 4 -online -predict\".split())\n", + "# args = parser.parse_args(\"-r 0.4 -cx 0.2 -cy 0.3 -o 4 -online -predict\".split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6533977a-df7b-478e-9e7c-347cebf0e1c9", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "if (myid == 0):\n", + " parser.print_options(args)\n", + "\n", + "precision = 8\n", + "save_dofs = args.save\n", + "basename = args.outputfile_name\n", + "offline = args.offline\n", + "online = args.online\n", + "pointwiseSnapshots = args.pw_snap\n", + "rdim = args.rdim\n", + "mesh_file = os.path.abspath(os.path.join('../data', args.mesh))\n", + "ode_solver_type = args.ode_solver\n", + "ser_ref_levels = args.refine_serial\n", + "par_ref_levels = args.refine_parallel\n", + "pwx, pwy, pwz = args.pwx, args.pwy, args.pwz\n", + "order = args.order\n", + "alpha = args.alpha\n", + "kappa = args.kappa\n", + "radius = args.radius\n", + "cx, cy = args.center_x, args.center_y\n", + "visit = args.visit_datafiles\n", + "adios2 = args.adios2_streams\n", + "visualization = args.visualization\n", + "csvFormat = args.csvFormat\n", + "dt = args.time_step\n", + "t_final = args.t_final\n", + "vis_steps = args.visualization_steps\n", + "closest_rbf_val = args.crv\n", + "predict = args.predict" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db080227-ae0b-4874-a196-9c6b88600ced", + "metadata": {}, + "outputs": [], + "source": [ + "outputPath = \".\"\n", + "if (save_dofs):\n", + " outputPath = \"run\"\n", + " if (basename != \"\"):\n", + " outputPath += \"/\" + basename\n", + "\n", + " if (myid == 0):\n", + " pathlib.Path(outputPath).mkdir(parents=True, exist_ok=True)\n", + "\n", + "check = pointwiseSnapshots ^ offline ^ online ^ save_dofs\n", + "if not check:\n", + " raise RuntimeError(\"Only one of offline, online, save, or pwsnap must be true!\")\n", + "\n", + "if (offline and (rdim <= 0)):\n", + " raise RuntimeError(\"rdim must be set.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa89c2e9-f374-4ec9-8b8a-ab2acce4aa81", + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Read the serial mesh from the given mesh file on all processors. We can\n", + "# handle triangular, quadrilateral, tetrahedral and hexahedral meshes\n", + "# with the same code.\n", + "mesh = mfem.Mesh(mesh_file, 1, 1)\n", + "dim = mesh.Dimension()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a04eb5b-fa77-4db8-9978-2164d9f71ed9", + "metadata": {}, + "outputs": [], + "source": [ + "# 4. Define the ODE solver used for time integration. Several implicit\n", + "# singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as\n", + "# explicit Runge-Kutta methods are available.\n", + "if ode_solver_type == 1:\n", + " ode_solver = mfem.BackwardEulerSolver()\n", + "elif ode_solver_type == 2:\n", + " ode_solver = mfem.SDIRK23Solver(2)\n", + "elif ode_solver_type == 3:\n", + " ode_solver = mfem.SDIRK33Solver()\n", + "elif ode_solver_type == 11:\n", + " ode_solver = mfem.ForwardEulerSolver()\n", + "elif ode_solver_type == 12:\n", + " ode_solver = mfem.RK2Solver(0.5)\n", + "elif ode_solver_type == 13:\n", + " ode_solver = mfem.RK3SSPSolver()\n", + "elif ode_solver_type == 14:\n", + " ode_solver = mfem.RK4Solver()\n", + "elif ode_solver_type == 22:\n", + " ode_solver = mfem.ImplicitMidpointSolver()\n", + "elif ode_solver_type == 23:\n", + " ode_solver = mfem.SDIRK23Solver()\n", + "elif ode_solver_type == 24:\n", + " ode_solver = mfem.SDIRK34Solver()\n", + "else:\n", + " print(\"Unknown ODE solver type: \" + str(ode_solver_type))\n", + " exit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a946ca25-10a1-44a4-bf81-172a0f2359a8", + "metadata": {}, + "outputs": [], + "source": [ + "# 5. Refine the mesh in serial to increase the resolution. In this example\n", + "# we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is\n", + "# a command-line parameter.\n", + "for lev in range(ser_ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfb75432-80c0-4016-a6f2-0c0ab2debe0e", + "metadata": {}, + "outputs": [], + "source": [ + "# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine\n", + "# this mesh further in parallel to increase the resolution. Once the\n", + "# parallel mesh is defined, the serial mesh can be deleted.\n", + "pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh)\n", + "del mesh\n", + "for x in range(par_ref_levels):\n", + " pmesh.UniformRefinement()\n", + "\n", + "# TODO(kevin): enforce user to install pymfem with gslib?\n", + "# #ifndef MFEM_USE_GSLIB\n", + "# if (pointwiseSnapshots) {\n", + "# cout << \"To use pointwise snapshots, compile with -mg option\" << endl;\n", + "# MFEM_ABORT(\"Pointwise snapshots aren't available, since the \"\n", + "# \"compilation is done without the -mg option\");\n", + "# }\n", + "pws = None\n", + "pwsnap = mfem.Vector()\n", + "pwsnap_CAROM = None\n", + "\n", + "if (pointwiseSnapshots):\n", + " pmesh.EnsureNodes()\n", + " dmdDim = [pwx, pwy, pwz]\n", + " pws = PointwiseSnapshot(dim, dmdDim)\n", + " pws.SetMesh(pmesh)\n", + "\n", + " snapshotSize = np.prod(dmdDim[:dim])\n", + "\n", + " pwsnap.SetSize(snapshotSize)\n", + " if (myid == 0):\n", + " pwsnap_CAROM = linalg.Vector(pwsnap.GetDataArray(), True, False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37df6a1d-c0e2-4df7-95a1-b24d6a2a0ae6", + "metadata": {}, + "outputs": [], + "source": [ + "# 7. Define the vector finite element space representing the current and the\n", + "# initial temperature, u_ref.\n", + "fe_coll = mfem.H1_FECollection(order, dim)\n", + "fespace = mfem.ParFiniteElementSpace(pmesh, fe_coll)\n", + "\n", + "fe_size = fespace.GlobalTrueVSize()\n", + "if (myid == 0):\n", + " print(\"Number of temperature unknowns: %d\" % fe_size)\n", + "\n", + "u_gf = mfem.ParGridFunction(fespace)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7ca8f6a-e6a3-45ef-8983-0cac8cfe851e", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Set the initial conditions for u. All boundaries are considered\n", + "# natural.\n", + "u_0 = InitialTemperature(radius, cx, cy)\n", + "u_gf.ProjectCoefficient(u_0)\n", + "u = mfem.Vector()\n", + "u_gf.GetTrueDofs(u)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ae64f02-d078-4a47-8e08-248ad9909520", + "metadata": {}, + "outputs": [], + "source": [ + "# 9. Initialize the conduction operator and the VisIt visualization.\n", + "oper = ConductionOperator(fespace, alpha, kappa, u)\n", + "u_gf.SetFromTrueDofs(u)\n", + "\n", + "mesh_name = \"%s/parametric_heat_conduction_%f_%f_%f_%f-mesh.%06d\" % (outputPath, radius, alpha, cx, cy, myid)\n", + "sol_name = \"%s/parametric_heat_conduction_%f_%f_%f_%f-init.%06d\" % (outputPath, radius, alpha, cx, cy, myid)\n", + "\n", + "pmesh.Print(mesh_name, precision)\n", + "\n", + "output = io.StringIO()\n", + "output.precision = precision\n", + "u_gf.Save(output)\n", + "fid = open(sol_name, 'w')\n", + "fid.write(output.getvalue())\n", + "fid.close()\n", + "\n", + "visit_name = \"%s/Parametric_Heat_Conduction_%f_%f_%f_%f\" % (outputPath, radius, alpha, cx, cy)\n", + "visit_dc = mfem.VisItDataCollection(visit_name, pmesh)\n", + "visit_dc.RegisterField(\"temperature\", u_gf)\n", + "if (visit):\n", + " visit_dc.SetCycle(0)\n", + " visit_dc.SetTime(0.0)\n", + " visit_dc.Save()\n", + "\n", + "# Optionally output a BP (binary pack) file using ADIOS2. This can be\n", + "# visualized with the ParaView VTX reader.\n", + "# TODO(kevin): enforce user to install pymfem with adios2?\n", + "#ifdef MFEM_USE_ADIOS2\n", + "if (adios2):\n", + " postfix = mesh_file[len(\"../data/\"):]\n", + " postfix += \"_o%d\" % order\n", + " postfix += \"_solver%d\" % ode_solver_type\n", + " collection_name = \"%s/parametric_heat_conduction-p-%s.bp\" % (outputPath, postfix)\n", + "\n", + " adios2_dc = mfem.ADIOS2DataCollection(MPI.COMM_WORLD, collection_name, pmesh)\n", + " adios2_dc.SetParameter(\"SubStreams\", \"%d\" % (num_procs/2) )\n", + " adios2_dc.RegisterField(\"temperature\", u_gf)\n", + " adios2_dc.SetCycle(0)\n", + " adios2_dc.SetTime(0.0)\n", + " adios2_dc.Save()\n", + "#endif" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91e180dd-ae83-4d48-9366-cdaadace4320", + "metadata": {}, + "outputs": [], + "source": [ + "if visualization:\n", + " sol_sock = mfem.socketstream(\"localhost\", 19916)\n", + " if not sol_sock.good():\n", + " print(\"Unable to connect to GLVis server at localhost:19916\")\n", + " visualization = False\n", + " print(\"GLVis visualization disabled.\")\n", + " else:\n", + " sol_sock << \"parallel \" << num_procs << \" \" << myid << \"\\n\"\n", + " sol_sock.precision(precision)\n", + " sol_sock << \"solution\\n\" << pmesh << u_gf\n", + " sol_sock << \"pause\\n\"\n", + " sol_sock.flush()\n", + " print(\n", + " \"GLVis visualization paused. Press space (in the GLVis window) to resume it.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22a9c517-0749-4a4a-a189-9cdc8a99cd77", + "metadata": {}, + "outputs": [], + "source": [ + "#ifdef MFEM_USE_GSLIB\n", + "# TODO(kevin): enforce user to install pymfem with gslib?\n", + "if (pointwiseSnapshots):\n", + " pws.GetSnapshot(u_gf, pwsnap)\n", + "\n", + " dmd_filename = \"snap_%f_%f_%f_%f_0\" % (radius, alpha, cx, cy)\n", + " if (myid == 0):\n", + " print(\"Writing DMD snapshot at step 0, time 0.0\")\n", + " pwsnap_CAROM.write(dmd_filename)\n", + "#endif\n", + "\n", + "fom_timer, dmd_training_timer, dmd_prediction_timer = StopWatch(), StopWatch(), StopWatch()\n", + "\n", + "fom_timer.Start()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41c1aaeb-8866-4bf3-b2d7-ba1808887f76", + "metadata": {}, + "outputs": [], + "source": [ + "# 10. Perform time-integration (looping over the time iterations, ti, with a\n", + "# time-step dt).\n", + "ode_solver.Init(oper)\n", + "t = 0.0\n", + "ts = []\n", + "# CAROM::Vector* init = NULL;\n", + "\n", + "# CAROM::Database *db = NULL;\n", + "if (csvFormat):\n", + " db = utils.CSVDatabase()\n", + "else:\n", + " db = utils.HDFDatabase()\n", + "\n", + "snap_list = []\n", + "\n", + "fom_timer.Stop()\n", + "\n", + "# CAROM::DMD* dmd_u = NULL;\n", + "\n", + "if (offline):\n", + " dmd_training_timer.Start()\n", + "\n", + " # 11. Create DMD object and take initial sample.\n", + " u_gf.SetFromTrueDofs(u)\n", + " dmd_u = algo.DMD(u.Size(), dt)\n", + " dmd_u.takeSample(u.GetDataArray(), t)\n", + "\n", + " if (myid == 0):\n", + " print(\"Taking snapshot at: %f\" % t)\n", + "\n", + " dmd_training_timer.Stop()\n", + "\n", + "if (online):\n", + " u_gf.SetFromTrueDofs(u)\n", + " init = linalg.Vector(u.GetDataArray(), u.Size(), True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d702dc3-96d4-4c03-89c8-dad75233d906", + "metadata": {}, + "outputs": [], + "source": [ + "if (save_dofs and (myid == 0)):\n", + " if (csvFormat):\n", + " pathlib.Path(\"%s/step0\" % outputPath).mkdir(parents=True, exist_ok=True)\n", + " db.putDoubleArray(outputPath + \"/step0/sol.csv\", u.GetDataArray(), u.Size())\n", + " else:\n", + " db.create(outputPath + \"/dmd_0.hdf\")\n", + " db.putDoubleArray(\"step0sol\", u.GetDataArray(), u.Size())\n", + "\n", + "ts += [t]\n", + "snap_list += [0]\n", + "\n", + "last_step = False\n", + "ti = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "feeb685c-a84f-413b-af6e-62020337fd9a", + "metadata": {}, + "outputs": [], + "source": [ + "while (not last_step):\n", + " fom_timer.Start()\n", + "\n", + " if (t + dt >= t_final - dt / 2.):\n", + " last_step = True\n", + "\n", + " t, dt = ode_solver.Step(u, t, dt)\n", + "\n", + " fom_timer.Stop()\n", + "\n", + " if (offline):\n", + " dmd_training_timer.Start()\n", + "\n", + " u_gf.SetFromTrueDofs(u)\n", + " dmd_u.takeSample(u.GetDataArray(), t)\n", + "\n", + " if (myid == 0):\n", + " print(\"Taking snapshot at: %f\" % t)\n", + "\n", + " dmd_training_timer.Stop()\n", + "\n", + " if (save_dofs and (myid == 0)):\n", + " if (csvFormat):\n", + " pathlib.Path(\"%s/step%d\" % (outputPath, ti)).mkdir(parents=True, exist_ok=True)\n", + " db.putDoubleArray(\"%s/step%d/sol.csv\" % (outputPath, ti), u.GetDataArray(), u.Size())\n", + " else:\n", + " db.putDoubleArray(\"step%dsol\" % ti, u.GetDataArray(), u.Size())\n", + "\n", + " ts += [t]\n", + " snap_list += [ti]\n", + "\n", + " if (last_step or ((ti % vis_steps) == 0)):\n", + " if (myid == 0):\n", + " print(\"step %d, t = %f\" % (ti, t))\n", + "\n", + " u_gf.SetFromTrueDofs(u)\n", + " if (visualization):\n", + " if sol_sock.good():\n", + " sol_sock << \"parallel \" << num_procs << \" \" << myid << \"\\n\"\n", + " sol_sock << \"solution\\n\" << pmesh << u_gf\n", + " # sol_sock << \"pause\\n\"\n", + " sol_sock.flush()\n", + "\n", + " if (visit):\n", + " visit_dc.SetCycle(ti)\n", + " visit_dc.SetTime(t)\n", + " visit_dc.Save()\n", + "\n", + "#ifdef MFEM_USE_ADIOS2\n", + " if (adios2):\n", + " adios2_dc.SetCycle(ti)\n", + " adios2_dc.SetTime(t)\n", + " adios2_dc.Save()\n", + "#endif\n", + "\n", + "#ifdef MFEM_USE_GSLIB\n", + " if (pointwiseSnapshots):\n", + " pws.GetSnapshot(u_gf, pwsnap)\n", + "\n", + " dmd_filename = \"snap_%f_%f_%f_%f_%d\" % (radius, alpha, cx, cy, ti)\n", + " if (myid == 0):\n", + " print(\"Writing DMD snapshot at step %d, time %f\" % (ti, t))\n", + " pwsnap_CAROM.write(dmd_filename)\n", + "#endif\n", + "\n", + " oper.SetParameters(u)\n", + "\n", + " ti += 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "316a5fad-b8ce-4d50-8088-b313f6870d6f", + "metadata": {}, + "outputs": [], + "source": [ + "if (save_dofs and (myid == 0)):\n", + " if (csvFormat):\n", + " db.putDoubleVector(outputPath + \"/tval.csv\", ts, len(ts))\n", + " db.putInteger(outputPath + \"/numsnap\", len(snap_list))\n", + " db.putIntegerArray(outputPath + \"/snap_list.csv\", snap_list, len(snap_list))\n", + " else:\n", + " db.putDoubleVector(\"tval\", ts, len(ts))\n", + " db.putInteger(\"numsnap\", len(snap_list))\n", + " db.putInteger(\"snap_bound_size\", 0)\n", + " db.putIntegerArray(\"snap_list\", snap_list, len(snap_list))\n", + "#ifdef MFEM_USE_ADIOS2\n", + "if (adios2):\n", + " del adios2_dc\n", + "#endif" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fecb068-fb2d-4c4b-9a1c-f73ed1a326bf", + "metadata": {}, + "outputs": [], + "source": [ + "# 12. Save the final solution in parallel. This output can be viewed later\n", + "# using GLVis: \"glvis -np -m parametric_heat_conduction-mesh -g parametric_heat_conduction-final\".\n", + "sol_name = \"%s/parametric_heat_conduction_%f_%f_%f_%f-final.%06d\" % (outputPath, radius, alpha, cx, cy, myid)\n", + "output = io.StringIO()\n", + "output.precision = precision\n", + "u_gf.Save(output)\n", + "fid = open(sol_name, 'w')\n", + "fid.write(output.getvalue())\n", + "fid.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8abf542b-e6ef-4637-9f7b-410ee00038b7", + "metadata": {}, + "outputs": [], + "source": [ + "# 13. Calculate the DMD modes.\n", + "if (offline or online):\n", + " if (offline):\n", + " if (myid == 0):\n", + " print(\"Creating DMD with rdim: %d\" % rdim)\n", + "\n", + " dmd_training_timer.Start()\n", + "\n", + " dmd_u.train(rdim)\n", + "\n", + " dmd_training_timer.Stop()\n", + "\n", + " dmd_u.save(\"%s/%f_%f_%f_%f\" % (outputPath, radius, alpha, cx, cy))\n", + "\n", + " if (myid == 0):\n", + " with open(\"parameters.txt\", \"ab\") as f:\n", + " np.savetxt(f, [radius, alpha, cx, cy])\n", + "\n", + " if (online):\n", + " if (myid == 0):\n", + " print(\"Creating DMD using the rdim of the offline phase\")\n", + "\n", + " param_hist = np.loadtxt(\"parameters.txt\")\n", + " param_hist = param_hist.reshape((int(param_hist.size / 4), 4))\n", + "\n", + " dmd_paths = []\n", + " param_vectors = []\n", + "\n", + " for curr_param in param_hist:\n", + " curr_radius, curr_alpha, curr_cx, curr_cy = curr_param\n", + "\n", + " dmd_paths += [\"%s/%f_%f_%f_%f\" % (outputPath, curr_radius, curr_alpha, curr_cx, curr_cy)]\n", + " param_vectors += [linalg.Vector([curr_radius, curr_alpha, curr_cx, curr_cy], False)]\n", + "\n", + " desired_param = linalg.Vector([radius, alpha, cx, cy], False)\n", + "\n", + " dmd_training_timer.Start()\n", + "\n", + " dmd_u = algo.getParametricDMD(algo.DMD, param_vectors, dmd_paths, desired_param,\n", + " \"G\", \"LS\", closest_rbf_val)\n", + "\n", + " dmd_u.projectInitialCondition(init)\n", + "\n", + " dmd_training_timer.Stop()\n", + " del desired_param\n", + "\n", + " if (predict):\n", + " true_solution_u = mfem.Vector(u.GetDataArray(), u.Size())\n", + "\n", + " dmd_prediction_timer.Start()\n", + "\n", + " # 14. Predict the state at t_final using DMD.\n", + " if (myid == 0):\n", + " print(\"Predicting temperature using DMD at: %f\" % ts[0])\n", + "\n", + " result_u = dmd_u.predict(ts[0])\n", + " initial_dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + " u_gf.SetFromTrueDofs(initial_dmd_solution_u)\n", + "\n", + " visit_name = \"%s/DMD_Parametric_Heat_Conduction_%f_%f_%f_%f\" % (outputPath, radius, alpha, cx, cy)\n", + " dmd_visit_dc = mfem.VisItDataCollection(visit_name, pmesh)\n", + " dmd_visit_dc.RegisterField(\"temperature\", u_gf)\n", + " if (visit):\n", + " dmd_visit_dc.SetCycle(0)\n", + " dmd_visit_dc.SetTime(0.0)\n", + " dmd_visit_dc.Save()\n", + "\n", + " del result_u\n", + "\n", + " if (visit):\n", + " for i, tsi in enumerate(ts):\n", + " if ((i == len(ts) - 1) or ((i % vis_steps) == 0)):\n", + " result_u = dmd_u.predict(tsi)\n", + " if (myid == 0):\n", + " print(\"Predicting temperature using DMD at: %f\" % tsi)\n", + "\n", + " dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + " u_gf.SetFromTrueDofs(dmd_solution_u)\n", + "\n", + " dmd_visit_dc.SetCycle(i)\n", + " dmd_visit_dc.SetTime(tsi)\n", + " dmd_visit_dc.Save()\n", + "\n", + " del result_u\n", + "\n", + " dmd_prediction_timer.Stop()\n", + "\n", + " result_u = dmd_u.predict(t_final)\n", + "\n", + " # 15. Calculate the relative error between the DMD final solution and the true solution.\n", + " dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + " diff_u = mfem.Vector(true_solution_u.Size())\n", + " mfem.subtract_vector(dmd_solution_u, true_solution_u, diff_u)\n", + "\n", + " tot_diff_norm_u = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff_u, diff_u))\n", + " tot_true_solution_u_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD,\n", + " true_solution_u, true_solution_u))\n", + "\n", + " if (myid == 0):\n", + " print(\"Relative error of DMD temperature (u) at t_final: %f is %f\" % (t_final,\n", + " tot_diff_norm_u / tot_true_solution_u_norm))\n", + "\n", + " print(\"Elapsed time for predicting DMD: %e second\" % dmd_prediction_timer.duration)\n", + "\n", + " del result_u\n", + "\n", + " if (myid == 0):\n", + " print(\"Elapsed time for training DMD: %e second\" % dmd_training_timer.duration)\n", + "\n", + "if (myid == 0):\n", + " print(\"Elapsed time for solving FOM: %e second\" % fom_timer.duration)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "874121dd-46eb-46d9-b33f-006272fc0ade", + "metadata": {}, + "outputs": [], + "source": [ + "# 16. Free the used memory.\n", + "del ode_solver\n", + "del pmesh\n", + "if (offline):\n", + " del dmd_u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76d7f673-2d84-453b-93e7-8117266c29d8", + "metadata": {}, + "outputs": [], + "source": [ + "#ifdef MFEM_USE_GSLIB\n", + "del pws\n", + "del pwsnap_CAROM\n", + "#endif\n", + "\n", + "MPI.Finalize()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/dmd/wave_equation.ipynb b/examples/dmd/wave_equation.ipynb new file mode 100644 index 0000000..00c3baf --- /dev/null +++ b/examples/dmd/wave_equation.ipynb @@ -0,0 +1,633 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "98f8cede-1466-4742-b9a4-dc742fb8890b", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import io\n", + "import pathlib\n", + "import sys\n", + "import time\n", + "try:\n", + " import mfem.par as mfem\n", + "except ModuleNotFoundError:\n", + " msg = \"PyMFEM is not installed yet. Install PyMFEM:\\n\"\n", + " msg += \"\\tgit clone https://github.com/mfem/PyMFEM.git\\n\"\n", + " msg += \"\\tcd PyMFEM\\n\"\n", + " msg += \"\\tpython3 setup.py install --with-parallel\\n\"\n", + " raise ModuleNotFoundError(msg)\n", + " \n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum\n", + "sys.path.append(\"../../build\")\n", + "import pylibROM.algo as algo\n", + "import pylibROM.linalg as linalg" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e1eba099-9048-4b39-a30b-2f3448c8e498", + "metadata": {}, + "outputs": [], + "source": [ + "# pyMFEM does not provide mfem::StopWatch.\n", + "# there is not really a similar stopwatch package in python.. (seriously?)\n", + "class StopWatch:\n", + " import time\n", + " duration = 0.0\n", + " start_time = 0.0\n", + " stop_time = 0.0\n", + " running = False\n", + "\n", + " def __init__(self):\n", + " self.Reset()\n", + " return\n", + " \n", + " def Start(self):\n", + " assert(not self.running)\n", + " self.start_time = time.time()\n", + " self.running = True\n", + " return\n", + " \n", + " def Stop(self):\n", + " assert(self.running)\n", + " self.stop_time = time.time()\n", + " self.duration += self.stop_time - self.start_time\n", + " self.running = False\n", + " return\n", + " \n", + " def Reset(self):\n", + " self.duration = 0.0\n", + " self.start_time = 0.0\n", + " self.stop_time = 0.0\n", + " self.running = False\n", + " return\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81554a20-12e5-4616-904e-ca2b8a568633", + "metadata": {}, + "outputs": [], + "source": [ + "class WaveOperator(mfem.SecondOrderTimeDependentOperator):\n", + " def __init__(self, fespace, ess_bdr, speed):\n", + " mfem.SecondOrderTimeDependentOperator.__init__(\n", + " self, fespace.GetTrueVSize(), 0.0)\n", + "\n", + " self.ess_tdof_list = mfem.intArray()\n", + "\n", + " rel_tol = 1e-8\n", + " fespace.GetEssentialTrueDofs(ess_bdr, self.ess_tdof_list)\n", + "\n", + " c2 = mfem.ConstantCoefficient(speed*speed)\n", + " K = mfem.BilinearForm(fespace)\n", + " K.AddDomainIntegrator(mfem.DiffusionIntegrator(c2))\n", + " K.Assemble()\n", + "\n", + " self.Kmat0 = mfem.SparseMatrix()\n", + " self.Kmat = mfem.SparseMatrix()\n", + " dummy = mfem.intArray()\n", + " K.FormSystemMatrix(dummy, self.Kmat0)\n", + " K.FormSystemMatrix(self.ess_tdof_list, self.Kmat)\n", + " self.K = K\n", + "\n", + " self.Mmat = mfem.SparseMatrix()\n", + " M = mfem.BilinearForm(fespace)\n", + " M.AddDomainIntegrator(mfem.MassIntegrator())\n", + " M.Assemble()\n", + " M.FormSystemMatrix(self.ess_tdof_list, self.Mmat)\n", + " self.M = M\n", + "\n", + " M_solver = mfem.CGSolver()\n", + " M_prec = mfem.DSmoother()\n", + " M_solver.iterative_mode = False\n", + " M_solver.SetRelTol(rel_tol)\n", + " M_solver.SetAbsTol(0.0)\n", + " M_solver.SetMaxIter(30)\n", + " M_solver.SetPrintLevel(0)\n", + " M_solver.SetPreconditioner(M_prec)\n", + " M_solver.SetOperator(self.Mmat)\n", + " self.M_prec = M_prec\n", + " self.M_solver = M_solver\n", + "\n", + " T_solver = mfem.CGSolver()\n", + " T_prec = mfem.DSmoother()\n", + " T_solver.iterative_mode = False\n", + " T_solver.SetRelTol(rel_tol)\n", + " T_solver.SetAbsTol(0.0)\n", + " T_solver.SetMaxIter(100)\n", + " T_solver.SetPrintLevel(0)\n", + " T_solver.SetPreconditioner(T_prec)\n", + " self.T_prec = T_prec\n", + " self.T_solver = T_solver\n", + " self.T = None\n", + "\n", + " def Mult(self, u, du_dt, d2udt2):\n", + " # Compute:\n", + " # d2udt2 = M^{-1}*-K(u)\n", + " # for d2udt2\n", + " z = mfem.Vector(u.Size())\n", + " self.Kmat.Mult(u, z)\n", + " z.Neg() # z = -z\n", + " self.M_solver.Mult(z, d2udt2)\n", + "\n", + " def ImplicitSolve(self, fac0, fac1, u, dudt, d2udt2):\n", + " # Solve the equation:\n", + " # d2udt2 = M^{-1}*[-K(u + fac0*d2udt2)]\n", + " # for d2udt2\n", + " if self.T is None:\n", + " self.T = mfem.Add(1.0, self.Mmat, fac0, self.Kmat)\n", + " self.T_solver.SetOperator(self.T)\n", + " z = mfem.Vector(u.Size())\n", + " self.Kmat0.Mult(u, z)\n", + " z.Neg()\n", + "\n", + " # iterate over Array :D\n", + " for j in self.ess_tdof_list:\n", + " z[j] = 0.0\n", + "\n", + " self.T_solver.Mult(z, d2udt2)\n", + "\n", + " def SetParameters(self, u):\n", + " self.T = None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29781759-957e-43c3-9975-31a6e5ab2c28", + "metadata": {}, + "outputs": [], + "source": [ + "class InitialSolution(mfem.PyCoefficient):\n", + " def EvalValue(self, x):\n", + " norm2 = sum(x**2)\n", + " return exp(-norm2*30)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60e4d06f-84ff-4de1-8231-613e7c3645a1", + "metadata": {}, + "outputs": [], + "source": [ + "class InitialRate(mfem.PyCoefficient):\n", + " def EvalValue(self, x):\n", + " return 0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1538376-b14a-41db-b69d-6542e2c5124a", + "metadata": {}, + "outputs": [], + "source": [ + "from mfem.common.arg_parser import ArgParser\n", + "parser = ArgParser(description=\"DMD - MFEM wave equation (ex23) example.\")\n", + "parser.add_argument('-m', '--mesh',\n", + " default='../data/star.mesh',\n", + " action='store', type=str,\n", + " help='Mesh file to use.')\n", + "parser.add_argument('-r', '--refine',\n", + " action='store', default=2, type=int,\n", + " help=\"Number of times to refine the mesh uniformly before parallel\")\n", + "parser.add_argument('-o', '--order',\n", + " action='store', default=2, type=int,\n", + " help=\"Finite element order (polynomial degree)\")\n", + "help_ode = '\\n'.join([\"ODE solver: [0--10] \\t- GeneralizedAlpha(0.1 * s),\",\n", + " \"11 \\t - Average Acceleration,\",\n", + " \"12 \\t - Linear Acceleration\",\n", + " \"13 \\t- CentralDifference\",\n", + " \"14 \\t- FoxGoodwin\"])\n", + "parser.add_argument('-s', '--ode-solver',\n", + " action='store', default=10, type=int,\n", + " help=help_ode)\n", + "parser.add_argument('-tf', '--t-final',\n", + " action='store', default=0.5, type=float,\n", + " help=\"Final time; start time is 0.\")\n", + "parser.add_argument('-dt', '--time-step',\n", + " action='store', default=1e-2, type=float,\n", + " help=\"Time step\")\n", + "parser.add_argument(\"-c\", \"--speed\",\n", + " action='store', default=1.0, type=float,\n", + " help=\"Wave speed.\")\n", + "parser.add_argument(\"-neu\", \"--neumann\",\n", + " action='store_true', default=False,\n", + " help=\"BC switch.\")\n", + "parser.add_argument('-vis', '--visualization',\n", + " action='store_true', default=False,\n", + " help='Enable GLVis visualization')\n", + "parser.add_argument('-visit', '--visit-datafiles',\n", + " action='store_true', default=False,\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-vs\", \"--visualization-steps\",\n", + " action='store', default=5, type=int,\n", + " help=\"Visualize every n-th timestep.\")\n", + "\n", + "parser.add_argument('-rd', '--reference',\n", + " default='', action='store', type=str,\n", + " help='Reference directory for checking final solution.')\n", + "parser.add_argument('-ef', '--energy_fraction',\n", + " action='store', default=0.9999, type=float,\n", + " help=\"Energy fraction for DMD.\")\n", + "parser.add_argument('-rdim', '--rdim',\n", + " action='store', default=-1, type=int,\n", + " help=\"Reduced dimension for DMD.\")\n", + "parser.add_argument('-nwinsamp', '--numwindowsamples',\n", + " action='store', default=sys.maxsize, type=int,\n", + " help=\"Number of samples in DMD windows.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02291ca2-7c87-4e2a-8d74-62aeaaaae31d", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run for Time-Windowing DMD:\n", + "# args = parser.parse_args(\"-o 4 -tf 5 -nwinsamp 25\".split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ba14cd8-5309-4178-97cb-21daa3b17107", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "parser.print_options(args)\n", + "mesh_file = expanduser(os.path.join('../data', args.mesh))\n", + "ref_levels = args.refine\n", + "order = args.order\n", + "ode_solver_type = args.ode_solver\n", + "t_final = args.t_final\n", + "dt = args.time_step\n", + "speed = args.speed\n", + "dirichlet = (not args.neumann)\n", + "visit = args.visit_datafiles\n", + "visualization = args.visualization\n", + "vis_steps = args.visualization_steps\n", + "ref_dir = args.reference\n", + "ef = args.energy_fraction\n", + "rdim = args.rdim\n", + "windowNumSamples = args.numwindowsamples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4586ee22-f647-436a-ba4d-176dbeec4cea", + "metadata": {}, + "outputs": [], + "source": [ + "if ((rdim <= 0) and (rdim != -1)):\n", + " raise ValueError(\"rdim is set to %d, rdim can only be a positive integer or -1\" % rdim)\n", + "\n", + "if (ef <= 0.0):\n", + " raise ValueError(\"ef must be a positive, it is %f\" % ef)\n", + "elif (rdim != -1):\n", + " print(\"rdim is set to %d\" % rdim)\n", + "\n", + "mesh = mfem.Mesh(mesh_file, 1, 1)\n", + "dim = mesh.Dimension()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7840e315-ff4c-4b1e-99a3-59c7909598c7", + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Define the ODE solver used for time integration. Several second order\n", + "# time integrators are available.\n", + "if ode_solver_type <= 10:\n", + " ode_solver = mfem.GeneralizedAlpha2Solver(ode_solver_type / 10.)\n", + "elif ode_solver_type == 11:\n", + " ode_solver = mfem.AverageAccelerationSolver()\n", + "elif ode_solver_type == 12:\n", + " ode_solver = mfem.LinearAccelerationSolver()\n", + "elif ode_solver_type == 13:\n", + " ode_solver = mfem.CentralDifferenceSolver()\n", + "elif ode_solver_type == 14:\n", + " ode_solver = mfem.FoxGoodwinSolver()\n", + "else:\n", + " print(\"Unknown ODE solver type: \" + str(ode_solver_type))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e676b3f5-5b40-46c4-bd32-25cc680cd27d", + "metadata": {}, + "outputs": [], + "source": [ + "# 4. Refine the mesh to increase the resolution. In this example we do\n", + "# 'ref_levels' of uniform refinement, where 'ref_levels' is a\n", + "# command-line parameter.\n", + "for lev in range(ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a14fcd51-7b98-46ef-a2b4-a532f9de3153", + "metadata": {}, + "outputs": [], + "source": [ + "# 5. Define the vector finite element space representing the current and the\n", + "# initial temperature, u_ref.\n", + "fe_coll = mfem.H1_FECollection(order, dim)\n", + "fespace = mfem.FiniteElementSpace(mesh, fe_coll)\n", + "\n", + "fe_size = fespace.GetTrueVSize()\n", + "print(\"Number of temperature unknowns: \" + str(fe_size))\n", + "\n", + "u_gf = mfem.GridFunction(fespace)\n", + "dudt_gf = mfem.GridFunction(fespace)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d199926-043d-45a1-9be1-df332cd9f8e0", + "metadata": {}, + "outputs": [], + "source": [ + "# 6. Set the initial conditions for u. All boundaries are considered\n", + "# natural.\n", + "u_0 = InitialSolution()\n", + "dudt_0 = InitialRate()\n", + "\n", + "u_gf.ProjectCoefficient(u_0)\n", + "u = mfem.Vector()\n", + "u_gf.GetTrueDofs(u)\n", + "\n", + "dudt_gf.ProjectCoefficient(dudt_0)\n", + "dudt = mfem.Vector()\n", + "dudt_gf.GetTrueDofs(dudt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68b27976-65d4-4c02-952d-4ff7076bfb92", + "metadata": {}, + "outputs": [], + "source": [ + "# 7. Initialize the conduction operator and the visualization.\n", + "ess_bdr = mfem.intArray()\n", + "if mesh.bdr_attributes.Size():\n", + " ess_bdr.SetSize(mesh.bdr_attributes.Max())\n", + " if (dirichlet):\n", + " ess_bdr.Assign(1)\n", + " else:\n", + " ess_bdr.Assigne(0)\n", + "\n", + "oper = WaveOperator(fespace, ess_bdr, speed)\n", + "\n", + "u_gf.SetFromTrueDofs(u)\n", + "\n", + "mesh.Print(\"wave_equation.mesh\", 8)\n", + "output = io.StringIO()\n", + "output.precision = 8\n", + "u_gf.Save(output)\n", + "dudt_gf.Save(output)\n", + "fid = open(\"wave_equation-init.gf\", 'w')\n", + "fid.write(output.getvalue())\n", + "fid.close()\n", + "\n", + "if visit:\n", + " visit_dc = mfem.VisItDataCollection(\"Wave_Equation\", mesh)\n", + " visit_dc.RegisterField(\"solution\", u_gf)\n", + " visit_dc.RegisterField(\"rate\", dudt_gf)\n", + " visit_dc.SetCycle(0)\n", + " visit_dc.SetTime(0.0)\n", + " visit_dc.Save()\n", + "\n", + "if visualization:\n", + " sout = mfem.socketstream(\"localhost\", 19916)\n", + " if not sout.good():\n", + " print(\"Unable to connect to GLVis server at localhost:19916\")\n", + " visualization = False\n", + " print(\"GLVis visualization disabled.\")\n", + " else:\n", + " sout.precision(precision)\n", + " sout << \"solution\\n\" << mesh << dudt_gf\n", + " sout << \"pause\\n\"\n", + " sout.flush()\n", + " print(\n", + " \"GLVis visualization paused. Press space (in the GLVis window) to resume it.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "342f630c-4684-404a-9e5b-1459a10c77f2", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Perform time-integration (looping over the time iterations, ti, with a\n", + "# time-step dt).\n", + "# mfem::StopWatch is not binded by pyMFEM.\n", + "fom_timer, dmd_training_timer, dmd_prediction_timer = StopWatch(), StopWatch(), StopWatch()\n", + "fom_timer.Start()\n", + "ode_solver.Init(oper)\n", + "t = 0.0\n", + "fom_timer.Stop()\n", + "dmd_training_timer.Start()\n", + "curr_window = 0\n", + "ts, dmd_u = [], []\n", + "dmd_u += [algo.DMD(u.Size(), dt)]\n", + "\n", + "# NOTE: mfem Vector::GetData returns a SWIG Object of type double *.\n", + "# To make it compatible with pybind11, we use ctypes to read data from the memory address.\n", + "from ctypes import *\n", + "uData = (c_double * u.Size()).from_address(int(u.GetData())) # this does not copy the data.\n", + "# uData = list(uData) # this copies the data.\n", + "uData = np.array(uData, copy=False)\n", + "\n", + "# Showing the memory address info\n", + "print(\"All of these memory addresses are different.\")\n", + "print(\"id(uData[0]): %d\" % id(uData[0]))\n", + "print(\"int(u.GetData()): %d\" % (int(u.GetData()))) # this is not the same as u[0], yet still points to the data.\n", + "print(\"id(uData): %d\" % id(uData)) # this is not the same as u[0], yet still points to the data.\n", + "\n", + "print(\"But uData[*] points to the right memory.\")\n", + "print(\"id(u[0]): %d =? id(uData[0]): %d\" % (id(u[0]), id(uData[0])))\n", + "print(\"id(u[1]): %d =? id(uData[1]): %d\" % (id(u[1]), id(uData[1])))\n", + "print(\"uData type: %s\" % type(uData))\n", + "\n", + "dmd_u[curr_window].takeSample(uData, t)\n", + "ts += [t]\n", + "dmd_training_timer.Stop()\n", + "\n", + "last_step = False\n", + "ti = 0\n", + "while not last_step:\n", + " ti += 1\n", + " if t + dt >= t_final - dt/2:\n", + " last_step = True\n", + "\n", + " fom_timer.Start()\n", + " t, dt = ode_solver.Step(u, dudt, t, dt)\n", + " fom_timer.Stop()\n", + "\n", + " dmd_training_timer.Start()\n", + " dmd_u[curr_window].takeSample(uData, t)\n", + " \n", + " if (last_step or (ti % windowNumSamples == 0)):\n", + " print(\"step %d, t= %f\" % (ti, t))\n", + "\n", + " if (rdim != -1):\n", + " print(\"Creating DMD with rdim %d at window index: %d\" % (rdim, curr_window))\n", + " dmd_u[curr_window].train(rdim)\n", + " else:\n", + " print(\"Creating DMD with energy fraction: %f at window index: %d\" % (ef, curr_window))\n", + " dmd_u[curr_window].train(ef)\n", + "\n", + " if (not last_step):\n", + " curr_window += 1\n", + " dmd_u += [algo.DMD(u.Size(), dt)]\n", + " dmd_u[curr_window].takeSample(uData, t)\n", + " ts += [t]\n", + " dmd_training_timer.Stop()\n", + "\n", + " if last_step or (ti % vis_steps == 0):\n", + " print(\"step \" + str(ti) + \", t = \" + \"{:g}\".format(t))\n", + "\n", + " u_gf.SetFromTrueDofs(u)\n", + " dudt_gf.SetFromTrueDofs(dudt)\n", + " if visualization:\n", + " sout << \"solution\\n\" << mesh << u_gf\n", + " sout.flush()\n", + "\n", + " if visit:\n", + " visit_dc.SetCycle(ti)\n", + " visit_dc.SetTime(t)\n", + " visit_dc.Save()\n", + "\n", + " oper.SetParameters(u)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb68e215-9541-49bf-8040-4dba942e9157", + "metadata": {}, + "outputs": [], + "source": [ + "# 9. Save the final solution. This output can be viewed later using GLVis:\n", + "# \"glvis -m wave_equation.mesh -g wave_equation-final.gf\".\n", + "output = io.StringIO()\n", + "output.precision = 8\n", + "u_gf.Save(output)\n", + "dudt_gf.Save(output)\n", + "fid = open(\"wave_equation-final.gf\", 'w')\n", + "fid.write(output.getvalue())\n", + "fid.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1cd37a0-c6fa-4b86-b87a-f14bc094e649", + "metadata": {}, + "outputs": [], + "source": [ + "# 10. Predict the state at t_final using DMD.\n", + "print(\"Predicting temperature using DMD\")\n", + "dmd_visit_dc = mfem.VisItDataCollection(\"DMD_Wave_Equation\", mesh)\n", + "dmd_visit_dc.RegisterField(\"solution\", u_gf)\n", + "curr_window = 0\n", + "if (visit):\n", + " dmd_prediction_timer.Start()\n", + " result_u = dmd_u[curr_window].predict(ts[0])\n", + " dmd_prediction_timer.Stop()\n", + "\n", + " # result_u.getData() returns a numpy array, which shares the memory buffer.\n", + " # result_u.getData() does not own the memory.\n", + " initial_dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + " u_gf.SetFromTrueDofs(initial_dmd_solution_u)\n", + " dmd_visit_dc.SetCycle(0)\n", + " dmd_visit_dc.SetTime(0.0)\n", + " dmd_visit_dc.Save()\n", + "\n", + "for i in range(1, len(ts)):\n", + " if ((i == len(ts) - 1) or (i % vis_steps == 0)):\n", + " if (visit):\n", + " dmd_prediction_timer.Start()\n", + " result_u = dmd_u[curr_window].predict(ts[i])\n", + " dmd_prediction_timer.Stop()\n", + "\n", + " dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + " u_gf.SetFromTrueDofs(dmd_solution_u)\n", + " dmd_visit_dc.SetCycle(i)\n", + " dmd_visit_dc.SetTime(ts[i])\n", + " dmd_visit_dc.Save()\n", + "\n", + " if ((i % windowNumSamples == 0) and (i < len(ts)-1)):\n", + " curr_window += 1\n", + "\n", + "dmd_prediction_timer.Start()\n", + "result_u = dmd_u[curr_window].predict(t_final)\n", + "dmd_prediction_timer.Stop()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f2e2fd2-9c9d-4329-a422-f129bde64b37", + "metadata": {}, + "outputs": [], + "source": [ + "# 11. Calculate the relative error between the DMD final solution and the true solution.\n", + "dmd_solution_u = mfem.Vector(result_u.getData(), result_u.dim())\n", + "\n", + "diff_u = mfem.Vector(u.Size())\n", + "mfem.subtract_vector(dmd_solution_u, u, diff_u)\n", + "tot_diff_norm_u = np.sqrt(mfem.InnerProduct(diff_u, diff_u))\n", + "tot_true_solution_u_norm = np.sqrt(mfem.InnerProduct(u, u))\n", + "\n", + "print(\"Relative error of DMD solution (u) at t_final: %f is %.3E\" % (t_final, tot_diff_norm_u / tot_true_solution_u_norm))\n", + "print(\"Elapsed time for solving FOM: %e second\\n\" % fom_timer.duration)\n", + "print(\"Elapsed time for training DMD: %e second\\n\" % dmd_training_timer.duration)\n", + "print(\"Elapsed time for predicting DMD: %e second\\n\" % dmd_prediction_timer.duration)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/prom/nonlinear_elasticity_global_rom.ipynb b/examples/prom/nonlinear_elasticity_global_rom.ipynb new file mode 100644 index 0000000..d2d6952 --- /dev/null +++ b/examples/prom/nonlinear_elasticity_global_rom.ipynb @@ -0,0 +1,1534 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "01fba0d2-8d5a-4a79-9fa3-8c9a8bc13087", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import io\n", + "import sys\n", + "import time\n", + "try:\n", + " import mfem.par as mfem\n", + "except ModuleNotFoundError:\n", + " msg = \"PyMFEM is not installed yet. Install PyMFEM:\\n\"\n", + " msg += \"\\tgit clone https://github.com/mfem/PyMFEM.git\\n\"\n", + " msg += \"\\tcd PyMFEM\\n\"\n", + " msg += \"\\tpython3 setup.py install --with-parallel\\n\"\n", + " raise ModuleNotFoundError(msg)\n", + "\n", + "from ctypes import c_double\n", + "from mfem.par import intArray, add_vector, subtract_vector\n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum\n", + "from scipy.special import erfc\n", + "\n", + "sys.path.append(\"../../build\")\n", + "import pylibROM.linalg as linalg\n", + "import pylibROM.hyperreduction as hyper\n", + "import pylibROM.mfem as mfem_support\n", + "from pylibROM.mfem import ComputeCtAB\n", + "from pylibROM.python_utils import StopWatch" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75401dd7-36fd-4efb-9ba9-74bcf154d10f", + "metadata": {}, + "outputs": [], + "source": [ + "class HyperelasticOperator(mfem.PyTimeDependentOperator):\n", + " def __init__(self, fespace, ess_tdof_list_, visc, mu, K):\n", + " mfem.PyTimeDependentOperator.__init__(self, 2*fespace.TrueVSize(), 0.0)\n", + "\n", + " rel_tol = 1e-8\n", + " skip_zero_entries = 0\n", + " ref_density = 1.0\n", + "\n", + " self.ess_tdof_list = ess_tdof_list_\n", + " self.z = mfem.Vector(self.Height() // 2)\n", + " self.z2 = mfem.Vector(self.Height() // 2)\n", + " self.H_sp = mfem.Vector(self.Height() // 2)\n", + " self.dvxdt_sp = mfem.Vector(self.Height() // 2)\n", + " self.fespace = fespace\n", + " self.viscosity = visc\n", + "\n", + " M = mfem.ParBilinearForm(fespace)\n", + " S = mfem.ParBilinearForm(fespace)\n", + " H = mfem.ParNonlinearForm(fespace)\n", + " self.M = M\n", + " self.H = H\n", + " self.S = S\n", + "\n", + " rho = mfem.ConstantCoefficient(ref_density)\n", + " M.AddDomainIntegrator(mfem.VectorMassIntegrator(rho))\n", + " M.Assemble(skip_zero_entries)\n", + " M.Finalize(skip_zero_entries)\n", + " self.Mmat = M.ParallelAssemble()\n", + " self.Mmat.EliminateRowsCols(self.ess_tdof_list)\n", + "\n", + " M_solver = mfem.CGSolver(fespace.GetComm())\n", + " M_prec = mfem.HypreSmoother()\n", + " M_solver.iterative_mode = False\n", + " M_solver.SetRelTol(rel_tol)\n", + " M_solver.SetAbsTol(0.0)\n", + " M_solver.SetMaxIter(30)\n", + " M_solver.SetPrintLevel(0)\n", + " M_prec.SetType(mfem.HypreSmoother.Jacobi)\n", + " M_solver.SetPreconditioner(M_prec)\n", + " M_solver.SetOperator(self.Mmat)\n", + "\n", + " self.M_solver = M_solver\n", + " self.M_prec = M_prec\n", + "\n", + " model = mfem.NeoHookeanModel(mu, K)\n", + " H.AddDomainIntegrator(mfem.HyperelasticNLFIntegrator(model))\n", + " H.SetEssentialTrueDofs(self.ess_tdof_list)\n", + " self.model = model\n", + "\n", + " visc_coeff = mfem.ConstantCoefficient(visc)\n", + " S.AddDomainIntegrator(mfem.VectorDiffusionIntegrator(visc_coeff))\n", + " S.Assemble(skip_zero_entries)\n", + " S.Finalize(skip_zero_entries)\n", + " self.Smat = mfem.HypreParMatrix()\n", + " S.FormSystemMatrix(self.ess_tdof_list, self.Smat)\n", + "\n", + " def Mult(self, vx, dvx_dt):\n", + " sc = self.Height() // 2\n", + " v = mfem.Vector(vx, 0, sc)\n", + " x = mfem.Vector(vx, sc, sc)\n", + " dv_dt = mfem.Vector(dvx_dt, 0, sc)\n", + " dx_dt = mfem.Vector(dvx_dt, sc, sc)\n", + "\n", + " self.H.Mult(x, self.z)\n", + " self.H_sp.Assign(self.z)\n", + "\n", + " if (self.viscosity != 0.0):\n", + " self.Smat.Mult(v, self.z2)\n", + " self.z += self.z2\n", + " \n", + " self.z.Neg()\n", + " self.M_solver.Mult(self.z, dv_dt)\n", + "\n", + " dx_dt.Assign(v) # this changes dvx_dt\n", + " self.dvxdt_sp.Assign(dvx_dt)\n", + "\n", + " def ElasticEnergy(self, x):\n", + " return self.H.GetEnergy(x)\n", + "\n", + " def KineticEnergy(self, v):\n", + " from mpi4py import MPI\n", + " local_energy = 0.5*self.M.InnerProduct(v, v)\n", + " energy = MPI.COMM_WORLD.allreduce(local_energy, op=MPI.SUM)\n", + " return energy\n", + "\n", + " def GetElasticEnergyDensity(self, x, w):\n", + " w_coeff = ElasticEnergyCoefficient(self.model, x)\n", + " w.ProjectCoefficient(w_coeff)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdce55b2-05b2-42ec-bc28-b6e546be509b", + "metadata": {}, + "outputs": [], + "source": [ + "class RomOperator(mfem.PyTimeDependentOperator):\n", + " # RomOperator::RomOperator(HyperelasticOperator* fom_,\n", + " # HyperelasticOperator* fomSp_, const int rvdim_, const int rxdim_,\n", + " # const int hdim_, CAROM::SampleMeshManager* smm_, const Vector* v0_,\n", + " # const Vector* x0_, const Vector v0_fom_, const CAROM::Matrix* V_v_,\n", + " # const CAROM::Matrix* V_x_, const CAROM::Matrix* U_H_,\n", + " # const CAROM::Matrix* Hsinv_, const int myid, const bool oversampling_,\n", + " # const bool hyperreduce_, const bool x_base_only_)\n", + " def __init__(self, fom_, fomSp_, rvdim_, rxdim_, hdim_, smm_, v0_,\n", + " x0_, v0_fom_, V_v_, V_x_, U_H_, Hsinv_, myid, oversampling_,\n", + " hyperreduce_, x_base_only_):\n", + " mfem.PyTimeDependentOperator.__init__(self, rxdim_ + rvdim_, 0.0)\n", + " self.fom = fom_\n", + " self.fomSp = fomSp_\n", + " self.rxdim, self.rvdim, self.hdim = rxdim_, rvdim_, hdim_\n", + " self.x0, self.v0 = x0_, v0_\n", + " self.v0_fom = v0_fom_\n", + " self.smm = smm_\n", + " self.nsamp_H = smm_.GetNumVarSamples(\"H\")\n", + " self.V_x, self.V_v = V_x_, V_v_\n", + " self.U_H, self.Hsinv = U_H_, Hsinv_\n", + " self.zN = linalg.Vector(max(self.nsamp_H, 1), False)\n", + " self.zX = linalg.Vector(max(self.nsamp_H, 1), False)\n", + " self.oversampling = oversampling_\n", + " self.M_hat_solver = mfem.CGSolver(fom_.fespace.GetComm())\n", + " self.z = mfem.Vector(self.Height() // 2)\n", + " self.hyperreduce = hyperreduce_\n", + " self.x_base_only = x_base_only_\n", + "\n", + " if (myid == 0):\n", + " self.V_v_sp = linalg.Matrix(self.fomSp.Height() // 2, self.rvdim, False)\n", + " self.V_x_sp = linalg.Matrix(self.fomSp.Height() // 2, self.rxdim, False)\n", + " # TODO(kevin): we might need to initialize for non-root processes as well.\n", + "\n", + " # Gather distributed vectors\n", + " if (self.x_base_only):\n", + " self.smm.GatherDistributedMatrixRows(\"X\", self.V_v, self.rvdim, self.V_v_sp)\n", + " else:\n", + " self.smm.GatherDistributedMatrixRows(\"V\", self.V_v, self.rvdim, self.V_v_sp)\n", + "\n", + " self.smm.GatherDistributedMatrixRows(\"X\", self.V_x, self.rxdim, self.V_x_sp)\n", + "\n", + " # Create V_vTU_H, for hyperreduction\n", + " self.V_vTU_H = self.V_v.transposeMult(self.U_H)\n", + "\n", + " self.S_hat = linalg.Matrix(self.rvdim, self.rvdim, False)\n", + " self.S_hat_v0 = linalg.Vector(self.rvdim, False)\n", + " self.S_hat_v0_temp = mfem.Vector(self.v0_fom.Size())\n", + " self.S_hat_v0_temp_librom = linalg.Vector(self.S_hat_v0_temp.GetDataArray(), True, False)\n", + " self.M_hat = linalg.Matrix(self.rvdim, self.rvdim, False)\n", + " self.M_hat_inv = linalg.Matrix(self.rvdim, self.rvdim, False)\n", + "\n", + " # Create S_hat\n", + " ComputeCtAB(self.fom.Smat, self.V_v, self.V_v, self.S_hat)\n", + "\n", + " # Apply S_hat to the initial velocity and store\n", + " self.fom.Smat.Mult(self.v0_fom, self.S_hat_v0_temp)\n", + " self.V_v.transposeMult(self.S_hat_v0_temp_librom, self.S_hat_v0)\n", + "\n", + " # Create M_hat\n", + " ComputeCtAB(self.fom.Mmat, self.V_v, self.V_v, self.M_hat)\n", + "\n", + " # Invert M_hat and store\n", + " self.M_hat.inverse(self.M_hat_inv)\n", + "\n", + " if (myid == 0):\n", + " self.spdim = self.fomSp.Height() # Reduced height\n", + "\n", + " self.zH = mfem.Vector(self.spdim // 2) # Samples of H\n", + "\n", + " # Allocate auxillary vectors\n", + " self.z.SetSize(self.spdim // 2)\n", + " self.z_v = mfem.Vector(self.spdim // 2)\n", + " self.z_x = mfem.Vector(self.spdim // 2)\n", + " self.z_librom = linalg.Vector(self.z.GetDataArray(), False, False)\n", + " self.z_v_librom = linalg.Vector(self.z_v.GetDataArray(), False, False)\n", + " self.z_x_librom = linalg.Vector(self.z_x.GetDataArray(), False, False)\n", + "\n", + " # This is for saving the recreated predictions\n", + " self.psp = mfem.Vector(self.spdim)\n", + " self.psp_librom = linalg.Vector(self.psp.GetDataArray(), False, False)\n", + "\n", + " # Define sub-vectors of psp.\n", + " self.psp_x = mfem.Vector(self.psp.GetData(), self.spdim // 2)\n", + " self.psp_v = mfem.Vector(self.psp, self.spdim // 2, self.spdim // 2)\n", + "\n", + " self.psp_v_librom = linalg.Vector(self.psp_v.GetDataArray(), False, False)\n", + "\n", + " if (not self.hyperreduce):\n", + " fdim = self.fom.Height() # Unreduced height\n", + "\n", + " self.z.SetSize(fdim // 2)\n", + " self.z_v.SetSize(fdim // 2)\n", + " self.z_x.SetSize(fdim // 2)\n", + " self.z_librom = linalg.Vector(self.z.GetDataArray(), False, False)\n", + " self.z_v_librom = linalg.Vector(self.z_v.GetDataArray(), True, False)\n", + " self.z_x_librom = linalg.Vector(self.z_x.GetDataArray(), True, False)\n", + "\n", + " # This is for saving the recreated predictions\n", + " self.pfom = mfem.Vector(fdim)\n", + " self.pfom_librom = linalg.Vector(self.pfom.GetDataArray(), False, False)\n", + "\n", + " # Define sub-vectors of pfom.\n", + " self.pfom_x = mfem.Vector(self.pfom.GetData(), fdim // 2)\n", + " self.pfom_v = mfem.Vector(self.pfom, fdim // 2, fdim // 2)\n", + " self.zfom_x = mfem.Vector(fdim / 2)\n", + " self.zfom_x_librom = linalg.Vector(self.zfom_x.GetDataArray(), True, False)\n", + "\n", + " self.pfom_v_librom = linalg.Vector(self.pfom_v.GetDataArray(), True, False)\n", + "\n", + " def Mult(self, vx, dvx_dt):\n", + " if (self.hyperreduce):\n", + " self.Mult_Hyperreduced(vx, dvx_dt)\n", + " else:\n", + " self.Mult_FullOrder(vx, dvx_dt)\n", + "\n", + " def Mult_Hyperreduced(self, vx, dvx_dt):\n", + " # Check that the sizes match\n", + " assert((vx.Size() == self.rvdim + self.rxdim) and (dvx_dt.Size() == self.rvdim + self.rxdim))\n", + "\n", + " # Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt\n", + " v = mfem.Vector(vx, 0, self.rvdim)\n", + " v_librom = linalg.Vector(vx.GetDataArray()[:self.rvdim], False, False)\n", + " x_librom = linalg.Vector(vx.GetDataArray()[self.rvdim:], False, False)\n", + " dv_dt = mfem.Vector(dvx_dt, 0, self.rvdim)\n", + " dx_dt = mfem.Vector(dvx_dt, self.rvdim, self.rxdim)\n", + " dv_dt_librom = linalg.Vector(dv_dt.GetDataArray(), False, False)\n", + " dx_dt_librom = linalg.Vector(dx_dt.GetDataArray(), False, False)\n", + "\n", + " # Lift the x-, and v-vectors\n", + " # I.e. perform v = v0 + V_v v^, where v^ is the input\n", + " self.V_v_sp.mult(v_librom, self.z_v_librom)\n", + " self.V_x_sp.mult(x_librom, self.z_x_librom)\n", + "\n", + " add_vector(self.z_v, self.v0, self.psp_v) # Store liftings\n", + " add_vector(self.z_x, self.x0, self.psp_x)\n", + "\n", + " # Hyperreduce H\n", + " # Apply H to x to get zH\n", + " self.fomSp.H.Mult(self.psp_x, self.zH)\n", + "\n", + " # Sample the values from zH\n", + " self.smm.GetSampledValues(\"H\", self.zH, self.zN)\n", + "\n", + " # Apply inverse H-basis\n", + " if (self.oversampling):\n", + " self.Hsinv.transposeMult(self.zN, self.zX)\n", + " else:\n", + " self.Hsinv.mult(self.zN, self.zX)\n", + "\n", + " # Multiply by V_v^T * U_H\n", + " self.V_vTU_H.mult(self.zX, self.z_librom)\n", + "\n", + " if (self.fomSp.viscosity != 0.0):\n", + " # Apply S^, the reduced S operator, to v\n", + " self.S_hat.multPlus(self.z_librom, v_librom, 1.0)\n", + " self.z_librom += self.S_hat_v0\n", + "\n", + " self.z.Neg() # z = -z\n", + " self.M_hat_inv.mult(self.z_librom, dv_dt_librom) # to invert reduced mass matrix operator.\n", + "\n", + " self.V_x_sp.transposeMult(self.psp_v_librom, dx_dt_librom)\n", + "\n", + " def Mult_FullOrder(self, vx, dvx_dt):\n", + " # Check that the sizes match\n", + " assert((vx.Size() == self.rvdim + self.rxdim) and (dvx_dt.Size() == self.rvdim + self.rxdim))\n", + "\n", + " # Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt\n", + " v = mfem.Vector(vx, 0, self.rvdim)\n", + " v_librom = linalg.Vector(vx.GetDataArray()[:self.rvdim], False, False)\n", + " x_librom = linalg.Vector(vx.GetDataArray()[self.rvdim:], False, False)\n", + " dv_dt = mfem.Vector(dvx_dt, 0, self.rvdim)\n", + " dx_dt = mfem.Vector(dvx_dt, self.rvdim, self.rxdim)\n", + " dv_dt_librom = linalg.Vector(dv_dt.GetDataArray(), False, False)\n", + " dx_dt_librom = linalg.Vector(dx_dt.GetDataArray(), False, False)\n", + "\n", + " # Lift the x-, and v-vectors\n", + " # I.e. perform v = v0 + V_v v^, where v^ is the input\n", + " self.V_x.mult(x_librom, self.z_x_librom)\n", + " self.V_v.mult(v_librom, self.z_v_librom)\n", + "\n", + " add_vector(self.z_x, self.x0, self.pfom_x) # Store liftings\n", + " add_vector(self.z_v, self.v0, self.pfom_v)\n", + "\n", + " # Apply H to x to get z\n", + " self.fom.H.Mult(self.pfom_x, self.zfom_x)\n", + "\n", + " self.V_v.transposeMult(self.zfom_x_librom, self.z_librom)\n", + "\n", + " if (self.fom.viscosity != 0.0):\n", + " # Apply S^, the reduced S operator, to v\n", + " self.S_hat.multPlus(self.z_librom, v_librom, 1.0)\n", + " self.z_librom += self.S_hat_v0\n", + "\n", + " self.z.Neg() # z = -z\n", + " self.M_hat_inv.mult(self.z_librom, dv_dt_librom) # to invert reduced mass matrix operator.\n", + "\n", + " self.V_x.transposeMult(self.pfom_v_librom, dx_dt_librom)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d840cd1-dd27-4b92-baab-7f6f04cda191", + "metadata": {}, + "outputs": [], + "source": [ + "# /** Function representing the elastic energy density for the given hyperelastic\n", + "# model+deformation. Used in HyperelasticOperator::GetElasticEnergyDensity. */\n", + "class ElasticEnergyCoefficient(mfem.PyCoefficient):\n", + " def __init__(self, model, x):\n", + " self.x = x\n", + " self.model = model\n", + " self.J = mfem.DenseMatrix()\n", + " mfem.PyCoefficient.__init__(self)\n", + "\n", + " def Eval(self, T, ip):\n", + " self.model.SetTransformation(T)\n", + " self.x.GetVectorGradient(T, self.J)\n", + " return self.model.EvalW(self.J)/(self.J.Det())\n", + "\n", + "class InitialDeformationIC1(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " from copy import deepcopy\n", + " y = deepcopy(x)\n", + " return y\n", + "\n", + "class InitialVelocityIC1(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " global s\n", + " s_eff = s / 80.0\n", + "\n", + " v = np.zeros(len(x))\n", + " v[-1] = -s_eff * sin(s * x[0])\n", + " return v\n", + "\n", + "class InitialDeformationIC2(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " global s\n", + " s_eff = s\n", + "\n", + " from copy import deepcopy\n", + " y = deepcopy(x)\n", + " y[-1] = x[-1] + 0.25 * x[0] * s_eff\n", + " return y\n", + "\n", + "class InitialVelocityIC2(mfem.VectorPyCoefficient):\n", + " def EvalValue(self, x):\n", + " v = np.zeros(len(x))\n", + " return v" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c75ebc53-7557-4051-89e7-433f76d93cda", + "metadata": {}, + "outputs": [], + "source": [ + "# def GetFirstColumns(N, A):\n", + "# S = linalg.Matrix(A.numRows(), min(N, A.numColumns()), A.distributed())\n", + "# for (int i = 0; i < S->numRows(); ++i)\n", + "# {\n", + "# for (int j = 0; j < S->numColumns(); ++j)\n", + "# (*S)(i, j) = (*A)(i, j);\n", + "# }\n", + "\n", + "# # delete A; # TODO: find a good solution for this.\n", + "# return S;\n", + "\n", + "# TODO: move this to the library?\n", + "def BasisGeneratorFinalSummary(bg, energyFraction, cutoff, cutoffOutputPath):\n", + " rom_dim = bg.getSpatialBasis().numColumns()\n", + " sing_vals = bg.getSingularValues()\n", + "\n", + " assert(rom_dim <= sing_vals.dim())\n", + "\n", + " sum = 0.0\n", + " for sv in range(sing_vals.dim()):\n", + " sum += sing_vals[sv]\n", + "\n", + " energy_fractions = [ 0.9999999, 0.999999, 0.99999, 0.9999, 0.999, 0.99, 0.9 ]\n", + " reached_cutoff = False\n", + "\n", + " outfile = open(cutoffOutputPath, 'w')\n", + "\n", + " partialSum = 0.0\n", + " for sv in range(sing_vals.dim()):\n", + " partialSum += sing_vals[sv]\n", + " for i in range(len(energy_fractions)-1, -1, -1):\n", + " if (partialSum / sum > energy_fractions[i]):\n", + " outfile.write(\"For energy fraction: %.5E, take first %d of %d basis vectors\" % (energy_fractions[i], sv+1, sing_vals.dim()))\n", + " energy_fractions.pop(-1)\n", + " else:\n", + " break\n", + "\n", + " if ((not reached_cutoff) and (partialSum / sum > energyFraction)):\n", + " cutoff = sv + 1\n", + " reached_cutoff = True\n", + "\n", + " if (not reached_cutoff): cutoff = sing_vals.dim()\n", + " outfile.write(\"Take first %d of %d basis vectors\" % (cutoff, sing_vals.dim()))\n", + " outfile.close()\n", + "\n", + " return cutoff" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca328166-308a-470e-89c4-cbe4a4780d18", + "metadata": {}, + "outputs": [], + "source": [ + "def MergeBasis(dimFOM, nparam, max_num_snapshots, name):\n", + " if (not (nparam > 0)):\n", + " raise RuntimeError(\"Must specify a positive number of parameter sets\")\n", + "\n", + " update_right_SV = False\n", + " isIncremental = False\n", + "\n", + " options = linalg.Options(dimFOM, nparam * max_num_snapshots, 1, update_right_SV)\n", + " generator = linalg.BasisGenerator(options, isIncremental, \"basis\" + name)\n", + "\n", + " for paramID in range(nparam):\n", + " snapshot_filename = \"basis%d_%s_snapshot\" % (paramID, name)\n", + " generator.loadSamples(snapshot_filename, \"snapshot\")\n", + "\n", + " generator.endSamples() # save the merged basis file\n", + "\n", + " cutoff = 0\n", + " cutoff = BasisGeneratorFinalSummary(generator, 0.9999999, cutoff, \"mergedSV_\" + name + \".txt\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "600ebc91-e18d-48c0-af2d-d746b80246c3", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: remove this by making online computation serial?\n", + "def BroadcastUndistributedRomVector(v):\n", + " N = v.dim()\n", + " assert(N > 0)\n", + "\n", + " from copy import deepcopy\n", + " d = deepcopy(v.getData())\n", + " \n", + " from mpi4py import MPI\n", + " MPI.COMM_WORLD.Bcast([d, MPI.DOUBLE], root=0)\n", + "\n", + " for i in range(N):\n", + " v[i] = d[i]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fa0e027-9445-41eb-8267-d331a1a0187d", + "metadata": {}, + "outputs": [], + "source": [ + "def visualize(out, mesh, deformed_nodes, field,\n", + " field_name='', init_vis=False):\n", + " nodes = deformed_nodes\n", + " owns_nodes = 0\n", + "\n", + " nodes, owns_nodes = mesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " out.send_text(\"parallel \" + str(mesh.GetNRanks()) + \" \" + str(mesh.GetMyRank()))\n", + " out.send_solution(mesh, field)\n", + "\n", + " nodes, owns_nodes = mesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " if (init_vis):\n", + " out.send_text(\"window_size 800 800\")\n", + " out.send_text(\"window_title '\" + field_name)\n", + " if (mesh.SpaceDimension() == 2):\n", + " out.send_text(\"view 0 0\") # view from top\n", + " out.send_text(\"keys jl\") # turn off perspective and light\n", + " out.send_text(\"keys cm\") # show colorbar and mesh\n", + " # update value-range; keep mesh-extents fixed\n", + " out.send_text(\"autoscale value\")\n", + " out.send_text(\"pause\")\n", + " out.flush()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0ada436-dd19-4def-964c-2c3213cef476", + "metadata": {}, + "outputs": [], + "source": [ + "# Scaling factor for parameterization\n", + "s = 1.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d202353-56af-4f58-8356-b2ea40168404", + "metadata": {}, + "outputs": [], + "source": [ + "# 1. Initialize MPI.\n", + "from mpi4py import MPI\n", + "comm = MPI.COMM_WORLD\n", + "myid = comm.Get_rank()\n", + "num_procs = comm.Get_size()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09a2ae2a-9c47-46e1-803f-0044e7c819fd", + "metadata": {}, + "outputs": [], + "source": [ + "# 2. Parse command-line options.\n", + "from mfem.common.arg_parser import ArgParser\n", + "parser = ArgParser(description=\"Projection ROM - MFEM nonlinear elasticity equation example.\")\n", + "parser.add_argument(\"-m\", \"--mesh\",\n", + " action='store', dest='mesh_file', default=\"../data/beam-quad.mesh\", type=str,\n", + " help=\"Mesh file to use.\")\n", + "parser.add_argument(\"-rs\", \"--refine-serial\",\n", + " action='store', dest='ser_ref_levels', default=2, type=int,\n", + " help=\"Number of times to refine the mesh uniformly in serial.\")\n", + "parser.add_argument(\"-rp\", \"--refine-parallel\",\n", + " action='store', dest='par_ref_levels', default=0, type=int,\n", + " help=\"Number of times to refine the mesh uniformly in parallel.\")\n", + "parser.add_argument(\"-o\", \"--order\",\n", + " action='store', default=2, type=int,\n", + " help=\"Order (degree) of the finite elements.\")\n", + "parser.add_argument(\"-s\", \"--ode-solver\",\n", + " action='store', dest='ode_solver_type', default=14, type=int,\n", + " help=\"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\\n\\t\"\n", + " \" 11 - Forward Euler, 12 - RK2,\\n\\t\"\n", + " \" 13 - RK3 SSP, 14 - RK4.\"\n", + " \" 22 - Implicit Midpoint Method,\\n\\t\"\n", + " \" 23 - SDIRK23 (A-stable), 24 - SDIRK34\")\n", + "parser.add_argument(\"-tf\", \"--t-final\",\n", + " action='store', default=15.0, type=float,\n", + " help=\"Final time; start time is 0.\")\n", + "parser.add_argument(\"-dt\", \"--time-step\",\n", + " action='store', dest='dt', default=0.03, type=float,\n", + " help=\"Time step.\")\n", + "parser.add_argument(\"-v\", \"--viscosity\",\n", + " action='store', dest='visc', default=1e-2, type=float,\n", + " help=\"Viscosity coefficient.\")\n", + "parser.add_argument(\"-mu\", \"--shear-modulus\",\n", + " action='store', dest='mu', default=0.25, type=float,\n", + " help=\"Shear modulus in the Neo-Hookean hyperelastic model.\")\n", + "parser.add_argument(\"-K\", \"--bulk-modulus\",\n", + " action='store', dest='K', default=5.0, type=float,\n", + " help=\"Bulk modulus in the Neo-Hookean hyperelastic model.\")\n", + "parser.add_argument(\"-alrtol\", \"--adaptive-lin-rtol\",\n", + " action='store_true', default=True, dest='adaptive_lin_rtol',\n", + " help=\"Enable adaptive linear solver rtol.\")\n", + "parser.add_argument(\"-no-alrtol\", \"--no-adaptive-lin-rtol\",\n", + " action='store_false', dest='adaptive_lin_rtol',\n", + " help=\"Disable adaptive linear solver rtol.\")\n", + "parser.add_argument(\"-vis\", \"--visualization\",\n", + " action='store_true', default=True, dest='visualization',\n", + " help=\"Enable GLVis visualization.\")\n", + "parser.add_argument(\"-no-vis\", \"--no-visualization\",\n", + " action='store_false', dest='visualization',\n", + " help=\"Disable GLVis visualization.\")\n", + "parser.add_argument(\"-visit\", \"--visit-datafiles\",\n", + " action='store_true', default=False, dest='visit',\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-no-visit\", \"--no-visit-datafiles\",\n", + " action='store_false', dest='visit',\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-vs\", \"--visualization-steps\",\n", + " action='store', dest='vis_steps', default=1, type=int,\n", + " help=\"Visualize every n-th timestep.\")\n", + "parser.add_argument(\"-ns\", \"--nset\",\n", + " action='store', dest='nsets', default=0, type=int,\n", + " help=\"Number of parametric snapshot sets\")\n", + "parser.add_argument(\"-offline\", \"--offline\", \n", + " action='store_true', dest='offline', default=False,\n", + " help=\"Enable the offline phase.\")\n", + "parser.add_argument(\"-no-offline\", \"--no-offline\",\n", + " action='store_false', dest='offline',\n", + " help=\"Disable the offline phase.\")\n", + "parser.add_argument(\"-online\", \"--online\", \n", + " action='store_true', dest='online', default=False,\n", + " help=\"Enable the online phase.\")\n", + "parser.add_argument(\"-no-online\", \"--no-online\",\n", + " action='store_false', dest='online',\n", + " help=\"Disable the online phase.\")\n", + "parser.add_argument(\"-merge\", \"--merge\", \n", + " action='store_true', dest='merge', default=False,\n", + " help=\"Enable the merge phase.\")\n", + "parser.add_argument(\"-no-merge\", \"--no-merge\",\n", + " action='store_false', dest='merge',\n", + " help=\"Disable the merge phase.\")\n", + "parser.add_argument(\"-sopt\", \"--sopt\", \n", + " action='store_true', dest='use_sopt', default=False,\n", + " help=\"Use S-OPT sampling instead of DEIM for the hyperreduction.\")\n", + "parser.add_argument(\"-no-sopt\", \"--no-sopt\",\n", + " action='store_false', dest='use_sopt',\n", + " help=\"(disable) Use S-OPT sampling instead of DEIM for the hyperreduction.\")\n", + "parser.add_argument(\"-nsr\", \"--nsr\",\n", + " action='store', dest='num_samples_req', default=-1, type=int,\n", + " help=\"number of samples we want to select for the sampling algorithm.\")\n", + "parser.add_argument(\"-rxdim\", \"--rxdim\",\n", + " action='store', default=-1, type=int,\n", + " help=\"Basis dimension for displacement solution space.\")\n", + "parser.add_argument(\"-rvdim\", \"--rvdim\",\n", + " action='store', default=-1, type=int,\n", + " help=\"Basis dimension for velocity solution space.\")\n", + "parser.add_argument(\"-hdim\", \"--hdim\",\n", + " action='store', default=-1, type=int,\n", + " help=\"Basis dimension for the nonlinear term.\")\n", + "parser.add_argument(\"-id\", \"--id\",\n", + " action='store', dest='id_param', default=0, type=int,\n", + " help=\"Parametric index\")\n", + "parser.add_argument(\"-hyp\", \"--hyperreduce\",\n", + " action='store_true', dest='hyperreduce', default=False,\n", + " help=\"Enable Hyper reduce nonlinear term\")\n", + "parser.add_argument(\"-no-hyp\", \"--no-hyperreduce\",\n", + " action='store_false', dest='hyperreduce',\n", + " help=\"Disable Hyper reduce nonlinear term\")\n", + "parser.add_argument(\"-xbo\", \"--xbase-only\",\n", + " action='store_true', dest='x_base_only', default=False,\n", + " help=\"Use the displacement (X) basis to approximate velocity.\")\n", + "parser.add_argument(\"-no-xbo\", \"--not-xbase-only\",\n", + " action='store_false', dest='x_base_only',\n", + " help=\"not use the displacement (X) basis to approximate velocity.\")\n", + "parser.add_argument(\"-def-ic\", \"--deformation-ic\",\n", + " action='store_true', dest='def_ic', default=False,\n", + " help=\"Use a deformation initial condition. Default is velocity IC.\")\n", + "parser.add_argument(\"-vel-ic\", \"--velocity-ic\",\n", + " action='store_false', dest='def_ic',\n", + " help=\"Use velocity initial condition. Default is velocity IC.\")\n", + "parser.add_argument(\"-sc\", \"--scaling\",\n", + " action='store', dest='s', default=1.0, type=float,\n", + " help=\"Scaling factor for initial condition.\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7914daf4-3fd7-4dfd-a301-54e336ad2a4b", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run:\n", + "\n", + "# Offline phase:\n", + "# args = parser.parse_args(\"-offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 0.90 -id 0\".split())\n", + "# args = parser.parse_args(\"-offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.10 -id 1\".split())\n", + "\n", + "# Merge phase:\n", + "# args = parser.parse_args(\"-merge -ns 2 -dt 0.01 -tf 5.0\".split())\n", + "\n", + "# Create FOM comparison data:\n", + "# args = parser.parse_args(\"-offline -dt 0.01 -tf 5.0 -s 14 -vs 100 -sc 1.00 -id 2\".split())\n", + "\n", + "# Online phase with full sampling:\n", + "# args = parser.parse_args(\"-online -dt 0.01 -tf 5.0 -s 14 -vs 100 -hyp -rvdim 40 -rxdim 10 -hdim 71 -nsr 1170 -sc 1.00\".split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4836af8-b05d-48d1-983f-78b5ed7f5bae", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "if (myid == 0): parser.print_options(args)\n", + "\n", + "mesh_file = args.mesh_file\n", + "ser_ref_levels = args.ser_ref_levels\n", + "par_ref_levels = args.par_ref_levels\n", + "order = args.order\n", + "ode_solver_type = args.ode_solver_type\n", + "vis_steps = args.vis_steps\n", + "t_final = args.t_final\n", + "dt = args.dt\n", + "visc = args.visc\n", + "mu = args.mu\n", + "K = args.K\n", + "adaptive_lin_rtol = args.adaptive_lin_rtol\n", + "visualization = args.visualization\n", + "visit = args.visit\n", + "def_ic = args.def_ic\n", + "global s\n", + "s = args.s\n", + "\n", + "# ROM parameters\n", + "offline = args.offline\n", + "merge = args.merge\n", + "online = args.online\n", + "use_sopt = args.use_sopt\n", + "hyperreduce = args.hyperreduce\n", + "x_base_only = args.x_base_only\n", + "num_samples_req = args.num_samples_req\n", + "nsets = args.nsets\n", + "id_param = args.id_param\n", + "\n", + "# number of basis vectors to use\n", + "rxdim = args.rxdim\n", + "rvdim = args.rvdim\n", + "hdim = args.hdim\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a6a0ac1-2fed-440b-9650-68c9a0f809fa", + "metadata": {}, + "outputs": [], + "source": [ + "check = (offline and (not merge) and (not online)) or ((not offline) and merge and (not online)) or ((not offline) and (not merge) and online)\n", + "if not check:\n", + " raise RuntimeError(\"only one of offline, merge, or online must be true!\")\n", + "\n", + "solveTimer, totalTimer = StopWatch(), StopWatch()\n", + "totalTimer.Start()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e2b3c1b-abd6-4582-8ee9-07b933bc0857", + "metadata": {}, + "outputs": [], + "source": [ + "# // 3. Read the serial mesh from the given mesh file on all processors. We can\n", + "# // handle triangular, quadrilateral, tetrahedral and hexahedral meshes\n", + "# // with the same code.\n", + "mesh = mfem.Mesh(mesh_file, 1, 1)\n", + "dim = mesh.Dimension()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a73719f3-db65-4c73-aee0-5c0fd0b42d49", + "metadata": {}, + "outputs": [], + "source": [ + "# // 4. Define the ODE solver used for time integration. Several implicit\n", + "# // singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as\n", + "# // explicit Runge-Kutta methods are available.\n", + "# if ode_solver_type == 1:\n", + "# ode_solver = BackwardEulerSolver()\n", + "# elif ode_solver_type == 2:\n", + "# ode_solver = mfem.SDIRK23Solver(2)\n", + "# elif ode_solver_type == 3:\n", + "# ode_solver = mfem.SDIRK33Solver()\n", + "if ode_solver_type == 11:\n", + " ode_solver = ForwardEulerSolver()\n", + "elif ode_solver_type == 12:\n", + " ode_solver = mfem.RK2Solver(0.5)\n", + "elif ode_solver_type == 13:\n", + " ode_solver = mfem.RK3SSPSolver()\n", + "elif ode_solver_type == 14:\n", + " ode_solver = mfem.RK4Solver()\n", + "elif ode_solver_type == 15:\n", + " ode_solver = mfem.GeneralizedAlphaSolver(0.5)\n", + "# elif ode_solver_type == 22:\n", + "# ode_solver = mfem.ImplicitMidpointSolver()\n", + "# elif ode_solver_type == 23:\n", + "# ode_solver = mfem.SDIRK23Solver()\n", + "# elif ode_solver_type == 24:\n", + "# ode_solver = mfem.SDIRK34Solver()\n", + "else:\n", + " if myid == 0:\n", + " print(\"Unknown ODE solver type: \" + str(ode_solver_type))\n", + " sys.exit()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2cba58c-75bf-4969-a307-62bcd03aa405", + "metadata": {}, + "outputs": [], + "source": [ + "# // 5. Refine the mesh in serial to increase the resolution. In this example\n", + "# // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is\n", + "# // a command-line parameter.\n", + "for lev in range(ser_ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "057129e0-34ab-4269-8b6e-446b1d3d3a69", + "metadata": {}, + "outputs": [], + "source": [ + "# // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine\n", + "# // this mesh further in parallel to increase the resolution. Once the\n", + "# // parallel mesh is defined, the serial mesh can be deleted.\n", + "pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh)\n", + "del mesh\n", + "\n", + "for lev in range(par_ref_levels):\n", + " pmesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd101ad2-234e-4afd-8bd3-f8683802bc79", + "metadata": {}, + "outputs": [], + "source": [ + "# // 7. Define the parallel vector finite element spaces representing the mesh\n", + "# // deformation x_gf, the velocity v_gf, and the initial configuration,\n", + "# // x_ref. Define also the elastic energy density, w_gf, which is in a\n", + "# // discontinuous higher-order space. Since x and v are integrated in time\n", + "# // as a system, we group them together in block vector vx, on the unique\n", + "# // parallel degrees of freedom, with offsets given by array true_offset.\n", + "fe_coll = mfem.H1_FECollection(order, dim)\n", + "fespace = mfem.ParFiniteElementSpace(pmesh, fe_coll, dim)\n", + "glob_size = fespace.GlobalTrueVSize()\n", + "if (myid == 0):\n", + " print('Number of velocity/deformation unknowns: ' + str(glob_size))\n", + "\n", + "true_size = fespace.TrueVSize()\n", + "true_offset = mfem.intArray(3)\n", + "true_offset[0] = 0\n", + "true_offset[1] = true_size\n", + "true_offset[2] = 2*true_size\n", + "\n", + "vx = mfem.BlockVector(true_offset)\n", + "\n", + "v_gf = mfem.ParGridFunction(fespace)\n", + "x_gf = mfem.ParGridFunction(fespace)\n", + "v_gf.MakeTRef(fespace, vx, true_offset[0])\n", + "x_gf.MakeTRef(fespace, vx, true_offset[1])\n", + "# v_gf.MakeTRef(&fespace, vx,\n", + "# true_offset[0]); // Associate a new FiniteElementSpace and new true-dof data with the GridFunction.\n", + "# x_gf.MakeTRef(&fespace, vx, true_offset[1]);\n", + "\n", + "x_ref = mfem.ParGridFunction(fespace)\n", + "pmesh.GetNodes(x_ref)\n", + "\n", + "w_fec = mfem.L2_FECollection(order + 1, dim)\n", + "w_fespace = mfem.ParFiniteElementSpace(pmesh, w_fec)\n", + "w_gf = mfem.ParGridFunction(w_fespace)\n", + "\n", + "# Basis params\n", + "update_right_SV = False\n", + "isIncremental = False\n", + "basisFileName = \"basis\" + str(id_param)\n", + "max_num_snapshots = int(t_final / dt) + 2\n", + "\n", + "# The merge phase\n", + "if (merge):\n", + " totalTimer.Reset()\n", + " totalTimer.Start()\n", + "\n", + " # Merge bases\n", + " if (not x_base_only):\n", + " MergeBasis(true_size, nsets, max_num_snapshots, \"V\")\n", + "\n", + " MergeBasis(true_size, nsets, max_num_snapshots, \"X\")\n", + " MergeBasis(true_size, nsets, max_num_snapshots, \"H\")\n", + "\n", + " totalTimer.Stop()\n", + " if (myid == 0):\n", + " print(\"Elapsed time for merging and building ROM basis: %e second\\n\" % totalTimer.duration)\n", + "\n", + " sys.exit(0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe196f09-8531-41b3-8086-78d829492d3c", + "metadata": {}, + "outputs": [], + "source": [ + "# // 8. Set the initial conditions for v_gf, x_gf and vx, and define the\n", + "# // boundary conditions on a beam-like mesh (see description above).\n", + "if (def_ic):\n", + " velo = InitialVelocityIC2(dim)\n", + "else:\n", + " velo = InitialVelocityIC1(dim)\n", + "\n", + "v_gf.ProjectCoefficient(velo)\n", + "v_gf.SetTrueVector()\n", + "\n", + "if (def_ic):\n", + " deform = InitialDeformationIC2(dim)\n", + "else:\n", + " deform = InitialDeformationIC1(dim)\n", + "\n", + "x_gf.ProjectCoefficient(deform)\n", + "x_gf.SetTrueVector()\n", + "\n", + "v_gf.SetFromTrueVector()\n", + "x_gf.SetFromTrueVector()\n", + "\n", + "v_gf.GetTrueDofs(vx.GetBlock(0))\n", + "x_gf.GetTrueDofs(vx.GetBlock(1))\n", + "\n", + "ess_bdr = mfem.intArray(fespace.GetMesh().bdr_attributes.Max())\n", + "ess_bdr.Assign(0)\n", + "ess_bdr[0] = 1 # boundary attribute 1 (index 0) is fixed\n", + "\n", + "ess_tdof_list = mfem.intArray()\n", + "fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list)\n", + "\n", + "# Store initial vx\n", + "vx0 = mfem.BlockVector(vx)\n", + "vx_diff = mfem.BlockVector(vx)\n", + "\n", + "# Reduced order solution\n", + "# Vector* wMFEM = 0;\n", + "# CAROM::Vector* w = 0;\n", + "# CAROM::Vector* w_v = 0;\n", + "# CAROM::Vector* w_x = 0;\n", + "\n", + "# Initialize reconstructed solution\n", + "v_rec = mfem.Vector(v_gf.GetTrueVector())\n", + "x_rec = mfem.Vector(x_gf.GetTrueVector())\n", + "\n", + "v_rec_librom = linalg.Vector(v_rec.GetDataArray(), True, False)\n", + "x_rec_librom = linalg.Vector(x_rec.GetDataArray(), True, False)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1985e0b5-273f-4378-90c2-227e434f040e", + "metadata": {}, + "outputs": [], + "source": [ + "# // 9. Initialize the hyperelastic operator, the GLVis visualization and print\n", + "# // the initial energies.\n", + "oper = HyperelasticOperator(fespace, ess_tdof_list, visc, mu, K)\n", + "soper = None\n", + "\n", + "# Fill dvdt and dxdt\n", + "dvxdt = mfem.Vector(true_size * 2)\n", + "dvdt = mfem.Vector(dvxdt, 0, true_size)\n", + "dxdt = mfem.Vector(dvxdt, true_size, true_size)\n", + "\n", + "if (visualization):\n", + " vishost = \"localhost\"\n", + " visport = 19916\n", + " vis_v = mfem.socketstream(vishost, visport)\n", + " vis_v.precision(8)\n", + " visualize(vis_v, pmesh, x_gf, v_gf, \"Velocity\", True)\n", + " # // Make sure all ranks have sent their 'v' solution before initiating\n", + " # // another set of GLVis connections (one from each rank):\n", + " pmesh.GetComm().Barrier()\n", + " vis_w = mfem.socketstream(vishost, visport)\n", + " if (vis_w.good()):\n", + " oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + " vis_w.precision(8)\n", + " visualize(vis_w, pmesh, x_gf, w_gf, \"Elastic energy density\", True)\n", + "\n", + "# // Create data collection for solution output: either VisItDataCollection for\n", + "# // ascii data files, or SidreDataCollection for binary data files.\n", + "dc = None\n", + "if (visit):\n", + " if (offline):\n", + " dc = mfem.VisItDataCollection(\"nlelast-fom\", pmesh)\n", + " else:\n", + " dc = mfem.VisItDataCollection(\"nlelast-rom\", pmesh)\n", + "\n", + " dc.SetPrecision(8)\n", + " # // To save the mesh using MFEM's parallel mesh format:\n", + " # // dc->SetFormat(DataCollection::PARALLEL_FORMAT);\n", + " dc.RegisterField(\"x\", x_gf)\n", + " dc.RegisterField(\"v\", v_gf)\n", + " dc.SetCycle(0)\n", + " dc.SetTime(0.0)\n", + " dc.Save()\n", + "\n", + "ee0 = oper.ElasticEnergy(x_gf)\n", + "ke0 = oper.KineticEnergy(v_gf)\n", + "\n", + "if (myid == 0):\n", + " print(\"initial elastic energy (EE) = %e\" % ee0)\n", + " print(\"initial kinetic energy (KE) = %e\" % ke0)\n", + " print(\"initial total energy (TE) = %e\" % (ee0 + ke0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ccd82203-79c3-4fa6-b7aa-8a18acb6e969", + "metadata": {}, + "outputs": [], + "source": [ + "# // 10. Create pROM object.\n", + "if (offline):\n", + " options = linalg.Options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV)\n", + "\n", + " if (not x_base_only):\n", + " basis_generator_v = linalg.BasisGenerator(options, isIncremental, basisFileName + \"_V\")\n", + "\n", + " basis_generator_x = linalg.BasisGenerator(options, isIncremental, basisFileName + \"_X\")\n", + "\n", + " basis_generator_H = linalg.BasisGenerator(options, isIncremental, basisFileName + \"_H\")\n", + "\n", + "# RomOperator* romop = 0;\n", + "\n", + "# const CAROM::Matrix* BV_librom = 0;\n", + "# const CAROM::Matrix* BX_librom = 0;\n", + "# const CAROM::Matrix* H_librom = 0;\n", + "# const CAROM::Matrix* Hsinv = 0;\n", + "\n", + "nsamp_H = -1\n", + "\n", + "# CAROM::SampleMeshManager* smm = nullptr;" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aeb3e590-cff0-4256-aa49-0b319a3eaa44", + "metadata": {}, + "outputs": [], + "source": [ + "# // 11. Initialize ROM operator\n", + "# // I guess most of this should be done on id =0\n", + "if (online):\n", + " # Read bases\n", + " if (x_base_only):\n", + " readerV = linalg.BasisReader(\"basisX\") # The basis for v uses the x basis instead.\n", + " rvdim = rxdim\n", + " else:\n", + " readerV = linalg.BasisReader(\"basisV\")\n", + "\n", + " BV_librom = readerV.getSpatialBasis(0.0)\n", + "\n", + " if (rvdim == -1): # Change rvdim\n", + " rvdim = BV_librom.numColumns()\n", + " else:\n", + " BV_librom = BV_librom.getFirstNColumns(rvdim)\n", + "\n", + " assert(BV_librom.numRows() == true_size)\n", + "\n", + " if (myid == 0):\n", + " print(\"reduced V dim = %d\\n\" % rvdim)\n", + "\n", + " readerX = linalg.BasisReader(\"basisX\")\n", + " BX_librom = readerX.getSpatialBasis(0.0)\n", + "\n", + " if (rxdim == -1): # Change rxdim\n", + " rxdim = BX_librom.numColumns()\n", + " else:\n", + " BX_librom = BX_librom.getFirstNColumns(rxdim)\n", + "\n", + " assert(BX_librom.numRows() == true_size)\n", + "\n", + " if (myid == 0):\n", + " print(\"reduced X dim = %d\\n\" % rxdim)\n", + "\n", + " # Hyper reduce H\n", + " readerH = linalg.BasisReader(\"basisH\")\n", + " H_librom = readerH.getSpatialBasis(0.0)\n", + "\n", + " # Compute sample points\n", + " if (hdim == -1):\n", + " hdim = H_librom.numColumns()\n", + "\n", + " assert(H_librom.numColumns() >= hdim)\n", + "\n", + " if (H_librom.numColumns() > hdim):\n", + " H_librom = H_librom.getFirstNColumns(hdim)\n", + "\n", + " if (myid == 0):\n", + " print(\"reduced H dim = %d\\n\" % hdim)\n", + "\n", + " # vector num_sample_dofs_per_proc(num_procs);\n", + "\n", + " if (num_samples_req != -1):\n", + " nsamp_H = num_samples_req\n", + " else:\n", + " nsamp_H = hdim\n", + "\n", + " Hsinv = linalg.Matrix(nsamp_H, hdim, False)\n", + " # vector sample_dofs(nsamp_H);\n", + " if (use_sopt):\n", + " if (myid == 0):\n", + " print(\"Using S_OPT sampling\\n\")\n", + " sample_dofs, num_sample_dofs_per_proc = hyper.S_OPT(H_librom, hdim, Hsinv,\n", + " myid, num_procs, nsamp_H)\n", + " elif (nsamp_H != hdim):\n", + " if (myid == 0):\n", + " print(\"Using GNAT sampling\\n\")\n", + " sample_dofs, num_sample_dofs_per_proc = hyper.GNAT(H_librom, hdim, Hsinv,\n", + " myid, num_procs, nsamp_H)\n", + " else:\n", + " if (myid == 0):\n", + " print(\"Using DEIM sampling\\n\")\n", + " sample_dofs, num_sample_dofs_per_proc = hyper.DEIM(H_librom, hdim, Hsinv,\n", + " myid, num_procs)\n", + "\n", + " # Construct sample mesh\n", + " nspaces = 1\n", + " spfespace = [None] * nspaces\n", + " spfespace[0] = fespace\n", + "\n", + " # ParFiniteElementSpace* sp_XV_space;\n", + " smm = mfem_support.SampleMeshManager(spfespace)\n", + "\n", + " # vector sample_dofs_empty;\n", + " num_sample_dofs_per_proc_empty = [0] * num_procs\n", + "\n", + " smm.RegisterSampledVariable(\"V\", 0, sample_dofs,\n", + " num_sample_dofs_per_proc)\n", + " smm.RegisterSampledVariable(\"X\", 0, sample_dofs,\n", + " num_sample_dofs_per_proc)\n", + " smm.RegisterSampledVariable(\"H\", 0, sample_dofs,\n", + " num_sample_dofs_per_proc)\n", + "\n", + " smm.ConstructSampleMesh()\n", + "\n", + " w = linalg.Vector(rxdim + rvdim, False)\n", + " w_v = linalg.Vector(rvdim, False)\n", + " w_x = linalg.Vector(rxdim, False)\n", + " w.fill(0.0)\n", + "\n", + " # Note that some of this could be done only on the ROM solver process, but it is tricky, since RomOperator assembles Bsp in parallel.\n", + " wMFEM = mfem.Vector(w.getData(), rxdim + rvdim)\n", + "\n", + " # Initial conditions\n", + " # Vector* w_v0 = 0;\n", + " # Vector* w_x0 = 0;\n", + "\n", + " sp_size = 0\n", + "\n", + " if (myid == 0):\n", + " # NOTE(kevin): SampleMeshManager::GetSampleFESpace returns a pointer to a ParFiniteElementSpace,\n", + " # which is binded via SWIG, not pybind11.\n", + " # We need a SWIG object 'shell' to wrap the c++ pointer we return from pybind11 side.\n", + " # The following instantiation is simply creating SWIG object,\n", + " # to which GetSampleFESpace will return the c++ pointer.\n", + " # Instantiation can be of any type, since we only need the 'shell'.\n", + " sp_XV_space = mfem.ParFiniteElementSpace(pmesh, fe_coll, dim)\n", + " # NOTE(kevin): Unlike c++ libROM, GetSampleFESpace returns a deep copy.\n", + " smm.GetSampleFESpace(0, sp_XV_space)\n", + "\n", + " sp_size = sp_XV_space.TrueVSize()\n", + " sp_offset = mfem.intArray(3)\n", + " sp_offset[0] = 0\n", + " sp_offset[1] = sp_size\n", + " sp_offset[2] = 2*sp_size\n", + "\n", + " # Initialize sp_p with initial conditions.\n", + " sp_vx = mfem.BlockVector(sp_offset)\n", + " sp_v_gf, sp_x_gf = mfem.ParGridFunction(), mfem.ParGridFunction()\n", + "\n", + " # // 12. Set the initial conditions for v_gf, x_gf and vx, and define the\n", + " # // boundary conditions on a beam-like mesh (see description above).\n", + "\n", + " sp_v_gf.MakeTRef(sp_XV_space, sp_vx,\n", + " sp_offset[0]) # Associate a new FiniteElementSpace and new true-dof data with the GridFunction.\n", + " sp_x_gf.MakeTRef(sp_XV_space, sp_vx, sp_offset[1])\n", + "\n", + " # VectorFunctionCoefficient* velo = 0;\n", + " # VectorFunctionCoefficient* deform = 0;\n", + "\n", + " if (def_ic):\n", + " velo = InitialVelocityIC2(dim)\n", + " else:\n", + " velo = InitialVelocityIC1(dim)\n", + "\n", + " sp_v_gf.ProjectCoefficient(velo)\n", + " sp_v_gf.SetTrueVector()\n", + "\n", + " if (def_ic):\n", + " deform = InitialDeformationIC2(dim)\n", + " else:\n", + " deform = InitialDeformationIC1(dim)\n", + "\n", + " sp_x_gf.ProjectCoefficient(deform)\n", + " sp_x_gf.SetTrueVector()\n", + "\n", + " sp_v_gf.SetFromTrueVector()\n", + " sp_x_gf.SetFromTrueVector()\n", + "\n", + " # Get initial conditions\n", + " w_v0 = mfem.Vector(sp_v_gf.GetTrueVector())\n", + " w_x0 = mfem.Vector(sp_x_gf.GetTrueVector())\n", + "\n", + " # Convert essential boundary list from FOM mesh to sample mesh\n", + " # Create binary list where 1 means essential boundary element, 0 means nonessential.\n", + " Ess_mat = linalg.Matrix(true_size, 1, True)\n", + " for i in range(true_size):\n", + " Ess_mat[i,0] = 0.\n", + " for j in range(ess_tdof_list.Size()):\n", + " if (ess_tdof_list[j] == i ):\n", + " Ess_mat[i,0] = 1.\n", + "\n", + " # Project binary FOM list onto sampling space\n", + " MPI.COMM_WORLD.bcast(sp_size, root=0)\n", + " Ess_mat_sp = linalg.Matrix(sp_size, 1, False)\n", + " smm.GatherDistributedMatrixRows(\"X\", Ess_mat, 1, Ess_mat_sp)\n", + "\n", + " # Count number of true elements in new matrix\n", + " num_ess_sp = 0\n", + "\n", + " for i in range(sp_size):\n", + " if (Ess_mat_sp[i,0] == 1):\n", + " num_ess_sp += 1\n", + "\n", + " # Initialize essential dof list in sampling space\n", + " ess_tdof_list_sp = mfem.intArray(num_ess_sp)\n", + "\n", + " # Add indices to list\n", + " ctr = 0\n", + " for i in range(sp_size):\n", + " if (Ess_mat_sp[i,0] == 1):\n", + " ess_tdof_list_sp[ctr] = i\n", + " ctr += 1\n", + "\n", + " if (myid == 0):\n", + " # Define operator in sample space\n", + " soper = HyperelasticOperator(sp_XV_space, ess_tdof_list_sp, visc, mu, K)\n", + "\n", + " if (hyperreduce):\n", + " romop = RomOperator(oper, soper, rvdim, rxdim, hdim, smm, w_v0, w_x0,\n", + " vx0.GetBlock(0), BV_librom, BX_librom, H_librom, Hsinv, myid,\n", + " (num_samples_req != -1), hyperreduce, x_base_only)\n", + " else:\n", + " romop = RomOperator(oper, soper, rvdim, rxdim, hdim, smm,\n", + " vx0.GetBlock(0), vx0.GetBlock(1), vx0.GetBlock(0),\n", + " BV_librom, BX_librom, H_librom, Hsinv, myid,\n", + " (num_samples_req != -1), hyperreduce, x_base_only)\n", + "\n", + " # Print lifted initial energies\n", + " BroadcastUndistributedRomVector(w)\n", + "\n", + " for i in range(rvdim):\n", + " w_v[i] = w[i]\n", + "\n", + " for i in range(rxdim):\n", + " w_x[i] = w[rvdim + i]\n", + "\n", + " romop.V_v.mult(w_v, v_rec_librom)\n", + " romop.V_x.mult(w_x, x_rec_librom)\n", + "\n", + " v_rec += vx0.GetBlock(0)\n", + " x_rec += vx0.GetBlock(1)\n", + "\n", + " v_gf.SetFromTrueDofs(v_rec)\n", + " x_gf.SetFromTrueDofs(x_rec)\n", + "\n", + " ee = oper.ElasticEnergy(x_gf)\n", + " ke = oper.KineticEnergy(v_gf)\n", + "\n", + " if (myid == 0):\n", + " print(\"Lifted initial energies, EE = %.5e, KE = %.5e, ΔTE = %.5e\" % (ee, ke, (ee + ke) - (ee0 + ke0)))\n", + "\n", + " ode_solver.Init(romop)\n", + "else:\n", + " # fom\n", + " ode_solver.Init(oper)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89a3719c-df0d-468c-823a-fa53eb80ca8e", + "metadata": {}, + "outputs": [], + "source": [ + "# 13. Perform time-integration\n", + "# (looping over the time iterations, ti, with a time-step dt).\n", + "# (taking samples and storing it into the pROM object)\n", + "\n", + "t = 0.0\n", + "ts = []\n", + "oper.SetTime(t)\n", + "\n", + "last_step = False\n", + "ti = 1\n", + "while (not last_step):\n", + " dt_real = min(dt, t_final - t)\n", + "\n", + " if (online):\n", + " if (myid == 0):\n", + " solveTimer.Start()\n", + " t, dt = ode_solver.Step(wMFEM, t, dt_real)\n", + " solveTimer.Stop()\n", + "\n", + " MPI.COMM_WORLD.bcast(t, root=0)\n", + " else:\n", + " solveTimer.Start()\n", + " t, dt = ode_solver.Step(vx, t, dt_real)\n", + " solveTimer.Stop()\n", + "\n", + " last_step = (t >= t_final - 1e-8 * dt)\n", + "\n", + " if (offline):\n", + " if (basis_generator_x.isNextSample(t) or (not x_base_only) and basis_generator_v.isNextSample(t)):\n", + " dvxdt = oper.dvxdt_sp.GetData()\n", + " vx_diff = mfem.BlockVector(vx)\n", + " vx_diff -= vx0\n", + "\n", + " # Take samples\n", + " if ((not x_base_only) and basis_generator_v.isNextSample(t)):\n", + " basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray(), t, dt)\n", + " basis_generator_v.computeNextSampleTime(vx_diff.GetBlock(0).GetDataArray(), dvdt.GetDataArray(), t)\n", + " basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)\n", + "\n", + " if (basis_generator_x.isNextSample(t)):\n", + " basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray(), t, dt)\n", + " basis_generator_x.computeNextSampleTime(vx_diff.GetBlock(1).GetDataArray(), dxdt.GetDataArray(), t)\n", + "\n", + " if (x_base_only):\n", + " basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)\n", + "\n", + " if (last_step or ((ti % vis_steps) == 0)):\n", + " if (online):\n", + " BroadcastUndistributedRomVector(w)\n", + "\n", + " for i in range(rvdim):\n", + " w_v[i] = w[i]\n", + "\n", + " for i in range(rxdim):\n", + " w_x[i] = w[rvdim + i]\n", + "\n", + " romop.V_v.mult(w_v, v_rec_librom)\n", + " romop.V_x.mult(w_x, x_rec_librom)\n", + "\n", + " v_rec += vx0.GetBlock(0)\n", + " x_rec += vx0.GetBlock(1)\n", + "\n", + " v_gf.SetFromTrueDofs(v_rec)\n", + " x_gf.SetFromTrueDofs(x_rec)\n", + "\n", + " else:\n", + " v_gf.SetFromTrueVector()\n", + " x_gf.SetFromTrueVector()\n", + "\n", + " ee = oper.ElasticEnergy(x_gf)\n", + " ke = oper.KineticEnergy(v_gf)\n", + "\n", + " if (myid == 0):\n", + " print(\"step %d, t = %f, EE = %.5e, KE = %.5e, ΔTE = %.5e\" % (ti, t, ee, ke, (ee + ke) - (ee0 + ke0)))\n", + "\n", + " if (visualization):\n", + " visualize(vis_v, pmesh, x_gf, v_gf)\n", + " if (vis_w):\n", + " oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + " visualize(vis_w, pmesh, x_gf, w_gf)\n", + "\n", + " if (visit):\n", + " nodes = x_gf\n", + " owns_nodes = 0\n", + " pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " dc.SetCycle(ti)\n", + " dc.SetTime(t)\n", + " dc.Save()\n", + "\n", + " ti += 1\n", + "# timestep loop\n", + "\n", + "if (myid == 0):\n", + " print(\"Elapsed time for time integration loop %.5e\" % solveTimer.duration)\n", + "\n", + "velo_name = \"velocity_s%f.%06d\" % (s, myid)\n", + "pos_name = \"position_s%f.%06d\" % (s, myid)\n", + "\n", + "if (offline):\n", + " # Sample final solution, to prevent extrapolation in ROM between the last sample and the end of the simulation.\n", + " dvxdt = oper.dvxdt_sp.GetData()\n", + " vx_diff = mfem.BlockVector(vx)\n", + " vx_diff -= vx0\n", + "\n", + " # Take samples\n", + " if (not x_base_only):\n", + " basis_generator_v.takeSample(vx_diff.GetBlock(0).GetDataArray(), t, dt)\n", + " basis_generator_v.writeSnapshot()\n", + " del basis_generator_v\n", + "\n", + " basis_generator_H.takeSample(oper.H_sp.GetDataArray(), t, dt)\n", + " basis_generator_H.writeSnapshot()\n", + " del basis_generator_H\n", + "\n", + " basis_generator_x.takeSample(vx_diff.GetBlock(1).GetDataArray(), t, dt)\n", + " basis_generator_x.writeSnapshot()\n", + " del basis_generator_x\n", + "\n", + " # 14. Save the displaced mesh, the velocity and elastic energy.\n", + " nodes = x_gf\n", + " owns_nodes = 0\n", + " pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " mesh_name = \"deformed.%06d\" % myid\n", + " ee_name = \"elastic_energy.%06d\" % myid\n", + "\n", + " pmesh.Print(mesh_name, 8)\n", + " pmesh.SwapNodes(nodes, owns_nodes)\n", + "\n", + " np.savetxt(velo_name, vx.GetBlock(0).GetDataArray(), fmt='%.16e')\n", + " np.savetxt(pos_name, vx.GetBlock(1).GetDataArray(), fmt='%.16e')\n", + "\n", + " with open(ee_name, 'w') as fid:\n", + " ee_ofs = io.StringIO()\n", + " ee_ofs.precision = 8\n", + " oper.GetElasticEnergyDensity(x_gf, w_gf)\n", + " w_gf.Save(ee_ofs)\n", + " fid.write(ee_ofs.getvalue())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86bfe4ce-b210-419e-94e1-8e7f6c61f08c", + "metadata": {}, + "outputs": [], + "source": [ + "# 15. Calculate the relative error between the ROM final solution and the true solution.\n", + "if (online):\n", + " # Initialize FOM solution\n", + " v_fom = mfem.Vector(v_rec.Size())\n", + " x_fom = mfem.Vector(x_rec.Size())\n", + "\n", + " # Open and load file\n", + " v_fom.Load(velo_name, v_rec.Size())\n", + "\n", + " x_fom.Load(pos_name, x_rec.Size())\n", + "\n", + " # Get difference vector\n", + " diff_v = mfem.Vector(v_rec.Size())\n", + " diff_x = mfem.Vector(x_rec.Size())\n", + "\n", + " subtract_vector(v_rec, v_fom, diff_v)\n", + " subtract_vector(x_rec, x_fom, diff_x)\n", + "\n", + " # Get norms\n", + " tot_diff_norm_v = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff_v, diff_v))\n", + " tot_diff_norm_x = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff_x, diff_x))\n", + "\n", + " tot_v_fom_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, v_fom, v_fom))\n", + " tot_x_fom_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, x_fom, x_fom))\n", + "\n", + " if (myid == 0):\n", + " print(\"Relative error of ROM position (x) at t_final: %f is %.8e\" % (t_final, tot_diff_norm_x / tot_x_fom_norm))\n", + " print(\"Relative error of ROM velocity (v) at t_final: %f is %.8e\" % (t_final, tot_diff_norm_v / tot_v_fom_norm))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e80c4ee5-b277-4ff2-8490-d7d8ad61cb55", + "metadata": {}, + "outputs": [], + "source": [ + "# 16. Free the used memory.\n", + "del ode_solver\n", + "del pmesh\n", + "\n", + "totalTimer.Stop()\n", + "if (myid == 0):\n", + " print(\"Elapsed time for entire simulation %.5e\" % totalTimer.duration)\n", + "\n", + "MPI.Finalize()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/prom/poisson_global_rom.ipynb b/examples/prom/poisson_global_rom.ipynb new file mode 100644 index 0000000..e565a8b --- /dev/null +++ b/examples/prom/poisson_global_rom.ipynb @@ -0,0 +1,692 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "eecdf3ba-95dd-456f-93d6-271b41e2fffc", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import io\n", + "import sys\n", + "import time\n", + "try:\n", + " import mfem.par as mfem\n", + "except ModuleNotFoundError:\n", + " msg = \"PyMFEM is not installed yet. Install PyMFEM:\\n\"\n", + " msg += \"\\tgit clone https://github.com/mfem/PyMFEM.git\\n\"\n", + " msg += \"\\tcd PyMFEM\\n\"\n", + " msg += \"\\tpython3 setup.py install --with-parallel\\n\"\n", + " raise ModuleNotFoundError(msg)\n", + "\n", + "from ctypes import c_double\n", + "from mfem.par import intArray\n", + "from os.path import expanduser, join, dirname\n", + "import numpy as np\n", + "from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "215cf768-c626-4b33-a92b-2e2b9942e926", + "metadata": {}, + "outputs": [], + "source": [ + "# pyMFEM does not provide mfem::StopWatch.\n", + "# there is not really a similar stopwatch package in python.. (seriously?)\n", + "class StopWatch:\n", + " import time\n", + " duration = 0.0\n", + " start_time = 0.0\n", + " stop_time = 0.0\n", + " running = False\n", + "\n", + " def __init__(self):\n", + " self.Reset()\n", + " return\n", + "\n", + " def Start(self):\n", + " assert(not self.running)\n", + " self.start_time = time.time()\n", + " self.running = True\n", + " return\n", + "\n", + " def Stop(self):\n", + " assert(self.running)\n", + " self.stop_time = time.time()\n", + " self.duration += self.stop_time - self.start_time\n", + " self.running = False\n", + " return\n", + "\n", + " def Reset(self):\n", + " self.duration = 0.0\n", + " self.start_time = 0.0\n", + " self.stop_time = 0.0\n", + " self.running = False\n", + " return\n", + "\n", + "sys.path.append(\"../../build\")\n", + "import pylibROM.linalg as libROM\n", + "from pylibROM.mfem import ComputeCtAB" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78bab678-1d1f-4dd1-b01f-200d9b131272", + "metadata": {}, + "outputs": [], + "source": [ + "from mpi4py import MPI\n", + "comm = MPI.COMM_WORLD\n", + "myid = comm.Get_rank()\n", + "num_procs = comm.Get_size()\n", + "\n", + "from mfem.common.arg_parser import ArgParser\n", + "parser = ArgParser(description=\"Projection ROM - MFEM Poisson equation example.\")\n", + "parser.add_argument('-m', '--mesh',\n", + " default='../data/star.mesh',\n", + " action='store', type=str,\n", + " help='Mesh file to use.')\n", + "parser.add_argument('-o', '--order',\n", + " action='store', default=1, type=int,\n", + " help=\"Finite element order (polynomial degree) or -1 for isoparametric space.\")\n", + "parser.add_argument(\"-id\", \"--id\",\n", + " action='store', default=0, type=int, help=\"Parametric id\")\n", + "parser.add_argument(\"-ns\", \"--nset\",\n", + " action='store', default=0, type=int, help=\"Number of parametric snapshot sets\")\n", + "parser.add_argument(\"-sc\", \"--static-condensation\",\n", + " action='store_true', default=False,\n", + " help=\"Enable static condensation.\")\n", + "parser.add_argument(\"-pa\", \"--partial-assembly\",\n", + " action='store_true', default=False,\n", + " help=\"Enable Partial Assembly.\")\n", + "parser.add_argument(\"-f\", \"--frequency\",\n", + " action='store', default=1.0, type=float,\n", + " help=\"Set the frequency for the exact solution.\")\n", + "parser.add_argument(\"-cf\", \"--coefficient\",\n", + " action='store', default=1.0, type=float,\n", + " help=\"Coefficient.\")\n", + "parser.add_argument(\"-d\", \"--device\",\n", + " action='store', default='cpu', type=str,\n", + " help=\"Device configuration string, see Device::Configure().\")\n", + "parser.add_argument(\"-visit\", \"--visit-datafiles\",\n", + " action='store_true', default=True,\n", + " help=\"Save data files for VisIt (visit.llnl.gov) visualization.\")\n", + "parser.add_argument(\"-vis\", \"--visualization\",\n", + " action='store_true', default=True,\n", + " help=\"Enable or disable GLVis visualization.\")\n", + "parser.add_argument(\"-fom\", \"--fom\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the fom phase.\")\n", + "parser.add_argument(\"-offline\", \"--offline\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the offline phase.\")\n", + "parser.add_argument(\"-online\", \"--online\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the online phase.\")\n", + "parser.add_argument(\"-merge\", \"--merge\",\n", + " action='store_true', default=False,\n", + " help=\"Enable or disable the merge phase.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69d9cc67-50f8-4bf2-88ba-f1423d650354", + "metadata": {}, + "outputs": [], + "source": [ + "# Sample run:\n", + "\n", + "# Offline phase:\n", + "# args = parser.parse_args(\"-offline -f 1.0 -id 0\".split())\n", + "# args = parser.parse_args(\"-offline -f 1.1 -id 1\".split())\n", + "# args = parser.parse_args(\"-offline -f 1.2 -id 2\".split()) \n", + "\n", + "# Merge phase: \n", + "# args = parser.parse_args(\"-merge -ns 3\".split())\n", + "\n", + "# FOM run for error calculation:\n", + "# args = parser.parse_args(\"-fom -f 1.15\".split())\n", + "\n", + "# Online phase:\n", + "# args = parser.parse_args(\"-online -f 1.15\".split()) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c2522ee-3a0c-4eda-8ac1-e035538c6d6d", + "metadata": {}, + "outputs": [], + "source": [ + "args = parser.parse_args([])\n", + "if (myid == 0): \n", + " parser.print_options(args)\n", + "\n", + "mesh_file = os.path.abspath(os.path.join('../data', args.mesh))\n", + "freq = args.frequency\n", + "fom = args.fom\n", + "offline = args.offline\n", + "online = args.online\n", + "merge = args.merge\n", + "device_config = args.device\n", + "id = args.id\n", + "order = args.order\n", + "nsets = args.nset\n", + "coef = args.coefficient\n", + "pa = args.partial_assembly\n", + "static_cond = args.static_condensation\n", + "visit = args.visit_datafiles\n", + "visualization = args.visualization\n", + "precision = 8\n", + "\n", + "kappa = freq * np.pi\n", + "if (fom):\n", + " if (not (fom and (not offline) and (not online))):\n", + " raise ValueError(\"offline and online must be turned off if fom is used.\")\n", + "else:\n", + " check = (offline and (not merge) and (not online)) \\\n", + " or ((not offline) and merge and (not online)) \\\n", + " or ((not offline) and (not merge) and online)\n", + " if (not check):\n", + " raise ValueError(\"only one of offline, merge, or online must be true!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1eb914f5-1d79-4eb9-8a05-fae57b117d96", + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Enable hardware devices such as GPUs, and programming models such as\n", + "# CUDA, OCCA, RAJA and OpenMP based on command line options.\n", + "device = mfem.Device(device_config)\n", + "if (myid == 0):\n", + " device.Print()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b00a6c0-d8be-41cd-b8fb-322947e4b245", + "metadata": {}, + "outputs": [], + "source": [ + "# 4. Read the (serial) mesh from the given mesh file on all processors. We\n", + "# can handle triangular, quadrilateral, tetrahedral, hexahedral, surface\n", + "# and volume meshes with the same code.\n", + "mesh = mfem.Mesh(mesh_file, 1, 1)\n", + "dim = mesh.Dimension()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a4c7f20-28dc-4d43-beb2-46e001bb7d3c", + "metadata": {}, + "outputs": [], + "source": [ + "# 5. Refine the serial mesh on all processors to increase the resolution. In\n", + "# this example we do 'ref_levels' of uniform refinement. We choose\n", + "# 'ref_levels' to be the largest number that gives a final mesh with no\n", + "# more than 10,000 elements.\n", + "ref_levels = int(np.floor(np.log(10000. / mesh.GetNE()) / log(2.) / dim))\n", + "for l in range(ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "99e8d104-2a29-4c8e-81b8-80eccd849f9d", + "metadata": {}, + "outputs": [], + "source": [ + "# 6. Define a parallel mesh by a partitioning of the serial mesh. Refine\n", + "# this mesh further in parallel to increase the resolution. Once the\n", + "# parallel mesh is defined, the serial mesh can be deleted.\n", + "pmesh = mfem.ParMesh(comm, mesh)\n", + "mesh.Clear()\n", + "par_ref_levels = 2\n", + "for l in range(par_ref_levels):\n", + " pmesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26a3b6a9-3467-454c-a389-8bb557946631", + "metadata": {}, + "outputs": [], + "source": [ + "# 7. Define a parallel finite element space on the parallel mesh. Here we\n", + "# use continuous Lagrange finite elements of the specified order. If\n", + "# order < 1, we instead use an isoparametric/isogeometric space.\n", + "if (order > 0):\n", + " fec = mfem.H1_FECollection(order, dim)\n", + " delete_fec = True\n", + "elif (pmesh.GetNodes()):\n", + " fec = pmesh.GetNodes().OwnFEC()\n", + " delete_fec = False\n", + " if (myid == 0):\n", + " print(\"Using isoparametric FEs: %s\" % fec.Name())\n", + "else:\n", + " fec = mfem.H1_FECollection(1, dim)\n", + " delete_fec = True\n", + "\n", + "fespace = mfem.ParFiniteElementSpace(pmesh, fec)\n", + "size = fespace.GlobalTrueVSize()\n", + "if (myid == 0):\n", + " print(\"Number of finite element unknowns: %d\" % size)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82ba06ab-bec0-43dc-9e4c-e4002bf930de", + "metadata": {}, + "outputs": [], + "source": [ + "# 8. Determine the list of true (i.e. parallel conforming) essential\n", + "# boundary dofs. In this example, the boundary conditions are defined\n", + "# by marking all the boundary attributes from the mesh as essential\n", + "# (Dirichlet) and converting them to a list of true dofs.\n", + "ess_tdof_list = mfem.intArray()\n", + "if (pmesh.bdr_attributes.Size() > 0):\n", + " ess_bdr = mfem.intArray(pmesh.bdr_attributes.Max())\n", + " ess_bdr.Assign(1)\n", + " fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2c811db-1cc3-4d3a-adac-375df80b1332", + "metadata": {}, + "outputs": [], + "source": [ + "# 9. Initiate ROM related variables\n", + "max_num_snapshots = 100\n", + "update_right_SV = False\n", + "isIncremental = False\n", + "basisName = \"basis\"\n", + "basisFileName = \"%s%d\" % (basisName, id)\n", + "solveTimer, assembleTimer, mergeTimer = StopWatch(), StopWatch(), StopWatch()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88cdd25b-5999-4c85-aa26-0e1c94a1798d", + "metadata": {}, + "outputs": [], + "source": [ + "# 10. Set BasisGenerator if offline\n", + "if (offline):\n", + " options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots, 1,\n", + " update_right_SV)\n", + " generator = libROM.BasisGenerator(options, isIncremental, basisFileName)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e426e37c-3c46-4a9d-9966-9dd1a8bbcaf9", + "metadata": {}, + "outputs": [], + "source": [ + "# 11. The merge phase\n", + "if (merge):\n", + " mergeTimer.Start()\n", + " options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots, 1,\n", + " update_right_SV)\n", + " generator = libROM.BasisGenerator(options, isIncremental, basisName)\n", + " for paramID in range(nsets):\n", + " snapshot_filename = \"%s%d_snapshot\" % (basisName, paramID)\n", + " generator.loadSamples(snapshot_filename,\"snapshot\", 5)\n", + "\n", + " generator.endSamples() # save the merged basis file\n", + " mergeTimer.Stop()\n", + " if (myid == 0):\n", + " print(\"Elapsed time for merging and building ROM basis: %e second\\n\" %\n", + " mergeTimer.duration)\n", + " del generator\n", + " del options\n", + " MPI.Finalize()\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a435a41-44dc-4f13-9190-865361a4668a", + "metadata": {}, + "outputs": [], + "source": [ + "# 12. Set up the parallel linear form b(.) which corresponds to the\n", + "# right-hand side of the FEM linear system, which in this case is\n", + "# (f,phi_i) where f is given by the function f_exact and phi_i are the\n", + "# basis functions in the finite element fespace.\n", + "assembleTimer.Start()\n", + "b = mfem.ParLinearForm(fespace)\n", + "class RightHandSide(mfem.PyCoefficient):\n", + " def EvalValue(self, x):\n", + " if (dim == 3):\n", + " return sin(kappa * (x[0] + x[1] + x[2]))\n", + " else:\n", + " return sin(kappa * (x[0] + x[1]))\n", + "f = RightHandSide()\n", + "b.AddDomainIntegrator(mfem.DomainLFIntegrator(f))\n", + "b.Assemble()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81f1a238-4462-40b6-b69b-7749df669858", + "metadata": {}, + "outputs": [], + "source": [ + "# 13. Define the solution vector x as a parallel finite element grid function\n", + "# corresponding to fespace. Initialize x with initial guess of zero,\n", + "# which satisfies the boundary conditions.\n", + "x = mfem.ParGridFunction(fespace)\n", + "x.Assign(0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea191b3a-2541-4502-905e-1eb7c11a1e9a", + "metadata": {}, + "outputs": [], + "source": [ + "# 14. Set up the parallel bilinear form a(.,.) on the finite element space\n", + "# corresponding to the Laplacian operator -Delta, by adding the Diffusion\n", + "# domain integrator.\n", + "a = mfem.ParBilinearForm(fespace)\n", + "one = mfem.ConstantCoefficient(coef)\n", + "if (pa):\n", + " a.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL)\n", + "a.AddDomainIntegrator(mfem.DiffusionIntegrator(one))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d65d3a52-38ea-48e4-a718-59015ea0836a", + "metadata": {}, + "outputs": [], + "source": [ + "# 15. Assemble the parallel bilinear form and the corresponding linear\n", + "# system, applying any necessary transformations such as: parallel\n", + "# assembly, eliminating boundary conditions, applying conforming\n", + "# constraints for non-conforming AMR, static condensation, etc.\n", + "if (static_cond):\n", + " a.EnableStaticCondensation()\n", + "a.Assemble()\n", + "\n", + "# A = mfem.OperatorPtr()\n", + "A = mfem.HypreParMatrix()\n", + "B = mfem.Vector()\n", + "X = mfem.Vector()\n", + "a.FormLinearSystem(ess_tdof_list, x, b, A, X, B)\n", + "assembleTimer.Stop()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ac0a3ce-9d1c-4cd9-9a8a-dc4457a6b342", + "metadata": {}, + "outputs": [], + "source": [ + "# 16. The offline phase\n", + "if (fom or offline):\n", + " # 17. Solve the full order linear system A X = B\n", + " prec = None\n", + " if pa:\n", + " if mfem.UsesTensorBasis(fespace):\n", + " prec = mfem.OperatorJacobiSmoother(a, ess_tdof_list)\n", + " else:\n", + " prec = mfem.HypreBoomerAMG(A)\n", + "\n", + " cg = mfem.CGSolver(comm)\n", + " cg.SetRelTol(1e-12)\n", + " cg.SetMaxIter(2000)\n", + " cg.SetPrintLevel(1)\n", + " if (prec is not None):\n", + " cg.SetPreconditioner(prec)\n", + " # cg.SetOperator(A.Ptr())\n", + " cg.SetOperator(A)\n", + " solveTimer.Start()\n", + " cg.Mult(B, X)\n", + " solveTimer.Stop()\n", + " if (prec is not None):\n", + " del prec\n", + "\n", + " # 18. take and write snapshot for ROM\n", + " if (offline):\n", + " # NOTE: mfem Vector::GetData returns a SWIG Object of type double *.\n", + " # To make it compatible with pybind11, we use ctypes to read data from the memory address.\n", + " xData = np.array((c_double * X.Size()).from_address(int(X.GetData())), copy=False) # this does not copy the data.\n", + " addSample = generator.takeSample(xData, 0.0, 0.01)\n", + " generator.writeSnapshot()\n", + " del generator\n", + " del options" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "feac8d64-4df4-4e5f-9eb6-90210fcfd3ba", + "metadata": {}, + "outputs": [], + "source": [ + "# 19. The online phase\n", + "if (online):\n", + " # 20. read the reduced basis\n", + " assembleTimer.Start()\n", + " reader = libROM.BasisReader(basisName)\n", + " spatialbasis = reader.getSpatialBasis(0.0)\n", + " numRowRB = spatialbasis.numRows()\n", + " numColumnRB = spatialbasis.numColumns()\n", + " if (myid == 0):\n", + " print(\"spatial basis dimension is %d x %d\\n\" % (numRowRB, numColumnRB))\n", + "\n", + " # libROM stores the matrix row-wise, so wrapping as a DenseMatrix in MFEM means it is transposed.\n", + " reducedBasisT = mfem.DenseMatrix(spatialbasis.getData())\n", + "\n", + " # 21. form inverse ROM operator\n", + " invReducedA = libROM.Matrix(numColumnRB, numColumnRB, False)\n", + " ComputeCtAB(A, spatialbasis, spatialbasis, invReducedA)\n", + " invReducedA.invert()\n", + "\n", + " bData = np.array((c_double * B.Size()).from_address(int(B.GetData())), copy=False)\n", + " B_carom = libROM.Vector(bData, True, False)\n", + " xData = np.array((c_double * X.Size()).from_address(int(X.GetData())), copy=False)\n", + " X_carom = libROM.Vector(xData, True, False)\n", + " reducedRHS = spatialbasis.transposeMult(B_carom)\n", + " reducedSol = libROM.Vector(numColumnRB, False)\n", + " assembleTimer.Stop()\n", + "\n", + " # 22. solve ROM\n", + " solveTimer.Start()\n", + " invReducedA.mult(reducedRHS, reducedSol)\n", + " solveTimer.Stop()\n", + "\n", + " # 23. reconstruct FOM state\n", + " spatialbasis.mult(reducedSol, X_carom)\n", + " del spatialbasis\n", + " del reducedRHS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28b56f02-5f41-48c1-998b-371d2e091f64", + "metadata": {}, + "outputs": [], + "source": [ + "# 24. Recover the parallel grid function corresponding to X. This is the\n", + "# local finite element solution on each processor.\n", + "a.RecoverFEMSolution(X, b, x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5ceb341-ee76-4611-8729-6ef57e6bc4fe", + "metadata": {}, + "outputs": [], + "source": [ + "# 25. Calculate the relative error of the ROM prediction compared to FOM\n", + "# ostringstream sol_dofs_name, sol_dofs_name_fom;\n", + "if (fom or offline):\n", + " sol_dofs_name = \"sol_dofs_fom.%06d\" % myid\n", + "if (online):\n", + " sol_dofs_name = \"sol_dofs.%06d\" % myid\n", + " sol_dofs_name_fom = \"sol_dofs_fom.%06d\" % myid\n", + "\n", + "if (online):\n", + " # Initialize FOM solution\n", + " x_fom = mfem.Vector(x.Size())\n", + "\n", + " # Open and load file\n", + " x_fom.Load(sol_dofs_name_fom, x_fom.Size())\n", + "\n", + " diff_x = mfem.Vector(x.Size())\n", + "\n", + " mfem.subtract_vector(x, x_fom, diff_x)\n", + "\n", + " # Get norms\n", + " tot_diff_norm = np.sqrt(mfem.InnerProduct(comm, diff_x, diff_x))\n", + " tot_fom_norm = np.sqrt(mfem.InnerProduct(comm, x_fom, x_fom))\n", + "\n", + " if (myid == 0):\n", + " print(\"Relative error of ROM solution = %.5E\" % (tot_diff_norm / tot_fom_norm))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17f2e60a-f9e9-40ad-8194-723afecc298e", + "metadata": {}, + "outputs": [], + "source": [ + "# 26. Save the refined mesh and the solution in parallel. This output can\n", + "# be viewed later using GLVis: \"glvis -np -m mesh -g sol\".\n", + "mesh_name = \"mesh.%06d\" % myid\n", + "sol_name = \"sol.%06d\" % myid\n", + "\n", + "pmesh.Print(mesh_name, precision)\n", + "\n", + "output = io.StringIO()\n", + "output.precision = precision\n", + "x.Save(output)\n", + "fid = open(sol_name, 'w')\n", + "fid.write(output.getvalue())\n", + "fid.close()\n", + "\n", + "xData = np.array((c_double * X.Size()).from_address(int(X.GetData())), copy=False)\n", + "np.savetxt(sol_dofs_name, xData, fmt='%.16f')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f261a03-bfed-4805-b68e-f9a0ea89ec84", + "metadata": {}, + "outputs": [], + "source": [ + "# 27. Save data in the VisIt format.\n", + "if (visit):\n", + " if (offline):\n", + " dc = mfem.VisItDataCollection(\"Example1\", pmesh)\n", + " elif (fom):\n", + " dc = mfem.VisItDataCollection(\"Example1_fom\", pmesh)\n", + " elif (online):\n", + " dc = mfem.VisItDataCollection(\"Example1_rom\", pmesh)\n", + " dc.SetPrecision(precision)\n", + " dc.RegisterField(\"solution\", x)\n", + " dc.Save()\n", + " del dc\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7f349c9-2543-4e14-98c3-d3f1af091b06", + "metadata": {}, + "outputs": [], + "source": [ + "# 28. Send the solution by socket to a GLVis server.\n", + "if visualization:\n", + " sol_sock = mfem.socketstream(\"localhost\", 19916)\n", + " if not sol_sock.good():\n", + " visualization = False\n", + " if (myid == 0):\n", + " print(\"Unable to connect to GLVis server at localhost:19916\")\n", + " print(\"GLVis visualization disabled.\")\n", + " else:\n", + " sol_sock << \"parallel \" << num_procs << \" \" << myid << \"\\n\"\n", + " sol_sock.precision(precision)\n", + " sol_sock << \"solution\\n\" << pmesh << x\n", + " sol_sock << \"pause\\n\"\n", + " sol_sock.flush()\n", + " if (myid == 0):\n", + " print(\n", + " \"GLVis visualization paused. Press space (in the GLVis window) to resume it.\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "389e3c01-6999-4fd1-b58e-b9a08bf17c79", + "metadata": {}, + "outputs": [], + "source": [ + "# 29. print timing info\n", + "if (myid == 0):\n", + " if (fom or offline):\n", + " print(\"Elapsed time for assembling FOM: %e second\\n\" % assembleTimer.duration)\n", + " print(\"Elapsed time for solving FOM: %e second\\n\" % solveTimer.duration)\n", + "\n", + " if(online):\n", + " print(\"Elapsed time for assembling ROM: %e second\\n\" % assembleTimer.duration)\n", + " print(\"Elapsed time for solving ROM: %e second\\n\" % solveTimer.duration)\n", + "\n", + "# 30. Free the used memory.\n", + "if (delete_fec):\n", + " del fec\n", + "MPI.Finalize()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}