diff --git a/examples/prom/elliptic_eigenproblem_global_rom.ipynb b/examples/prom/elliptic_eigenproblem_global_rom.ipynb new file mode 100644 index 0000000..ed3bc03 --- /dev/null +++ b/examples/prom/elliptic_eigenproblem_global_rom.ipynb @@ -0,0 +1,1220 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "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", + "import ctypes\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\n", + "from pylibROM.python_utils.StopWatch import StopWatch" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "sys.path.append(\"../../build\")\n", + "import pylibROM.linalg as libROM\n", + "from pylibROM.mfem import ComputeCtAB" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "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": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from mfem.common.arg_parser import ArgParser\n", + "\n", + "# Create the parser\n", + "parser = ArgParser(description='elliptic_eigenproblem_global_rom')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Default values\n", + "# mesh_file = \"/notebooks/star.mesh\"\n", + "mesh_file = \"\"\n", + "ser_ref_levels = 2\n", + "par_ref_levels = 1\n", + "dirichlet = True\n", + "order = 1\n", + "nev = 4\n", + "seed = 75\n", + "prescribe_init = False\n", + "lobpcg_niter = 200\n", + "lobpcg_tol = 1e-8\n", + "eig_tol = 1e-6\n", + "slu_solver = False\n", + "sp_solver = False\n", + "cpardiso_solver = False\n", + "\n", + "visualization = True\n", + "visit = False\n", + "vis_steps = 5\n", + "fom = False\n", + "offline = False\n", + "merge = False\n", + "online = False\n", + "id = 0\n", + "nsets = 0\n", + "ef = 1.0\n", + "rdim = -1\n", + "verbose_level = 1\n", + "\n", + "precision = 8\n", + "\n", + "problem = 1\n", + "amplitude = 0\n", + "relative_center = 0.0\n", + "L = 1.0 # Scaling Factor\n", + "potential_well_switch = 0\n", + "psedo_time = 0.0\n", + "\n", + "\n", + "# Add arguments\n", + "parser.add_argument('-m', '--mesh', type=str, default=mesh_file, help='Mesh file to use.')\n", + "parser.add_argument('-p', '--problem', type=int, default=problem, help='Problem setup to use.')\n", + "\n", + "\n", + "# Add the Dirichlet argument with a default of True\n", + "parser.add_argument('-dir', '--dirichlet', action='store_false', default=True, help='Use Dirichlet boundary condition (default: True).')\n", + "# Add the Neumann argument with a default of False\n", + "parser.add_argument('-neu', '--neumann', action='store_true', default=False, help='Use Neumann boundary condition (default: False).')\n", + "\n", + "parser.add_argument('-rs', '--refine-serial', type=int, default=ser_ref_levels, help='Number of times to refine the mesh uniformly in serial.')\n", + "parser.add_argument('-rp', '--refine-parallel', type=int, default=par_ref_levels, help='Number of times to refine the mesh uniformly in parallel.')\n", + "parser.add_argument('-o', '--order', type=int, default=order, help='Order (degree) of the finite elements.')\n", + "parser.add_argument('-n', '--num-eigs', type=int, default=nev, help='Number of desired eigenmodes.')\n", + "parser.add_argument('-s', '--seed', type=int, default=seed, help='Random seed used to initialize LOBPCG.')\n", + "parser.add_argument('-a', '--amplitude', type=float, default=amplitude, help='Amplitude of coefficient fields.')\n", + "parser.add_argument('-t', '--pseudo-time', type=float, default=psedo_time, help='Pseudo-time of the motion of the center.')\n", + "parser.add_argument('-c', '--center', type=float, default=relative_center, help='Number of grid elements to which center is shifted.')\n", + "parser.add_argument('-id', '--id', type=int, default=id, help='Parametric id.')\n", + "parser.add_argument('-ns', '--nset', type=int, default=nsets, help='Number of parametric snapshot sets.')\n", + "parser.add_argument('-vis', '--visualization', action='store_true', help='Enable GLVis visualization.')\n", + "parser.add_argument('-visit', '--visit-datafiles', action='store_true', help='Save data files for VisIt visualization.')\n", + "parser.add_argument('-vs', '--visualization-steps', type=int, default=vis_steps, help='Visualize every n-th timestep.')\n", + "parser.add_argument('-fom', '--fom', action='store_true', help='Enable the fom phase.')\n", + "parser.add_argument('-offline', '--offline', action='store_true', help='Enable the offline phase.')\n", + "parser.add_argument('-online', '--online', action='store_true', help='Enable the online phase.')\n", + "parser.add_argument('-merge', '--merge', action='store_true', help='Enable the merge phase.')\n", + "parser.add_argument('-ef', '--energy_fraction', type=float, default=ef, help='Energy fraction.')\n", + "parser.add_argument('-rdim', '--rdim', type=int, default=rdim, help='Reduced dimension.')\n", + "parser.add_argument('-v', '--verbose', type=int, default=verbose_level, help='Set the verbosity level of the LOBPCG solver and preconditioner.')\n", + "parser.add_argument('-fi', '--fom-iter', type=int, default=lobpcg_niter, help='Number of iterations for the LOBPCG solver.')\n", + "parser.add_argument('-tol', '--fom-tol', type=float, default=lobpcg_tol, help='Tolerance for the LOBPCG solver.')\n", + "parser.add_argument('-etol', '--eig-tol', type=float, default=eig_tol, help='Tolerance for eigenvalues to be considered equal.')\n", + "\n", + "parser.add_argument(\"-d\", \"--device\",\n", + " action='store', default='cpu', type=str,\n", + " help=\"Device configuration string, see Device::Configure().\")\n", + "\n", + "parser.add_argument('-slu', '--superlu', default=slu_solver, help='Use the SuperLU Solver.')\n", + "parser.add_argument('-sp', '--strumpack', default=sp_solver, help='Use the STRUMPACK Solver.')\n", + "parser.add_argument('-cpardiso', '--cpardiso', default = cpardiso_solver, help='Use the MKL CPardiso Solver.')\n", + "\n", + "\n", + "\n", + "# Set precision\n", + "import sys\n", + "sys.stdout.write(f\"{precision:.{precision}f}\\n\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run the following commands in sequence:\n", + "\n", + "- First run the **offline** command. For every training parameter (\"a\"), run the offline command once. \n", + "\n", + "- Use different ids for each offline run.\n", + "\n", + "- Once the offline/training dataset is generated run the **merge** command. The merge command builds the tools necessary for the reduced order model, including the reduced basis.\n", + "\n", + "- Finally, run the **online** command to evaluate the reduced order model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Dirichlet BC" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "## run the offline phase 2 times to generate 2 snapshots\n", + "args = parser.parse_args([\"-offline\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\",\"-id\", \"0\", \"-a\", \"0\"])\n", + "# args = parser.parse_args([\"-offline\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\",\"-id\", \"1\", \"-a\", \"1\"])\n", + "\n", + "\n", + "# args = parser.parse_args([\"-merge\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-ns\", \"2\"])\n", + "\n", + "# args = parser.parse_args([\"-fom\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-a\", \"0.5\"])\n", + "\n", + "\n", + "# args = parser.parse_args([\"-online\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-a\", \"0.5\", \"-ef\", \"1.0\"])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Neumann BC" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "## run the offline phase 2 times to generate 2 snapshots\n", + "# args = parser.parse_args([\"-offline\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\",\"-id\", \"0\", \"-a\", \"0\",\"-neu\"])\n", + "# args = parser.parse_args([\"-offline\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\",\"-id\", \"1\", \"-a\", \"1\",\"-neu\"])\n", + "\n", + "\n", + "# args = parser.parse_args([\"-merge\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-ns\", \"2\", \"-neu\"])\n", + "\n", + "# args = parser.parse_args([\"-fom\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-a\", \"0.5\",\"-neu\"])\n", + "\n", + "\n", + "# args = parser.parse_args([\"-online\", \"-p\", \"2\", \"-rs\", \"2\", \"-n\", \"4\", \"-a\", \"0.5\", \"-ef\", \"1.0\", \"-neu\"])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Parse the arguments\n", + "# args = parser.parse_args([\"-p 2 -rs 2 -n 4 -id 0 -a 0\"])\n", + "if (myid == 0): \n", + " parser.print_options(args)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "mesh_file = args.mesh\n", + "ser_ref_levels = args.refine_serial\n", + "par_ref_levels = args.refine_parallel\n", + "dirichlet = args.dirichlet\n", + "neumann = args.neumann\n", + "order = args.order\n", + "nev = args.num_eigs\n", + "seed = args.seed\n", + "slu_solver = args.superlu\n", + "sp_solver = args.strumpack\n", + "cpardiso_solver = args.cpardiso\n", + "\n", + "visualization = args.visualization\n", + "visit = args.visit_datafiles\n", + "vis_steps = args.visualization_steps\n", + "fom = args.fom\n", + "offline = args.offline\n", + "merge = args.merge\n", + "online = args.online\n", + "id = args.id\n", + "nsets = args.nset\n", + "ef = args.energy_fraction\n", + "rdim = args.rdim\n", + "verbose_level = args.verbose\n", + "problem = args.problem\n", + "amplitude = args.amplitude\n", + "relative_center = args.center\n", + "\n", + "\n", + "\n", + "# Check for conflicting solvers\n", + "if args.superlu and args.strumpack:\n", + " print(\"WARNING: Both SuperLU and STRUMPACK have been selected, please choose either one.\")\n", + " print(\" Defaulting to SuperLU.\")\n", + " args.strumpack = False\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "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, + "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(args.device)\n", + "if (myid == 0):\n", + " device.Print()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "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", + "\n", + "# Create a mesh\n", + "if args.mesh == \"\":\n", + " mesh = mfem.Mesh(mfem.Mesh.MakeCartesian2D(2, 2, mfem.Element.QUADRILATERAL))\n", + "else:\n", + " mesh = mfem.Mesh(args.mesh, 1, 1)\n", + "\n", + "dim = mesh.Dimension()" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the bounding box\n", + "bb_min, bb_max = mesh.GetBoundingBox(max(args.order, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# Vh = mfem.Vector()\n", + "# Vk z= mfem.Vector()\n", + "\n", + "h_min = mfem.doublep(); \n", + "h_max = mfem.doublep(); \n", + "kappa_min = mfem.doublep(); \n", + "kappa_max = mfem.doublep()\n", + "\n", + "mesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max)\n", + "\n", + "h_min = h_min.value();\n", + "h_max = h_max.value(); \n", + "kappa_min = kappa_min.value(); \n", + "kappa_max = kappa_max.value()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "# Update h_max\n", + "h_max *= np.power(0.5, ser_ref_levels + par_ref_levels)\n", + "\n", + "# Calculate the center\n", + "center = mfem.Vector(dim)\n", + "for i in range(dim):\n", + " center[i] = h_max * relative_center + 0.5 * (bb_min[i] + bb_max[i])" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "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", + "for l in range(ser_ref_levels):\n", + " mesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "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", + "for l in range(par_ref_levels):\n", + " pmesh.UniformRefinement()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "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", + "\n", + "size = fespace.GlobalTrueVSize()\n", + "if (myid == 0):\n", + " print(\"Number of finite element unknowns: %d\" % size)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# Initiate ROM related variables\n", + "max_num_snapshots = 100\n", + "update_right_SV = False\n", + "isIncremental = False\n", + "data_dir = os.getcwd()+\"/data/\"\n", + "baseName = \"elliptic_eigenproblem_\"\n", + "basis_filename = data_dir+baseName+\"basis\"\n", + "solveTimer, assembleTimer, mergeTimer = StopWatch(), StopWatch(), StopWatch()" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "# 10. Set BasisGenerator if offline\n", + "if (offline):\n", + " options = libROM.Options(fespace.GetTrueVSize(), nev, nev,\n", + " update_right_SV)\n", + " snapshot_basename = data_dir+baseName + \"par\" + f\"{id}\"\n", + " generator = libROM.BasisGenerator(options, isIncremental, snapshot_basename)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# 11. The merge phase\n", + "if (merge):\n", + " mergeTimer.Start()\n", + " options = libROM.Options(fespace.GetTrueVSize(), max_num_snapshots, nev,\n", + " update_right_SV)\n", + " generator = libROM.BasisGenerator(options, isIncremental, basis_filename)\n", + " for paramID in range(nsets):\n", + " snapshot_filename = \"%s%d_snapshot\" % (baseName+\"par\", paramID)\n", + " generator.loadSamples(data_dir+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()" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "class Conductivity(mfem.PyCoefficient):\n", + " \"\"\"\n", + " Defines the conductivity as a function of position.\n", + " \n", + " This class extends mfem.PyCoefficient and specifies the conductivity based on the problem and position vector x.\n", + " \"\"\"\n", + "\n", + " def __init__(self, problem, amplitude, center, bb_min, bb_max):\n", + " \"\"\"\n", + " Initializes the Conductivity object.\n", + " \n", + " Parameters:\n", + " - problem (int): The problem setup identifier.\n", + " - amplitude (float): The amplitude value for the conductivity calculation.\n", + " - center (array): The center position as a NumPy array.\n", + " - bb_min (array): The minimum bounding box coordinates as a NumPy array.\n", + " - bb_max (array): The maximum bounding box coordinates as a NumPy array.\n", + " \"\"\"\n", + " super().__init__()\n", + " self.problem = problem\n", + " self.amplitude = amplitude\n", + " self.center = center\n", + " self.bb_min = bb_min\n", + " self.bb_max = bb_max\n", + "\n", + " def EvalValue(self, x):\n", + " \"\"\"\n", + " Evaluates the conductivity at a given position x.\n", + " \n", + " Parameters:\n", + " - x (array): Position at which to evaluate the conductivity, given as a NumPy array.\n", + " \n", + " Returns:\n", + " - float: The evaluated conductivity value at position x.\n", + " \"\"\"\n", + " cx = 1.0\n", + " if self.problem == 1:\n", + " return 1.0\n", + " elif self.problem == 2:\n", + " cx += self.amplitude\n", + " for i in range(len(x)):\n", + " if 8 * abs(x[i] - self.center[i]) > (self.bb_max[i] - self.bb_min[i]):\n", + " cx = 1.0\n", + " return cx\n", + " elif self.problem in [3, 4]:\n", + " return 1.0\n", + " return 0.0" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "class Potential(mfem.PyCoefficient):\n", + " \"\"\"\n", + " Defines the potential as a function of position.\n", + " \n", + " This class extends mfem.PyCoefficient and specifies the potential based on the problem and position vector x.\n", + " \"\"\"\n", + "\n", + " def __init__(self, problem, amplitude, center, h_max):\n", + " \"\"\"\n", + " Initializes the Potential object.\n", + " \n", + " Parameters:\n", + " - problem (int): The problem setup identifier.\n", + " - amplitude (float): The amplitude value for the potential calculation.\n", + " - center (array): The center position as a NumPy array.\n", + " - h_max (float): The maximum mesh size.\n", + " \"\"\"\n", + " super().__init__()\n", + " self.problem = problem\n", + " self.amplitude = amplitude\n", + " self.center = center\n", + " self.h_max = h_max\n", + "\n", + "\n", + " def EvalValue(self, x):\n", + " \"\"\"\n", + " Evaluates the potential at a given position x.\n", + " \n", + " Parameters:\n", + " - x (array): Position at which to evaluate the potential, given as a NumPy array.\n", + " \n", + " Returns:\n", + " - float: The evaluated potential value at position x.\n", + " \"\"\"\n", + " radius = 5.0 * self.h_max\n", + " # x_np = x.GetDataArray() # Convert MFEM Vector to NumPy array\n", + " # print(type(self.center))\n", + " d_sq = np.sqrt(np.sum((x-self.center.GetDataArray())**2))\n", + " if self.problem in [1, 2]:\n", + " return 0.0\n", + " elif self.problem == 3:\n", + " return self.amplitude * np.exp(-d_sq / np.power(radius, 2.0))\n", + " elif self.problem == 4:\n", + " return self.amplitude * d_sq\n", + " return 0.0" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "assembleTimer.Start()\n", + "one = mfem.ConstantCoefficient(1.0)\n", + "# Create a parallel grid function\n", + "\n", + "kappa_0 = Conductivity(problem, amplitude, center, bb_min, bb_max)\n", + "v_0 = Potential(problem, amplitude, center, h_max)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "c_gf = mfem.ParGridFunction(fespace)\n", + "c_gf.ProjectCoefficient(kappa_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "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", + "\n", + "p_gf = mfem.ParGridFunction(fespace)\n", + "p_gf.ProjectCoefficient(v_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "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", + "\n", + "if (pmesh.bdr_attributes.Size() > 0):\n", + " ess_bdr = mfem.intArray(pmesh.bdr_attributes.Max())\n", + " ess_bdr.Assign(1 if dirichlet else 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "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", + "\n", + "a = mfem.ParBilinearForm(fespace)\n", + "a.AddDomainIntegrator(mfem.DiffusionIntegrator(kappa_0))\n", + "a.AddDomainIntegrator(mfem.MassIntegrator(v_0))\n", + "a.Assemble()\n", + "a.EliminateEssentialBCDiag(ess_bdr, 1.0)\n", + "a.Finalize()" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "one = mfem.ConstantCoefficient(1.0)\n", + "m = mfem.ParBilinearForm(fespace)\n", + "m.AddDomainIntegrator(mfem.MassIntegrator(one))\n", + "m.Assemble()\n", + "m.EliminateEssentialBCDiag(ess_bdr, 1/np.finfo(np.float64).min)\n", + "m.Finalize()" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "A = mfem.HypreParMatrix()\n", + "M = mfem.HypreParMatrix()\n", + "\n", + "# Parallel assemble the bilinear forms\n", + "A = a.ParallelAssemble()\n", + "M = m.ParallelAssemble()" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "# Assuming `A` is an instance of HypreParMatrix and is already defined.\n", + "# Arow = mfem.Operator()\n", + "Arow = None\n", + "\n", + "# Check if SuperLU solver is available and selected\n", + "try:\n", + " from mfem.common import SuperLURowLocMatrix\n", + " if slu_solver: # Assuming slu_solver is a boolean variable\n", + " Arow = mfem.SuperLURowLocMatrix(A)\n", + "except ImportError:\n", + " pass\n", + "\n", + "# Check if STRUMPACK solver is available and selected\n", + "try:\n", + " from mfem.common import STRUMPACKRowLocMatrix\n", + " if sp_solver: # Assuming sp_solver is a boolean variable\n", + " Arow = mfem.STRUMPACKRowLocMatrix(A)\n", + "except ImportError:\n", + " pass\n", + "\n", + "# Stop the assembly timer\n", + "assembleTimer.Stop() # You need to have the `assembleTimer` defined and started somewhere in your code.\n", + "\n", + "# Delete the bilinear forms\n", + "del a\n", + "del m\n", + "\n", + "# Create a parallel grid function\n", + "# x = mfem.ParGridFunction(fespace)\n", + "eigenfunction_i = mfem.ParGridFunction(fespace)\n", + "eigenvector_i = mfem.Vector(size)\n", + "# Create a HypreLOBPCG solver\n", + "lobpcg = mfem.HypreLOBPCG(comm)\n", + "\n", + "# Create arrays for eigenvalues and eigenvectors\n", + "eigenvalues = mfem.doubleArray()\n", + "# eigenvectors = mfem.DenseMatrix()\n", + "eigenvectors = []" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "prescribe_init = False" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if fom or offline:\n", + " # Define and configure the LOBPCG eigensolver and the BoomerAMG preconditioner for A\n", + " precond = None\n", + "\n", + " if not slu_solver and not sp_solver and not cpardiso_solver:\n", + " amg = mfem.HypreBoomerAMG(A)\n", + " amg.SetPrintLevel(verbose_level)\n", + " precond = amg\n", + " else:\n", + " if slu_solver:\n", + " try:\n", + " from mfem import SuperLUSolver\n", + " superlu = mfem.SuperLUSolver(MPI.COMM_WORLD)\n", + " superlu.SetPrintStatistics(verbose_level > 0)\n", + " superlu.SetSymmetricPattern(True)\n", + " superlu.SetColumnPermutation(\"PARMETIS\")\n", + " superlu.SetOperator(Arow)\n", + " precond = superlu\n", + " except ImportError:\n", + " pass\n", + "\n", + " if sp_solver:\n", + " try:\n", + " from mfem import STRUMPACKSolver\n", + " strumpack = STRUMPACKSolver(MPI.COMM_WORLD)\n", + " strumpack.SetPrintFactorStatistics(True)\n", + " strumpack.SetPrintSolveStatistics(verbose_level > 0)\n", + " strumpack.SetKrylovSolver(\"DIRECT\")\n", + " strumpack.SetReorderingStrategy(\"METIS\")\n", + " strumpack.SetMatching(\"NONE\")\n", + " strumpack.SetCompression(\"NONE\")\n", + " strumpack.SetOperator(Arow)\n", + " strumpack.SetFromCommandLine()\n", + " precond = strumpack\n", + " except ImportError:\n", + " pass\n", + "\n", + " if cpardiso_solver:\n", + " try:\n", + " from mfem import CPardisoSolver\n", + " cpardiso = CPardisoSolver(A.GetComm())\n", + " cpardiso.SetMatrixType(\"REAL_STRUCTURE_SYMMETRIC\")\n", + " cpardiso.SetPrintLevel(verbose_level)\n", + " cpardiso.SetOperator(A)\n", + " precond = cpardiso\n", + " except ImportError:\n", + " pass\n", + "\n", + " # Create and configure the HypreLOBPCG solver\n", + " lobpcg = mfem.HypreLOBPCG(MPI.COMM_WORLD)\n", + " lobpcg.SetNumModes(nev)\n", + " lobpcg.SetRandomSeed(seed)\n", + " lobpcg.SetPreconditioner(precond)\n", + " lobpcg.SetMaxIter(lobpcg_niter)\n", + " lobpcg.SetTol(1e-8)\n", + " lobpcg.SetPrecondUsageMode(1)\n", + " lobpcg.SetPrintLevel(verbose_level)\n", + " lobpcg.SetMassMatrix(M)\n", + " lobpcg.SetOperator(A)\n", + "\n", + "\n", + "\n", + " if prescribe_init and (fom or (offline and id > 0)):\n", + " snapshot_vecs = [np.loadtxt(data_dir+f\"{baseName}ref_snapshot_{i}.{myid}\") for i in range(nev)]\n", + " lobpcg.SetInitialVectors(np.array(snapshot_vecs))\n", + " if myid == 0:\n", + " print(\"LOBPCG initial vectors set\")\n", + "\n", + "\n", + " solveTimer.Start()\n", + " lobpcg.Solve()\n", + " solveTimer.Stop()\n", + " \n", + " # Extract eigenvalues\n", + " eigenvalues = mfem.doubleArray(nev)\n", + " lobpcg.GetEigenvalues(eigenvalues)\n", + "\n", + " for i in range(nev):\n", + " if myid == 0:\n", + " print(f\"Eigenvalue {i}: {eigenvalues[i]}\")\n", + " if offline:\n", + " eigenfunction_i.Assign(lobpcg.GetEigenvector(i))\n", + " eigenfunction_i /= sqrt(mfem.InnerProduct(eigenfunction_i, eigenfunction_i))\n", + " # generator.takeSample(eigenfunction_i.GetDataArray(),0,0)\n", + " generator.takeSample(eigenfunction_i.GetDataArray())\n", + "\n", + " if prescribe_init and id == 0:\n", + "\n", + " snapshot_filename = data_dir+f\"{baseName}ref_snapshot_{i}\"\n", + " snapshot_vec = lobpcg.GetEigenvector(i)\n", + " np.savetxt(snapshot_filename, snapshot_vec) # Assuming snapshot_vec is a NumPy array or similar\n", + " if myid == 0:\n", + " print(f\"Saved {snapshot_filename}\")\n", + "\n", + "\n", + " if offline:\n", + " print(\"here\")\n", + " generator.writeSnapshot()\n", + " del generator\n", + " del options\n", + "\n", + " # del precond" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [], + "source": [ + "if online:\n", + " # 16. Read the reduced basis\n", + " assembleTimer.Start()\n", + " reader = libROM.BasisReader(basis_filename)\n", + "\n", + " if rdim != -1:\n", + " spatialbasis = reader.getSpatialBasis(rdim)\n", + " else:\n", + " spatialbasis = reader.getSpatialBasis(ef)\n", + "\n", + " numRowRB = spatialbasis.numRows()\n", + " numColumnRB = spatialbasis.numColumns()\n", + " if myid == 0:\n", + " print(f\"Spatial basis dimension is {numRowRB} x {numColumnRB}\")\n", + "\n", + " # 17. Form ROM operator\n", + " ReducedA = libROM.Matrix()\n", + " ComputeCtAB(A, spatialbasis, spatialbasis, ReducedA)\n", + " A_mat = mfem.DenseMatrix(ReducedA.numRows(), ReducedA.numColumns())\n", + "\n", + " rows, cols = ReducedA.numRows(), ReducedA.numColumns() # Get matrix dimensions\n", + " ReducedA_data = ReducedA.getData().reshape((rows, cols)) # Reshape the flat data into a matrix form\n", + "\n", + " # Create a DenseMatrix from the reshaped data\n", + " A_mat.Set(1.0, mfem.DenseMatrix(ReducedA_data))\n", + "\n", + " ReducedM = libROM.Matrix()\n", + " ComputeCtAB(M, spatialbasis, spatialbasis, ReducedM)\n", + " M_mat = mfem.DenseMatrix(ReducedM.numRows(), ReducedM.numColumns())\n", + " \n", + " rows_M, cols_M = ReducedM.numRows(), ReducedM.numColumns() # Get matrix dimensions\n", + " ReducedM_data = ReducedM.getData().reshape((rows_M, cols_M)) # Reshape the flat data into a matrix form\n", + "\n", + "\n", + " M_mat.Set(1.0, mfem.DenseMatrix(ReducedM_data))\n", + " ### I am here \n", + "\n", + " assembleTimer.Stop()\n", + "\n", + " eigenvalues_rom = mfem.Vector() # Initialize with correct size\n", + " reduced_eigenvectors = mfem.DenseMatrix() # Initialize with correct size\n", + "\n", + "\n", + " # 18. Solve ROM\n", + " solveTimer.Start()\n", + " A_mat.Eigenvalues(M_mat, eigenvalues_rom, reduced_eigenvectors)\n", + " solveTimer.Stop()\n", + "\n", + " if myid == 0:\n", + "\n", + " size = eigenvalues_rom.Size()\n", + " data_ptr = eigenvalues_rom.GetDataArray()\n", + " eigenvalues = mfem.doubleArray(size)\n", + "\n", + " for i in range(size):\n", + " eigenvalues[i] = data_ptr[i]\n", + "\n", + " mode_rom = mfem.Vector(numRowRB)\n", + " modes_rom = []\n", + "\n", + " for i in range(min(eigenvalues_rom.Size(), nev)):\n", + " print(f\"Eigenvalue {i}: = {eigenvalues[i]}\")\n", + "\n", + " # tmp = mfem.DenseMatrix(eigenvectors)\n", + " # eigenvectors = mfem.DenseMatrix(nev, numRowRB)\n", + "\n", + " for j in range(min(eigenvalues_rom.Size(), nev)):\n", + " reduced_eigenvector_j = mfem.Vector()\n", + " reduced_eigenvectors.GetColumn(j, reduced_eigenvector_j)\n", + " reduced_eigenvector_j_carom = libROM.Vector(reduced_eigenvector_j.GetDataArray(), False, False)\n", + " eigenvector_j_carom = spatialbasis.mult(reduced_eigenvector_j_carom)\n", + "\n", + " mode_rom.Assign(eigenvector_j_carom.getData())\n", + "\n", + " mode_rom /= sqrt(mfem.InnerProduct(mode_rom,mode_rom))\n", + " modes_rom.append(mfem.Vector(mode_rom))\n", + " \n", + " del eigenvector_j_carom\n", + "\n", + " sol_eigenvalue_name_fom = data_dir+f\"sol_eigenvalues_fom.{myid:06d}\"\n", + " eigenvalues_fom = mfem.Vector(nev)\n", + " # Pass file path directly\n", + " eigenvalues_fom.Load(sol_eigenvalue_name_fom,nev)\n", + "\n", + "\n", + " # Assuming 'eigenvalues_rom' has been computed already (e.g., another vector of same size)\n", + " diff_eigenvalues = mfem.Vector(nev)\n", + "\n", + " for i in range(nev):\n", + " diff_eigenvalues[i] = eigenvalues_fom[i] - eigenvalues_rom[i]\n", + " if myid == 0:\n", + " print(f\"FOM solution for eigenvalue {i} = {eigenvalues_fom[i]}\")\n", + " print(f\"ROM solution for eigenvalue {i} = {eigenvalues_rom[i]}\")\n", + " print(f\"Absolute error of ROM solution for eigenvalue {i} = {abs(diff_eigenvalues[i])}\")\n", + " print(f\"Relative error of ROM solution for eigenvalue {i} = {abs(diff_eigenvalues[i]) / abs(eigenvalues_fom[i])}\")\n", + "\n", + "\n", + " # Load and calculate errors for eigenvectors\n", + " eigenvectors = []\n", + " for i in range(nev):\n", + " mode_name_fom = f\"eigenfunction_fom_{i:02d}.{myid:06d}\"\n", + " \n", + " # Load FOM eigenvector from file\n", + " mode_fom = mfem.Vector(numRowRB)\n", + " mode_fom.Load(data_dir+mode_name_fom, numRowRB)\n", + "\n", + " eigenvector_i = mfem.Vector(numRowRB) # Initialize eigenvector_i to zero\n", + " eigenvector_i.Assign(0.0)\n", + " \n", + " for j in range(nev):\n", + " if myid == 0 and verbose_level > 0:\n", + " print(f\"correlation_matrix({j+1},{i+1}) = {mfem.InnerProduct(mode_fom, modes_rom[j])};\")\n", + " \n", + " # If eigenvalues are nearly identical, add the corresponding ROM mode\n", + " if abs(eigenvalues_fom[j] - eigenvalues_fom[i]) < 1e-6:\n", + " eigenvector_i.Add(mfem.InnerProduct(mode_fom, modes_rom[j]), modes_rom[j])\n", + "\n", + " # Store the computed eigenvector\n", + " eigenvectors.append(eigenvector_i)\n", + " \n", + " # Compute the difference between FOM and ROM eigenvectors\n", + " mode_fom.Add(-1.0, eigenvector_i)\n", + " diff_norm = sqrt(mfem.InnerProduct(mode_fom, mode_fom))\n", + " \n", + " if myid == 0:\n", + " print(f\"Relative l2 error of ROM eigenvector {i} = {diff_norm}\")\n", + "\n", + "\n", + "\n", + " del spatialbasis\n", + " del A_mat\n", + " del M_mat" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the filename format helper\n", + "def get_filename(prefix, myid, digits=6):\n", + " return f\"{prefix}{myid:0{digits}}\"\n", + "\n", + "\n", + "## Print Mesh\n", + "mesh_name = get_filename(data_dir+\"elliptic_eigenproblem-mesh.\", myid)\n", + "pmesh.Print(mesh_name, precision)\n", + "\n", + "\n", + "## Print Eigenvalue\n", + "if fom or offline:\n", + " sol_eigenvalue_name = get_filename(data_dir+\"sol_eigenvalues_fom.\", myid)\n", + "\n", + "if online:\n", + " sol_eigenvalue_name = get_filename(data_dir+\"sol_eigenvalues.\", myid)\n", + "\n", + "# pmesh.Print(sol_eigenvalue_name, precision)" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [], + "source": [ + "# 19. Save the refined mesh and the modes in parallel.\n", + "sign_eigenvectors = mfem.Vector(nev) ## ?? no idea what is this!\n", + "\n", + "mode_prefix = data_dir+\"eigenfunction_\"\n", + "\n", + "if fom or offline:\n", + " mode_prefix +=\"fom_\"\n", + "\n", + "elif online:\n", + " mode_prefix += \"rom_\"\n", + "\n", + "\n", + "# Save the eigenvalues\n", + "with open(sol_eigenvalue_name, 'w') as eigvals:\n", + " for val in eigenvalues:\n", + " eigvals.write(f\"{val}\\n\")\n", + "\n", + "\n", + "# Save the eigenfunctions\n", + "for i in range(nev):\n", + "\n", + " if fom or offline:\n", + " eigenfunction_i = lobpcg.GetEigenvector(i)\n", + " eigenfunction_i /= sqrt(mfem.InnerProduct(eigenfunction_i, eigenfunction_i))\n", + " else:\n", + " eigenfunction_i = eigenvectors[i]\n", + "\n", + " # mode_ref = mfem.Vector(eigenvector_i.Size())\n", + " mode_name = f\"{mode_prefix}{i:02}.{myid:06}\"\n", + "\n", + "\n", + " with open(mode_name, 'w') as mode_ofs:\n", + " for val in eigenfunction_i:\n", + " mode_ofs.write(f\"{val}\\n\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Visualization" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [], + "source": [ + "visit_dc = None\n", + "if visit:\n", + " if offline:\n", + " visit_dc = mfem.VisItDataCollection(f\"{baseName}offline_par{id}\", pmesh)\n", + " elif fom:\n", + " visit_dc = mfem.VisItDataCollection(f\"{baseName}fom\", pmesh)\n", + " elif online:\n", + " visit_dc = mfem.VisItDataCollection(f\"{baseName}rom\", pmesh)\n", + " \n", + " visit_dc.RegisterField(\"Conductivity\", c_gf)\n", + " visit_dc.RegisterField(\"Potential\", p_gf)\n", + " visit_eigenvectors = []\n", + "\n", + " for i in range(min(nev, eigenvalues.Size())):\n", + " if fom or offline:\n", + " eigenvector_i = lobpcg.GetEigenvector(i)\n", + " else:\n", + " eigenvector_i = mfem.Vector()\n", + " eigenvectors.GetRow(i, eigenvector_i)\n", + "\n", + " eigenvector_i *= sign_eigenvectors[i] / np.sqrt(mfem.InnerProduct(eigenvector_i, eigenvector_i))\n", + "\n", + " x.Assign(eigenvector_i)\n", + " eigenvector_gf = mfem.ParGridFunction(x)\n", + " visit_eigenvectors.append(mfem.Vector(eigenvector_gf))\n", + " visit_dc.RegisterField(f\"Eigenmode_{i}\", eigenvector_gf)\n", + "\n", + " visit_dc.SetCycle(0)\n", + " visit_dc.SetTime(0.0)\n", + " visit_dc.Save()\n", + "\n", + " visit_eigenvectors.clear()\n", + "\n", + "\n", + "sout = None\n", + "if visualization:\n", + " vishost = \"localhost\"\n", + " visport = 19916\n", + " sout = mfem.socketstream(vishost, visport)\n", + " sout << f\"parallel {num_procs} {myid}\\n\"\n", + " \n", + " good = sout.good()\n", + " all_good = MPI.COMM_WORLD.allreduce(good, op=MPI.MIN)\n", + " \n", + " if not all_good:\n", + " sout.close()\n", + " visualization = False\n", + " if myid == 0:\n", + " print(f\"Unable to connect to GLVis server at {vishost}:{visport}\")\n", + " print(\"GLVis visualization disabled.\")\n", + " else:\n", + " for i in range(min(nev, eigenvalues.Size())):\n", + " if myid == 0:\n", + " print(f\"Eigenmode {i + 1}/{nev}, Lambda = {eigenvalues[i]}\")\n", + "\n", + " if fom or offline:\n", + " eigenvector_i = lobpcg.GetEigenvector(i)\n", + " else:\n", + " eigenvector_i = mfem.Vector()\n", + " eigenvectors.GetRow(i, eigenvector_i)\n", + " \n", + " eigenvector_i *= sign_eigenvectors[i] / np.sqrt(mfem.InnerProduct(eigenvector_i, eigenvector_i))\n", + " x.Assign(eigenvector_i)\n", + "\n", + " sout << f\"parallel {num_procs} {myid}\\nsolution\\n\" << pmesh << x << \"\\n\"\n", + " sout << f\"window_title 'Eigenmode {i + 1}/{nev}, Lambda = {eigenvalues[i]}'\\n\"\n", + "\n", + " if myid == 0:\n", + " c = input(\"press (q)uit or (c)ontinue --> \").strip()\n", + " else:\n", + " c = None\n", + " c = MPI.COMM_WORLD.bcast(c, root=0)\n", + "\n", + " if c != 'c':\n", + " break\n", + " sout.close()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 20. Print timing info\n", + "if myid == 0:\n", + " if fom or offline:\n", + " print(f\"Elapsed time for assembling FOM: {assembleTimer.duration:e} seconds\")\n", + " print(f\"Elapsed time for solving FOM: {solveTimer.duration:e} seconds\")\n", + " if online:\n", + " print(f\"Elapsed time for assembling ROM: {assembleTimer.duration:e} seconds\")\n", + " print(f\"Elapsed time for solving ROM: {solveTimer.duration:e} seconds\")\n", + "\n", + "# 21. Free the used memory.\n", + "if fom or offline:\n", + " del lobpcg\n", + "\n", + "del M\n", + "del A\n", + "\n", + "del fespace\n", + "if order > 0:\n", + " del fec\n", + "del pmesh\n", + "\n", + "# Free Arow if it was used\n", + "if 'Arow' in locals():\n", + " del Arow\n", + "\n", + "MPI.Finalize()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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": 2 +}