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/test_workflow.cpp

@@ -553,6 +553,73 @@ 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_["model_reduction"]["visualization"]["enabled"] = true;
+   config.dict_["model_reduction"]["visualization"]["prefix"] = "basis_paraview";
+   config.dict_["sample_generation"]["parameters"][0]["sample_size"] = 2;
+   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 < threshold);
+
+   return;
+}
+
+TEST(LinElast_Workflow, MFEMUniversalTest)
+{
+   config = InputParser("inputs/linelast.base.yml");
+   config.dict_["visualization"]["enabled"] = true;
+
+   config.dict_["model_reduction"]["rom_handler_type"] = "mfem";
+   config.dict_["model_reduction"]["visualization"]["enabled"] = true;
+   config.dict_["model_reduction"]["visualization"]["prefix"] = "basis_paraview";
+
+   config.dict_["single_run"]["linelast_disp"]["rdisp_f"] = 1.0;
+   config.dict_["sample_generation"]["parameters"][0]["sample_size"] = 1;
+   config.dict_["model_reduction"]["subdomain_training"] = "universal";
+   config.dict_["basis"]["number_of_basis"] = 2; 
+
+   // 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 < threshold); 
+
+   return;
+}
+
 int main(int argc, char* argv[])
 {
    ::testing::InitGoogleTest(&argc, argv);