Poisson rom (#257)

* rom workflow for buildROMPoissonOperator

* change from total to hartree potential

* test routine for Poisson ROM operator.

* POD reconstruction

* test Poisson ROM operator verified.

* buildROMPoissonOperator updated.

* ci test for Poisson ROM operator

.github/workflows/ci.yml

@@ -35,15 +35,21 @@ jobs:
             submodules: 'true'
       - name: cmake
         run: |
-            mkdir build
-            cd build
+            mkdir ${GITHUB_WORKSPACE}/build
+            cd ${GITHUB_WORKSPACE}/build
             cmake .. -DCMAKE_CXX_COMPILER=mpic++ -DCMAKE_Fortran_COMPILER=mpif90 -DMPIEXEC_PREFLAGS="--oversubscribe" -DUSE_LIBROM=On -DLIBROM_PATH=/env/dependencies/libROM
       - name: make
         run: |
-            cd build && make -j 4
+            cd ${GITHUB_WORKSPACE}/build && make -j 4
       - name: test
         run: |
-            cd build && ctest --no-compress-output -V -T Test -I 1,20,1
+            cd ${GITHUB_WORKSPACE}/build && ctest --no-compress-output -V -T Test -I 1,20,1
+      - name: test ROM Poisson operator
+        run: |
+            cd ${GITHUB_WORKSPACE}/tests/ROM/test_rom_poisson
+            ln -s ${GITHUB_WORKSPACE}/build/src/mgmol-rom .
+            ln -s ${GITHUB_WORKSPACE}/potentials/* .
+            mpirun -n 3 --oversubscribe ./mgmol-rom -c carbyne.cfg -i carbyne.in
   # code-style:
   #   runs-on: ubuntu-latest
   #   needs: [docker-image]

src/Control.cc

@@ -2069,6 +2069,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm)
             rom_pri_option.rom_stage = ROMStage::ONLINE;
         else if (str.compare("build") == 0)
             rom_pri_option.rom_stage = ROMStage::BUILD;
+        else if (str.compare("test_poisson") == 0)
+            rom_pri_option.rom_stage = ROMStage::TEST_POISSON;
         else if (str.compare("none") == 0)
             rom_pri_option.rom_stage = ROMStage::UNSUPPORTED;
 
@@ -2086,6 +2088,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm)
             rom_pri_option.variable = ROMVariable::NONE;
 
         rom_pri_option.save_librom_snapshot = vm["ROM.offline.save_librom_snapshot"].as<bool>();
+
+        rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as<int>();
     }  // onpe0
 
     // synchronize all processors
@@ -2130,4 +2134,7 @@ void Control::syncROMOptions()
 
     rom_pri_option.rom_stage = static_cast<ROMStage>(rom_stage);
     rom_pri_option.variable = static_cast<ROMVariable>(rom_var);
+
+    mpirc = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_);
+    bcast_check(mpirc);
 }

src/Electrostatic.h

@@ -46,6 +46,8 @@ class Electrostatic
     ~Electrostatic();
     static Timer solve_tm() { return solve_tm_; }
 
+    Poisson* getPoissonSolver() { return poisson_solver_; }
+
     void setup(const short max_sweeps);
     void setupPB(const double e0, const double rho0, const double drho0,
         Potentials& pot);

src/Hartree_CG.h

@@ -40,6 +40,13 @@ class Hartree_CG : public Poisson
 
     void solve(const pb::GridFunc<RHODTYPE>& rho,
         const pb::GridFunc<RHODTYPE>& rhoc) override;
+
+    void applyOperator(pb::GridFunc<POTDTYPE> &vh,
+        pb::GridFunc<POTDTYPE> &lhs) override
+    {
+        T *oper = poisson_solver_->getOperator();
+        oper->apply(vh, lhs);
+    }
 };
 
 #endif

src/PCGSolver.h

@@ -97,6 +97,8 @@ class PCGSolver
     double getFinalResidual() const { return final_residual_; }
     double getResidualReduction() const { return residual_reduction_; }
 
+    T* getOperator() { return &oper_; }
+
     // Destructor
     ~PCGSolver() { clear(); }
 };

src/Poisson.h

@@ -88,6 +88,12 @@ class Poisson : public PoissonInterface
         Int_vhrho_       = vel * vh_->gdot(rho);
         Int_vhrhoc_      = vel * vh_->gdot(rhoc);
     }
+
+    virtual void applyOperator(pb::GridFunc<POTDTYPE> &vh, pb::GridFunc<POTDTYPE> &lhs)
+    {
+        std::cerr << "ERROR: Abstract method Poisson::applyOperator()" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 0);
+    }
 };
 
 #endif

src/read_config.cc

@@ -427,6 +427,8 @@ void setupROMConfigOption(po::options_description &rom_cfg)
         ("ROM.offline.save_librom_snapshot", po::value<bool>()->default_value(false),
             "Save libROM snapshot file at FOM simulation.")
         ("ROM.offline.variable", po::value<std::string>()->default_value(""),
-            "FOM variable to perform POD: either orbitals or potential.");
+            "FOM variable to perform POD: either orbitals or potential.")
+        ("ROM.basis.number_of_potential_basis", po::value<int>()->default_value(-1),
+            "Number of potential POD basis to build Hartree potential ROM operator.");
 }
 #endif

src/rom_Control.h

@@ -22,6 +22,7 @@ enum class ROMStage
     ONLINE,
     RESTORE,    // TODO(kevin): what stage is this?
     BUILD,
+    TEST_POISSON,
     UNSUPPORTED
 };
 
@@ -45,6 +46,9 @@ struct ROMPrivateOptions
 
     /* save librom snapshot matrix at FOM simulation. */
     bool save_librom_snapshot = false;
+
+    /* options for ROM building */
+    int num_potbasis = -1;
 };
 
 #endif  // ROM_CONTROL_H

src/rom_main.cc

@@ -86,6 +86,10 @@ int main(int argc, char** argv)
     mmpi.bcastGlobal(input_filename);
     mmpi.bcastGlobal(lrs_filename);
 
+    /* Get ROM driver mode */
+    ROMPrivateOptions rom_options = ct.getROMOptions();
+    ROMStage rom_stage = rom_options.rom_stage;
+
     // Enter main scope
     {
         MGmolInterface* mgmol;
@@ -98,10 +102,36 @@ int main(int argc, char** argv)
 
         mgmol->setup();
 
-        if (ct.isLocMode())
-            readRestartFiles<LocGridOrbitals>(mgmol);
-        else
-            readRestartFiles<ExtendedGridOrbitals>(mgmol);
+        switch (rom_stage)
+        {
+        case (ROMStage::OFFLINE):
+            if (ct.isLocMode())
+                readRestartFiles<LocGridOrbitals>(mgmol);
+            else
+                readRestartFiles<ExtendedGridOrbitals>(mgmol);
+        break;
+
+        case (ROMStage::BUILD):
+            if (ct.isLocMode())
+                buildROMPoissonOperator<LocGridOrbitals>(mgmol);
+            else
+                buildROMPoissonOperator<ExtendedGridOrbitals>(mgmol);
+        break;
+
+        case (ROMStage::TEST_POISSON):
+            if (ct.isLocMode())
+                testROMPoissonOperator<LocGridOrbitals>(mgmol);
+            else
+                testROMPoissonOperator<ExtendedGridOrbitals>(mgmol);
+        break;
+
+        default:
+            std::cerr << "rom_main error: Unknown ROM stage" << std::endl;
+            MPI_Abort(MPI_COMM_WORLD, 0);
+        break;
+        }
+
+
 
         delete mgmol;