Reset and precomputation

LLNL · Jul 20, 2023 · 4a77ef0 · 4a77ef0
1 parent aebb5c7
commit 4a77ef0
Showing 1 changed file with 18 additions and 1 deletion.
diff --git a/examples/prom/nonlinear_elasticity_global_rom.cpp b/examples/prom/nonlinear_elasticity_global_rom.cpp
@@ -153,6 +153,8 @@ class RomOperator : public TimeDependentOperator
     mutable Vector *zfom_v;
     CAROM::Vector *zfom_x_librom;
 
+    CAROM::Vector *v0_fom_librom;
+
     CAROM::SampleMeshManager *smm;
 
     CAROM::Vector *z_v_librom;
@@ -203,7 +205,8 @@ class RomOperator : public TimeDependentOperator
     void Mult_FullOrder(const Vector &y, Vector &dy_dt) const;
     void SetEQP(CAROM::Vector *eqpSol);
 
-    CAROM::Matrix V_v, V_x, V_vTU_H;
+    CAROM::Matrix V_v, V_x, V_vTU_H, V_xTV_v;
+    CAROM::Vector V_xTV_vv0;
     const Vector *x0;
     const Vector *v0;
     Vector v0_fom;
@@ -1655,6 +1658,13 @@ RomOperator::RomOperator(HyperelasticOperator *fom_,
     // Invert M_hat and store
     M_hat->inverse(*M_hat_inv);
 
+    // Basis projection
+    V_x.transposeMult(V_v, V_xTV_v);
+    // Initial velocity
+    v0_fom_librom = new CAROM::Vector(v0_fom.GetData(),
+                                             v0_fom.Size(), true, false);
+    V_xTV_v.mult(*v0_fom_librom, V_xTV_vv0);
+
     if (myid == 0 && hyperreduce)
     {
         if (!eqp)
@@ -1768,6 +1778,11 @@ void RomOperator::Mult_Hyperreduced(const Vector &vx, Vector &dvx_dt) const
         V_v.mult(v_librom, *z_v_librom);
         add(z_v, *v0, *pfom_v);
         V_x.transposeMult(*pfom_v_librom, dx_dt_librom);
+
+        // TODO: test this when build succeeds
+        //V_xTV_v.mult(v_librom, *z_v_librom);
+        //V_xTV_vv0.add(*z_v_librom, dx_dt_librom);
+
         Vector resEQP;
         if (fastIntegration)
         {
@@ -1891,6 +1906,8 @@ void RomOperator::Mult(const Vector &vx, Vector &dvx_dt) const
         Mult_FullOrder(vx, dvx_dt);
 }
 
+// Set EQP coefficients of the ROM operator. 
+// This must be done at every time window change.
 void RomOperator::SetEQP(CAROM::Vector *eqpSol)
 {
     std::set<int> elements;