Merge pull request #87 from Pressio/workflow_code_85

Workflow code

tpls/pressio/commit.txt

@@ -0,0 +1,8 @@
+commit 020e9e9454b3b772f8630221e17223cc5e270c29 (HEAD -> develop, origin/develop)
+Merge: 54fd7884 9e13fe03
+Author: Francesco Rizzi <[email protected]>
+Date:   Wed Feb 2 19:28:05 2022 +0100
+
+    Merge pull request #335 from Pressio/observe_rhs_for_explicit_steppers_333
+
+    ode: support RHS observer for explicit steppers

tpls/pressio/include/pressio/ode_steppers/impl/ode_explicit_stepper.hpp

@@ -128,48 +128,73 @@ class ExplicitStepper
 
   template<class StepCountType>
   void operator()(StateType & odeState,
-		  const ScalarType & time,
+		  const ScalarType & currentTime,
 		  const ScalarType & dt,
-		  const StepCountType & step)
+		  const StepCountType & stepNumber)
+  {
+    auto dummyRhsObserver = [](const StepCountType &,
+			       const ScalarType &,
+			       const VelocityType &) { /*no op*/ };
+    (*this)(odeState, currentTime, dt, stepNumber, dummyRhsObserver);
+  }
+
+  template<class StepCountType, class RhsObserverType>
+  void operator()(StateType & odeState,
+		  const ScalarType & currentTime,
+		  const ScalarType & dt,
+		  const StepCountType & step,
+		  RhsObserverType & rhsObs)
   {
     if (name_ == ode::StepScheme::ForwardEuler){
-      doStepImpl(ode::ForwardEuler(), odeState, time, dt, step);
+      doStepImpl(ode::ForwardEuler(), odeState,
+		 currentTime, dt, step, rhsObs);
     }
+
     else if (name_ == ode::StepScheme::RungeKutta4){
-      doStepImpl(ode::RungeKutta4(), odeState, time, dt, step);
+      doStepImpl(ode::RungeKutta4(), odeState,
+		 currentTime, dt, step, rhsObs);
     }
+
     else if (name_ == ode::StepScheme::AdamsBashforth2){
-      doStepImpl(ode::AdamsBashforth2(), odeState, time, dt, step);
+      doStepImpl(ode::AdamsBashforth2(), odeState,
+		 currentTime, dt, step, rhsObs);
     }
+
     else if (name_ == ode::StepScheme::SSPRungeKutta3){
-      doStepImpl(ode::SSPRungeKutta3(), odeState, time, dt, step);
+      doStepImpl(ode::SSPRungeKutta3(), odeState,
+		 currentTime, dt, step, rhsObs);
     }
   }
 
 private:
-  template<class StepCountType>
+  template<class StepCountType, class RhsObserverType>
   void doStepImpl(ode::ForwardEuler,
 		  StateType & odeState,
-		  const ScalarType & time,
+		  const ScalarType & currentTime,
 		  const ScalarType & dt,
-		  const StepCountType & step)
+		  const StepCountType & stepNumber,
+		  RhsObserverType & rhsObs)
   {
     PRESSIOLOG_DEBUG("euler forward stepper: do step");
 
     auto & rhs = velocities_[0];
-    //eval RHS
-    systemObj_.get().velocity(odeState, time, rhs);
+
+    //eval and observe RHS
+    systemObj_.get().velocity(odeState, currentTime, rhs);
+    rhsObs(stepNumber, currentTime, rhs);
+
     // y = y + dt * rhs
     constexpr auto one  = ::pressio::utils::Constants<ScalarType>::one();
     ::pressio::ops::update(odeState, one, rhs, dt);
   }
 
-  template<class StepCountType>
+  template<class StepCountType, class RhsObserverType>
   void doStepImpl(ode::AdamsBashforth2,
 		  StateType & odeState,
 		  const ScalarType & currentTime,
 		  const ScalarType & dt,
-		  const StepCountType & stepNumber)
+		  const StepCountType & stepNumber,
+		  RhsObserverType & rhsObs)
   {
     PRESSIOLOG_DEBUG("adams-bashforth2 stepper: do step");
 
@@ -182,6 +207,7 @@ class ExplicitStepper
       // use Euler forward or we could use something else here maybe RK4
       auto & rhs = velocities_[0];
       systemObj_.get().velocity(odeState, currentTime, rhs);
+      rhsObs(stepNumber, currentTime, rhs);
 
       constexpr auto one   = ::pressio::utils::Constants<ScalarType>::one();
       ::pressio::ops::update(odeState, one, rhs, dt);
@@ -193,17 +219,19 @@ class ExplicitStepper
       ::pressio::ops::deep_copy(fnm1, fn);
 
       systemObj_.get().velocity(odeState, currentTime, fn);
+      rhsObs(stepNumber, currentTime, fn);
       constexpr auto one   = ::pressio::utils::Constants<ScalarType>::one();
       ::pressio::ops::update(odeState, one, fn, cfn, fnm1, cfnm1);
     }
   }
 
-  template<class StepCountType>
+  template<class StepCountType, class RhsObserverType>
   void doStepImpl(ode::SSPRungeKutta3,
 		  StateType & odeState,
-		  const ScalarType & time,
+		  const ScalarType & currentTime,
 		  const ScalarType & dt,
-		  const StepCountType & step)
+		  const StepCountType & stepNumber,
+		  RhsObserverType & rhsObs)
   {
     PRESSIOLOG_DEBUG("ssprk3 stepper: do step");
 
@@ -223,33 +251,35 @@ class ExplicitStepper
     // see e.g. https://gkeyll.readthedocs.io/en/latest/dev/ssp-rk.html
 
     // rhs(u_n, t_n)
-    systemObj_.get().velocity(odeState, time, rhs0);
+    systemObj_.get().velocity(odeState, currentTime, rhs0);
     // u_1 = u_n + dt * rhs(u_n, t_n)
     ::pressio::ops::update(auxiliaryState_, zero,
-                           odeState,     one,
+                           odeState,        one,
                            rhs0,            dt);
+    rhsObs(stepNumber, currentTime, rhs0);
 
     // rhs(u_1, t_n+dt)
-    systemObj_.get().velocity(auxiliaryState_, time+dt, rhs0);
+    systemObj_.get().velocity(auxiliaryState_, currentTime+dt, rhs0);
     // u_2 = 3/4*u_n + 1/4*u_1 + 1/4*dt*rhs(u_1, t_n+dt)
     ::pressio::ops::update(auxiliaryState_, fourInv,
-			   odeState,     threeOvFour,
+			   odeState,        threeOvFour,
                            rhs0,            fourInv*dt);
 
     // rhs(u_2, t_n + 0.5*dt)
-    systemObj_.get().velocity(auxiliaryState_, time + oneOvTwo*dt, rhs0);
+    systemObj_.get().velocity(auxiliaryState_, currentTime + oneOvTwo*dt, rhs0);
     // u_n+1 = 1/3*u_n + 2/3*u_2 + 2/3*dt*rhs(u_2, t_n+0.5*dt)
-    ::pressio::ops::update(odeState,     oneOvThree,
+    ::pressio::ops::update(odeState,        oneOvThree,
 			   auxiliaryState_, twoOvThree,
                            rhs0,            twoOvThree*dt);
   }
 
-  template<class StepCountType>
+  template<class StepCountType, class RhsObserverType>
   void doStepImpl(ode::RungeKutta4,
 		  StateType & odeState,
-		  const ScalarType & t,
+		  const ScalarType & currentTime,
 		  const ScalarType & dt,
-		  const StepCountType & step)
+		  const StepCountType & stepNumber,
+		  RhsObserverType & rhsObs)
   {
     PRESSIOLOG_DEBUG("rk4 stepper: do step");
 
@@ -263,12 +293,13 @@ class ExplicitStepper
     constexpr auto six  = two * three;
 
     const ScalarType dt_half = dt / two;
-    const ScalarType t_phalf = t + dt_half;
+    const ScalarType t_phalf = currentTime + dt_half;
     const ScalarType dt6 = dt / six;
     const ScalarType dt3 = dt / three;
 
     // stage 1: ytmp = y + rhs0*dt_half;
-    systemObj_.get().velocity(odeState, t, rhs0);
+    systemObj_.get().velocity(odeState, currentTime, rhs0);
+    rhsObs(stepNumber, currentTime, rhs0);
     this->rk4_stage_update_impl(auxiliaryState_, odeState, rhs0, dt_half);
 
     // stage 2: ytmp = y + rhs1*dt_half;
@@ -280,7 +311,7 @@ class ExplicitStepper
     this->rk4_stage_update_impl(auxiliaryState_, odeState, rhs2, dt);
 
     // stage 4: y_n += dt/6 * ( k1 + 2 * k2 + 2 * k3 + k4 )
-    systemObj_.get().velocity(auxiliaryState_,  t + dt, rhs3);
+    systemObj_.get().velocity(auxiliaryState_, currentTime + dt, rhs3);
     this->rk4_stage_update_impl(odeState, rhs0, rhs1, rhs2, rhs3, dt6, dt3);
   }
 

workflow_codes/2d_gs/__init__.py

@@ -0,0 +1,10 @@
+
+from .constants import *
+from .scenarios import *
+
+from .specialize_execute_foms import *
+from .specialize_extract_fom_states import *
+from .specialize_compute_pod import *
+from .specialize_compute_sample_meshes import *
+from .specialize_execute_galerkin_roms import *
+from .specialize_rom_reconstruction_and_error import *

workflow_codes/2d_gs/constants.py

@@ -0,0 +1,4 @@
+
+problemName    = "2d_gs"
+dimensionality = 2
+numDofsPerCell = 2

workflow_codes/2d_gs/impl.py

@@ -0,0 +1,111 @@
+#!/usr/bin/env python
+
+#==============================================================
+# imports
+#==============================================================
+
+import copy
+
+# LOCAL imports
+from .scenarios import \
+  base_dic, test_points, train_points
+from .constants import \
+  problemName, numDofsPerCell
+
+
+#==============================================================
+# functions
+#==============================================================
+
+def _replace_dic_field(scenario, dic, val):
+  if scenario == 1:
+    dic['physicalCoefficients']['feedRate'] = float(val)
+
+  # !!!!!
+  # if you add another scenario, add specifics here
+  # !!!!!
+
+  else:
+    sys.exit("2d_gs: scenario = {} not supported yet".format(scenario))
+
+
+def _create_fom_input_dics_impl(scenario, meshPath, kind):
+  assert(meshPath != None)
+
+  values = train_points[scenario] if kind=="train" else test_points[scenario]
+
+  result = []
+  for runNumber, v in values.items():
+    # make a copy of the base dic
+    thisDic = copy.deepcopy(base_dic[scenario])
+
+    thisDic['runId'] = runNumber
+
+    # we are making dics for the FOM, so remove the rom key
+    if 'rom' in thisDic: del thisDic['rom']
+
+    # always need to set mesh and problem
+    thisDic['general']['meshDir'] = meshPath
+    thisDic['general']['problem'] = problemName
+    thisDic['general']['numDofsPerCell'] = numDofsPerCell
+
+    _replace_dic_field(scenario, thisDic, v)
+    result.append(thisDic)
+
+  return result
+
+#============================================================
+def _create_rom_input_dics_impl(scenario, meshPath, algoName, \
+                                numModes, lsvFile, refStateFile, \
+                                projectorFile = None,\
+                                # we need an optional sampleMeshPath because
+                                # for masked galerkin the mesh to use is the FULL mesh
+                                # but we still need to use gids from the sample mesh
+                                sampleMeshPath = None):
+
+  # we only care about test points here
+  values = test_points[scenario]
+
+  result = []
+  for runNumber, v in values.items():
+    # make a copy of the base dic
+    thisDic = copy.deepcopy(base_dic[scenario])
+
+    # the run ID MUST be the runNumber here because it has to
+    # match the enumeration of the FOM
+    thisDic['runId'] = runNumber
+
+    # we are making dics for ROM, so remove the fom key
+    if 'fom' in thisDic: del thisDic['fom']
+
+    # always need to set mesh and problem
+    thisDic['general']['meshDir'] = meshPath
+    thisDic['general']['problem'] = problemName
+    thisDic['general']['numDofsPerCell'] = numDofsPerCell
+
+    thisDic['rom']['numModes'] = numModes
+    thisDic['rom']['podFile']  = lsvFile
+    thisDic['rom']['algoName'] = algoName
+    thisDic['rom']['referenceStateFile'] = refStateFile
+
+    if projectorFile != None:
+      # if the projector is nontrivial
+      thisDic['rom']['projectorFile'] = projectorFile
+
+      if "Masked" not in algoName:
+        # for real hyper-reduction, the mesh IS a sample mesh that contains gids
+        thisDic['rom']['sampleMeshGidsFile']  = meshPath+"/sample_mesh_gids.dat"
+        thisDic['rom']['stencilMeshGidsFile'] = meshPath+"/stencil_mesh_gids.dat"
+
+      elif "Masked" in algoName:
+        # if we are doing MASKED, then the mesh to use is a FULL mesh
+        # so the GIDs must be taken from the sampleMeshPath argument
+        assert( sampleMeshPath != None )
+        # gappy Galerkin needs to know about gids files
+        thisDic['rom']['sampleMeshGidsFile']  = sampleMeshPath+"/sample_mesh_gids.dat"
+        thisDic['rom']['stencilMeshGidsFile'] = sampleMeshPath+"/stencil_mesh_gids.dat"
+
+    _replace_dic_field(scenario, thisDic, v)
+    result.append(thisDic)
+
+  return result