Global ROM for linear elasticity

include/linelast_solver.hpp

@@ -75,7 +75,7 @@ class LinElastSolver : public MultiBlockSolver
    virtual void AssembleOperator();
    // For bilinear case.
    // system-specific.
-   virtual void AssembleInterfaceMatrixes();
+   virtual void AssembleInterfaceMatrices();
 
    virtual bool Solve();
 

src/linelast_solver.cpp

@@ -131,8 +131,8 @@ void LinElastSolver::InitVariables()
       // Does this make any difference?
       us[m]->SetTrueVector();
    }
-   // if (use_rom)  //Off for now
-   //   MultiBlockSolver::InitROMHandler();
+   if (use_rom)
+     MultiBlockSolver::InitROMHandler();
 }
 
 void LinElastSolver::BuildOperators()
@@ -245,7 +245,7 @@ void LinElastSolver::AssembleOperator()
          }
       }
    }
-   AssembleInterfaceMatrixes();
+   AssembleInterfaceMatrices();
    for (int m = 0; m < numSub; m++)
       as[m]->Finalize();
 
@@ -281,10 +281,10 @@ void LinElastSolver::AssembleOperator()
    }
 }
 
-void LinElastSolver::AssembleInterfaceMatrixes()
+void LinElastSolver::AssembleInterfaceMatrices()
 {
    assert(a_itf);
-   a_itf->AssembleInterfaceMatrixes(mats);
+   a_itf->AssembleInterfaceMatrices(mats);
 }
 
 bool LinElastSolver::Solve()
@@ -467,11 +467,18 @@ void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
    }
 }
 
+void LinElastSolver::ProjectOperatorOnReducedBasis()
+{ 
+   Array2D<Operator *> tmp(mats.NumRows(), mats.NumCols());
+   for (int i = 0; i < tmp.NumRows(); i++)
+      for (int j = 0; j < tmp.NumCols(); j++)
+         tmp(i, j) = mats(i, j);
+
+   rom_handler->ProjectOperatorOnReducedBasis(tmp);
+}
+
 // Component-wise assembly
 void LinElastSolver::BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) { "LinElastSolver::BuildCompROMElement is not implemented yet!\n"; }
 void LinElastSolver::BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp) { "LinElastSolver::BuildBdrROMElement is not implemented yet!\n"; }
 void LinElastSolver::BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp) { "LinElastSolver::BuildInterfaceROMElement is not implemented yet!\n"; }
-
-void LinElastSolver::ProjectOperatorOnReducedBasis() { "LinElastSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n"; }
-
 void LinElastSolver::SanityCheckOnCoeffs() { "LinElastSolver::SanityCheckOnCoeffs is not implemented yet!\n"; }

src/main_workflow.cpp

@@ -337,7 +337,6 @@ void BuildROM(MPI_Comm comm)
    // TODO: there are skippable operations depending on rom/fom mode.
    test->BuildOperators();
    test->SetupBCOperators();
-
    test->LoadReducedBasis();
 
    TopologyHandlerMode topol_mode = test->GetTopologyMode();
@@ -346,7 +345,6 @@ void BuildROM(MPI_Comm comm)
    // NOTE(kevin): global operator required only for global rom operator.
    if (save_operator == ROMBuildingLevel::GLOBAL)
       test->Assemble();
-
    switch (save_operator)
    {
       case ROMBuildingLevel::COMPONENT:

src/rom_handler.cpp

@@ -222,7 +222,6 @@ void ROMHandlerBase::LoadReducedBasis()
       else
          basis_name += GetBasisTagForComponent(k, train_mode, topol_handler);
       basis_reader = new CAROM::BasisReader(basis_name);
-
       carom_ref_basis[k] = basis_reader->getSpatialBasis(0.0, num_ref_basis[k]);
       numRowRB = carom_ref_basis[k]->numRows();
       numColumnRB = carom_ref_basis[k]->numColumns();

test/CMakeLists.txt

@@ -41,6 +41,9 @@ file(COPY inputs/steady_ns.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/in
 file(COPY meshes/test.1x1.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)
 file(COPY meshes/test.global.h5 DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)
 
+file(COPY inputs/linelast.base.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
+file(COPY meshes/beam-tri.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes)
+
 add_executable(test_hdf5 test_hdf5_utils.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
 add_executable(test_topol test_topol.cpp $<TARGET_OBJECTS:scaleupROMObj>)

test/dg_integ_mms.cpp

@@ -116,7 +116,7 @@ TEST(InterfaceDGElasticityIntegrator, Test_Quad)
    for (int i = 0; i < topol_data.numSub; i++)
       for (int j = 0; j < topol_data.numSub; j++)
          mats(i, j) = new SparseMatrix(fes[i]->GetTrueVSize(), fes[j]->GetTrueVSize());
-   a_itf.AssembleInterfaceMatrixes(mats);
+   a_itf.AssembleInterfaceMatrices(mats);
    for (int i = 0; i < topol_data.numSub; i++)
       for (int j = 0; j < topol_data.numSub; j++)
          mats(i, j)->Finalize();

test/inputs/linelast.base.yml

@@ -0,0 +1,64 @@
+main:
+  mode: single_run
+  use_rom: true
+  solver: linelast
+
+mesh:
+  filename: meshes/beam-tri.mesh
+  uniform_refinement: 1
+
+domain-decomposition:
+  type: interior_penalty
+
+discretization:
+  order: 1
+  full-discrete-galerkin: true
+
+visualization:
+  enabled: false
+  unified_paraview: true
+  file_path:
+    prefix: sample_gen_output
+
+parameterized_problem:
+  name: linelast_disp
+
+solver:
+  direct_solve: true
+
+single_run:
+  linelast_disp:
+    rdisp_f: 1.0
+
+sample_generation:
+  maximum_number_of_snapshots: 400
+  component_sampling: false
+  file_path:
+    prefix: "linelast_disp"
+  parameters:
+    - key: single_run/linelast_disp/rdisp_f
+      type: double
+      sample_size: 2
+      minimum: 0.9
+      maximum: 1.1
+
+basis:
+  prefix: "linelast_disp"
+  number_of_basis: 3
+  tags:
+    - name: comp0
+  svd:
+    save_spectrum: true
+    update_right_sv: false
+  visualization:
+      enabled: false
+
+model_reduction:
+  rom_handler_type: base
+  # individual/universal
+  subdomain_training: individual
+  save_operator:
+    level: global
+    prefix: "proj_inv"
+  compare_solution:
+    enabled: true

test/meshes/beam-tri.mesh

@@ -0,0 +1,77 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+16
+1 2 10 0 1
+1 2 0 10 9
+1 2 11 1 2
+1 2 1 11 10
+1 2 12 2 3
+1 2 2 12 11
+1 2 13 3 4
+1 2 3 13 12
+2 2 14 4 5
+2 2 4 14 13
+2 2 15 5 6
+2 2 5 15 14
+2 2 16 6 7
+2 2 6 16 15
+2 2 17 7 8
+2 2 7 17 16
+
+boundary
+18
+3 1 0 1
+3 1 1 2
+3 1 2 3
+3 1 3 4
+3 1 4 5
+3 1 5 6
+3 1 6 7
+3 1 7 8
+3 1 10 9
+3 1 11 10
+3 1 12 11
+3 1 13 12
+3 1 14 13
+3 1 15 14
+3 1 16 15
+3 1 17 16
+1 1 9 0
+2 1 8 17
+
+vertices
+18
+2
+0 0
+1 0
+2 0
+3 0
+4 0
+5 0
+6 0
+7 0
+8 0
+0 1
+1 1
+2 1
+3 1
+4 1
+5 1
+6 1
+7 1
+8 1

test/test_workflow.cpp

@@ -11,6 +11,7 @@ using namespace mfem;
 
 static const double threshold = 1.0e-14;
 static const double stokes_threshold = 1.0e-12;
+static const double linelast_threshold = 1.0e-13;
 
 /**
  * Simple smoke test to make sure Google Test is properly linked
@@ -553,6 +554,68 @@ TEST(SteadyNS_Workflow, ComponentSeparateVariable)
    return;
 }
 
+TEST(LinElast_Workflow, MFEMIndividualTest)
+{
+   config = InputParser("inputs/linelast.base.yml");
+
+   config.dict_["model_reduction"]["rom_handler_type"] = "mfem";
+   config.dict_["sample_generation"]["parameters"][0]["sample_size"] = 4;
+   config.dict_["basis"]["number_of_basis"] = 2; 
+   for (int k = 0; k < 2; k++)
+      config.dict_["basis"]["tags"][k]["name"] = "dom" + std::to_string(k);
+
+   config.dict_["main"]["mode"] = "sample_generation";
+   GenerateSamples(MPI_COMM_WORLD);
+
+   config.dict_["main"]["mode"] = "train_rom";
+   TrainROM(MPI_COMM_WORLD);
+   config.dict_["main"]["mode"] = "build_rom";
+   BuildROM(MPI_COMM_WORLD);
+   config.dict_["main"]["mode"] = "single_run";
+   double error = SingleRun(MPI_COMM_WORLD);
+
+   // This reproductive case must have a very small error at the level of finite-precision.
+   printf("Error: %.15E\n", error);
+   EXPECT_TRUE(error < linelast_threshold);
+
+   return;
+}
+
+TEST(LinElast_Workflow, MFEMUniversalTest)
+{
+   config = InputParser("inputs/linelast.base.yml");
+
+   config.dict_["model_reduction"]["rom_handler_type"] = "mfem";
+
+   config.dict_["single_run"]["linelast_disp"]["rdisp_f"] = 0.9;
+   config.dict_["sample_generation"]["parameters"][0]["sample_size"] = 2;
+   config.dict_["model_reduction"]["subdomain_training"] = "universal";
+   config.dict_["basis"]["number_of_basis"] = 4; 
+
+   // Test save/loadSolution as well.
+   config.dict_["save_solution"]["enabled"] = true;
+   config.dict_["model_reduction"]["compare_solution"]["load_solution"] = true;
+   config.dict_["model_reduction"]["compare_solution"]["fom_solution_file"] = "./sample0_solution.h5";
+
+   config.dict_["main"]["mode"] = "sample_generation";
+   GenerateSamples(MPI_COMM_WORLD); 
+
+   config.dict_["main"]["mode"] = "train_rom";
+   TrainROM(MPI_COMM_WORLD);
+   config.dict_["main"]["mode"] = "build_rom";
+   BuildROM(MPI_COMM_WORLD);
+
+   config.dict_["save_solution"]["enabled"] = false;
+   config.dict_["main"]["mode"] = "single_run";
+   double error = SingleRun(MPI_COMM_WORLD);
+
+   // This reproductive case must have a very small error at the level of finite-precision.
+   printf("Error: %.15E\n", error);
+   EXPECT_TRUE(error < linelast_threshold); 
+
+   return;
+}
+
 int main(int argc, char* argv[])
 {
    ::testing::InitGoogleTest(&argc, argv);