diff --git a/examples/dmd/dg_euler.ipynb b/examples/dmd/dg_euler.ipynb index 37f4b63..aaa38a9 100644 --- a/examples/dmd/dg_euler.ipynb +++ b/examples/dmd/dg_euler.ipynb @@ -31,7 +31,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 2, "id": "cb74c4b1-9055-4b1e-af71-6317b423c0f2", "metadata": {}, "outputs": [], @@ -92,7 +92,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 3, "id": "05c0a673-0043-4696-8ec2-7c88f456fc64", "metadata": {}, "outputs": [], @@ -118,7 +118,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 4, "id": "63aa9975-f4f7-4154-b088-72506a5f7162", "metadata": {}, "outputs": [ @@ -170,7 +170,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "id": "1be6c7af-bebc-4594-a846-ad04697251b0", "metadata": {}, "outputs": [], @@ -201,7 +201,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "id": "6734fecb-26ca-47ec-8bc3-2c6ca8d2a153", "metadata": {}, "outputs": [], @@ -215,7 +215,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 7, "id": "701e36df-2b80-4bc2-81e9-297230768ea8", "metadata": {}, "outputs": [], @@ -232,7 +232,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 8, "id": "b75b4260-72f1-4ead-9310-d10701dce84b", "metadata": {}, "outputs": [ @@ -264,7 +264,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 9, "id": "fe303417-3440-4bf2-a3eb-c1f606e1a3d0", "metadata": {}, "outputs": [], @@ -296,7 +296,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, "id": "5513d593-574e-4af6-8ed6-f07f34588cce", "metadata": {}, "outputs": [], @@ -315,7 +315,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 11, "id": "05c371cc-375d-4d44-84f6-da698f45b444", "metadata": {}, "outputs": [ @@ -442,7 +442,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 12, "id": "726039c6-f8ce-425e-8849-c90144441b0a", "metadata": {}, "outputs": [], @@ -458,7 +458,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 13, "id": "e9c61b98-49f2-4248-ba7b-a54a508140c1", "metadata": {}, "outputs": [ @@ -480,7 +480,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 14, "id": "17ee664a-0c3e-4698-9e0d-eca03d8099f2", "metadata": {}, "outputs": [ @@ -2285,7 +2285,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 15, "id": "7c55cd0d-e0e7-45c1-a2f6-8cec295c2b25", "metadata": {}, "outputs": [ @@ -2334,7 +2334,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 16, "id": "a90ddf1b-c15a-4053-9705-4142cf747dcd", "metadata": {}, "outputs": [ @@ -2346,9 +2346,9 @@ "Relative error of DMD x_mom at t_final: 2.0 is 0.0002069844\n", "Relative error of DMD y_mom at t_final: 2.0 is 0.0012537761\n", "Relative error of DMD e at t_final: 2.0 is 0.0001376798\n", - "Elapsed time for solving FOM: 3.655925e+02\n", - "Elapsed time for training DMD: 4.660985e+01\n", - "Elapsed time for predicting DMD: 1.512005e-01\n" + "Elapsed time for solving FOM: 3.629454e+02\n", + "Elapsed time for training DMD: 4.648910e+01\n", + "Elapsed time for predicting DMD: 1.452651e-01\n" ] } ], diff --git a/examples/dmd/parametric_heat_conduction.ipynb b/examples/dmd/parametric_heat_conduction.ipynb new file mode 100644 index 0000000..bd68c91 --- /dev/null +++ b/examples/dmd/parametric_heat_conduction.ipynb @@ -0,0 +1,1041 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "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": 2, + "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": 3, + "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": 4, + "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": 5, + "id": "6533977a-df7b-478e-9e7c-347cebf0e1c9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Options used:\n", + " --mesh ../data/star.mesh\n", + " --refine_serial 2\n", + " --refine_parallel 1\n", + " --order 2\n", + " --ode_solver 3\n", + " --t_final 0.5\n", + " --time_step 0.01\n", + " --alpha 0.01\n", + " --kappa 0.5\n", + " --radius 0.5\n", + " --center_x 0.0\n", + " --center_y 0.0\n", + " --crv 0.9\n", + " --visualization True\n", + " --visit_datafiles False\n", + " --visualization_steps 5\n", + " --rdim -1\n", + " --offline False\n", + " --online False\n", + " --predict False\n", + " --adios2_streams False\n", + " --save False\n", + " --outputfile_name \n", + " --pw_snap False\n", + " --pwx 0\n", + " --pwy 0\n", + " --pwz 0\n" + ] + } + ], + "source": [ + "args = parser.parse_args([])\n", + "if (myid == 0):\n", + " parser.print_options(args)\n", + "\n", + "precision = 8\n", + "\n", + "save_dofs = args.save\n", + "basename = args.outputfile_name\n", + "offline = True\n", + "online = args.online\n", + "pointwiseSnapshots = args.pw_snap\n", + "rdim = 16\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 = 4\n", + "alpha = args.alpha\n", + "kappa = args.kappa\n", + "radius = 0.1\n", + "cx, cy = 0.1, 0.1\n", + "visit = True\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": 6, + "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": 7, + "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": 8, + "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": 9, + "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": 10, + "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": 11, + "id": "37df6a1d-c0e2-4df7-95a1-b24d6a2a0ae6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of temperature unknowns: 20801\n" + ] + } + ], + "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": 12, + "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": 13, + "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": 14, + "id": "91e180dd-ae83-4d48-9366-cdaadace4320", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Unable to connect to GLVis server at localhost:19916\n", + "GLVis visualization disabled.\n" + ] + } + ], + "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": 15, + "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": 16, + "id": "41c1aaeb-8866-4bf3-b2d7-ba1808887f76", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Taking snapshot at: 0.000000\n" + ] + } + ], + "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": 17, + "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": 18, + "id": "feeb685c-a84f-413b-af6e-62020337fd9a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Taking snapshot at: 0.010000\n", + "Taking snapshot at: 0.020000\n", + "Taking snapshot at: 0.030000\n", + "Taking snapshot at: 0.040000\n", + "Taking snapshot at: 0.050000\n", + "step 5, t = 0.050000\n", + "Taking snapshot at: 0.060000\n", + "Taking snapshot at: 0.070000\n", + "Taking snapshot at: 0.080000\n", + "Taking snapshot at: 0.090000\n", + "Taking snapshot at: 0.100000\n", + "step 10, t = 0.100000\n", + "Taking snapshot at: 0.110000\n", + "Taking snapshot at: 0.120000\n", + "Taking snapshot at: 0.130000\n", + "Taking snapshot at: 0.140000\n", + "Taking snapshot at: 0.150000\n", + "step 15, t = 0.150000\n", + "Taking snapshot at: 0.160000\n", + "Taking snapshot at: 0.170000\n", + "Taking snapshot at: 0.180000\n", + "Taking snapshot at: 0.190000\n", + "Taking snapshot at: 0.200000\n", + "step 20, t = 0.200000\n", + "Taking snapshot at: 0.210000\n", + "Taking snapshot at: 0.220000\n", + "Taking snapshot at: 0.230000\n", + "Taking snapshot at: 0.240000\n", + "Taking snapshot at: 0.250000\n", + "step 25, t = 0.250000\n", + "Taking snapshot at: 0.260000\n", + "Taking snapshot at: 0.270000\n", + "Taking snapshot at: 0.280000\n", + "Taking snapshot at: 0.290000\n", + "Taking snapshot at: 0.300000\n", + "step 30, t = 0.300000\n", + "Taking snapshot at: 0.310000\n", + "Taking snapshot at: 0.320000\n", + "Taking snapshot at: 0.330000\n", + "Taking snapshot at: 0.340000\n", + "Taking snapshot at: 0.350000\n", + "step 35, t = 0.350000\n", + "Taking snapshot at: 0.360000\n", + "Taking snapshot at: 0.370000\n", + "Taking snapshot at: 0.380000\n", + "Taking snapshot at: 0.390000\n", + "Taking snapshot at: 0.400000\n", + "step 40, t = 0.400000\n", + "Taking snapshot at: 0.410000\n", + "Taking snapshot at: 0.420000\n", + "Taking snapshot at: 0.430000\n", + "Taking snapshot at: 0.440000\n", + "Taking snapshot at: 0.450000\n", + "step 45, t = 0.450000\n", + "Taking snapshot at: 0.460000\n", + "Taking snapshot at: 0.470000\n", + "Taking snapshot at: 0.480000\n", + "Taking snapshot at: 0.490000\n", + "Taking snapshot at: 0.500000\n", + "step 50, t = 0.500000\n" + ] + } + ], + "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": 19, + "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": 20, + "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": 21, + "id": "8abf542b-e6ef-4637-9f7b-410ee00038b7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating DMD with rdim: 16\n", + "Using 16 basis vectors out of 50.\n", + "Elapsed time for training DMD: 4.050960e+00 second\n", + "Elapsed time for solving FOM: 1.030787e+02 second\n" + ] + } + ], + "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": 22, + "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": 23, + "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..0f3737a --- /dev/null +++ b/examples/dmd/wave_equation.ipynb @@ -0,0 +1,703 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "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": 2, + "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": 3, + "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": 4, + "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": 5, + "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": 6, + "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": 7, + "id": "7ba14cd8-5309-4178-97cb-21daa3b17107", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Options used:\n", + " --mesh ../data/star.mesh\n", + " --refine 2\n", + " --order 2\n", + " --ode_solver 10\n", + " --t_final 0.5\n", + " --time_step 0.01\n", + " --speed 1.0\n", + " --neumann False\n", + " --visualization False\n", + " --visit_datafiles False\n", + " --visualization_steps 5\n", + " --reference \n", + " --energy_fraction 0.9999\n", + " --rdim -1\n", + " --numwindowsamples 9223372036854775807\n" + ] + } + ], + "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": 8, + "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": 9, + "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": 10, + "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": 11, + "id": "a14fcd51-7b98-46ef-a2b4-a532f9de3153", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of temperature unknowns: 1361\n" + ] + } + ], + "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": 12, + "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": 13, + "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": 14, + "id": "342f630c-4684-404a-9e5b-1459a10c77f2", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "All of these memory addresses are different.\n", + "id(uData[0]): 277216868816\n", + "int(u.GetData()): 277217190880\n", + "id(uData): 276715163344\n", + "But uData[*] points to the right memory.\n", + "id(u[0]): 276707586864 =? id(uData[0]): 277216869104\n", + "id(u[1]): 276707586864 =? id(uData[1]): 277216868816\n", + "uData type: \n", + "step 5, t = 0.05\n", + "step 10, t = 0.1\n", + "step 15, t = 0.15\n", + "step 20, t = 0.2\n", + "step 25, t = 0.25\n", + "step 30, t = 0.3\n", + "step 35, t = 0.35\n", + "step 40, t = 0.4\n", + "step 45, t = 0.45\n", + "step 50, t= 0.500000\n", + "Creating DMD with energy fraction: 0.999900 at window index: 0\n", + "Using 15 basis vectors out of 50.\n", + "step 50, t = 0.5\n" + ] + } + ], + "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": 15, + "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": 16, + "id": "c1cd37a0-c6fa-4b86-b87a-f14bc094e649", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Predicting temperature using DMD\n" + ] + } + ], + "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": 17, + "id": "4f2e2fd2-9c9d-4329-a422-f129bde64b37", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Relative error of DMD solution (u) at t_final: 0.500000 is 7.230E-03\n", + "Elapsed time for solving FOM: 2.736976e-01 second\n", + "\n", + "Elapsed time for training DMD: 3.120172e-01 second\n", + "\n", + "Elapsed time for predicting DMD: 6.806588e-02 second\n", + "\n" + ] + } + ], + "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/poisson_global_rom.ipynb b/examples/prom/poisson_global_rom.ipynb index da2f848..ef27c68 100644 --- a/examples/prom/poisson_global_rom.ipynb +++ b/examples/prom/poisson_global_rom.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 26, "id": "eecdf3ba-95dd-456f-93d6-271b41e2fffc", "metadata": {}, "outputs": [], @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 27, "id": "215cf768-c626-4b33-a92b-2e2b9942e926", "metadata": {}, "outputs": [], @@ -411,7 +411,7 @@ { "data": { "text/plain": [ - " >" + " >" ] }, "execution_count": 15, @@ -795,9 +795,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "Elapsed time for assembling FOM: 7.356629e+00 second\n", + "Elapsed time for assembling FOM: 7.358198e+00 second\n", "\n", - "Elapsed time for solving FOM: 1.179325e+00 second\n", + "Elapsed time for solving FOM: 1.184508e+00 second\n", "\n" ] }