From 03194998151c56a1c4ce70e5d5b5e440dd6330f9 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Wed, 3 Nov 2021 10:59:32 +0100
Subject: [PATCH] added: save anasols to HDF5

---
 src/SIM/SIMbase.C           | 62 +++++++++++++++++++++++++++++++++++++
 src/SIM/SIMbase.h           |  3 ++
 src/Utility/ExprFunctions.C |  2 +-
 src/Utility/ExprFunctions.h |  9 ++++++
 src/Utility/HDF5Writer.C    | 41 ++++++++++++++++++++++++
 src/Utility/HDF5Writer.h    |  6 ++++
 6 files changed, 122 insertions(+), 1 deletion(-)

diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index 9460747df..cdef7176e 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -24,6 +24,7 @@
 #include "AlgEqSystem.h"
 #include "LinSolParams.h"
 #include "EigSolver.h"
+#include "ExprFunctions.h"
 #include "GlbNorm.h"
 #include "GlbL2projector.h"
 #include "ElmNorm.h"
@@ -2400,3 +2401,64 @@ void SIMbase::dumpEqSys ()
       utl::zero_print_tol = old_tol;
     }
 }
+
+
+std::map<std::string, std::string> SIMbase::getAnaSolDefinitions () const
+{
+  std::map<std::string, std::string> result;
+
+  if (mySol) {
+    size_t nf = this->getNoFields(0);
+
+    if (nf == 1) // Scalar solution
+    {
+      const EvalFunction* pSol = dynamic_cast<const EvalFunction*>(mySol->getScalarSol());
+      if (pSol)
+        result.insert(std::make_pair(myProblem->getField1Name(0,"Exact"), pSol->asString()));
+
+      const VecFuncExpr* sSol = dynamic_cast<const VecFuncExpr*>(mySol->getScalarSecSol());
+      if (sSol)
+        for (size_t i = 1; i <= this->getNoSpaceDim(); ++i)
+          result.insert(std::make_pair(myProblem->getField2Name(i-1,"Exact"),sSol->asString(i)));
+    }
+    else
+    {
+      const VecFuncExpr* vSol = dynamic_cast<const VecFuncExpr*>(mySol->getVectorSol());
+      size_t idx = 0;
+      if (vSol) {
+        for (; idx < this->getNoSpaceDim(); ++idx)
+          result.insert(std::make_pair(myProblem->getField1Name(idx,"Exact"),vSol->asString(idx+1)));
+      }
+      size_t scalSol = 0;
+      while (mySol->getScalarSol(scalSol)) {
+        const EvalFunction* pSol = dynamic_cast<const EvalFunction*>(mySol->getScalarSol(scalSol));
+        if (pSol)
+          result.insert(std::make_pair(myProblem->getField1Name(idx++,"Exact"), pSol->asString()));
+        ++scalSol;
+      }
+
+      idx = 0;
+      const STensorFuncExpr* stSol = dynamic_cast<const STensorFuncExpr*>(mySol->getStressSol());
+      if (stSol)
+        for (size_t i = 1; i <= stSol->dim(); ++i, ++idx)
+          result.insert(std::make_pair(myProblem->getField2Name(idx,"Exact"),stSol->asString(i)));
+
+      const TensorFuncExpr* tSol = dynamic_cast<const TensorFuncExpr*>(mySol->getVectorSecSol());
+      if (tSol)
+        for (size_t i = 1; i <= tSol->dim(); ++i, ++idx)
+          result.insert(std::make_pair(myProblem->getField2Name(idx,"Exact"),tSol->asString(i)));
+
+      scalSol = 0;
+      while (mySol->getScalarSecSol(scalSol)) {
+        const VecFuncExpr* sSol = dynamic_cast<const VecFuncExpr*>(mySol->getScalarSecSol(scalSol));
+        if (sSol) {
+          for (size_t i = 1; i <= this->getNoSpaceDim(); ++i, ++idx)
+            result.insert(std::make_pair(myProblem->getField2Name(idx,"Exact"),sSol->asString(i)));
+        }
+        ++scalSol;
+      }
+    }
+  }
+
+  return result;
+}
diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h
index 5737865c6..d492acf50 100644
--- a/src/SIM/SIMbase.h
+++ b/src/SIM/SIMbase.h
@@ -540,6 +540,9 @@ class SIMbase : public SIMadmin, public SIMdependency
   //! \brief Returns whether a dual solution is available or not.
   virtual bool haveDualSol() const { return dualField ? true : false; }
 
+  //! \brief Returns analytic solutions as their string definitions if possible.
+  std::map<std::string, std::string> getAnaSolDefinitions() const;
+
   //! \brief Returns a pointer to a norm integrand object for this simulator.
   //! \note The object is allocated dynamically and has therefore to be
   //! manually deleted before the variable receiving the pointer value goes
diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C
index 02d1e79e6..af3fc9dc2 100644
--- a/src/Utility/ExprFunctions.C
+++ b/src/Utility/ExprFunctions.C
@@ -181,7 +181,7 @@ Real EvalFunc::deriv (Real x) const
 }
 
 
-EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
+EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}, definition(function)
 {
   try {
 #ifdef USE_OPENMP
diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h
index 5757e6804..c24bf79e6 100644
--- a/src/Utility/ExprFunctions.h
+++ b/src/Utility/ExprFunctions.h
@@ -123,6 +123,9 @@ class EvalFunction : public RealFunc
   //! \brief Set an additional parameter in the function.
   void setParam(const std::string& name, double value);
 
+  //! \brief Returns function string.
+  const std::string& asString() const { return definition; }
+
 protected:
   //! \brief Non-implemented copy constructor to disallow copying.
   EvalFunction(const EvalFunction&) = delete;
@@ -133,6 +136,8 @@ class EvalFunction : public RealFunc
 
   //! \brief Cleans up the allocated data.
   void cleanup();
+
+  std::string definition; //!< Function definition, used for serialization purposes.
 };
 
 
@@ -201,6 +206,10 @@ class EvalMultiFunction : public ParentFunc, public EvalFunctions
       func->setParam(name, value);
   }
 
+  //! \brief Returns function string for component.
+  //! \param idx 1-based index for component to return.
+  const std::string& asString(size_t idx) const { return p[idx-1]->asString(); }
+
 protected:
   //! \brief Sets the number of spatial dimensions (default implementation).
   void setNoDims() { ParentFunc::ncmp = nsd = p.size(); }
diff --git a/src/Utility/HDF5Writer.C b/src/Utility/HDF5Writer.C
index 7b3f65ac1..afa4208a4 100644
--- a/src/Utility/HDF5Writer.C
+++ b/src/Utility/HDF5Writer.C
@@ -511,6 +511,9 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
     }
   }
 
+  if (sim->haveAnaSol())
+    this->writeAnaSol(sim->getAnaSolDefinitions(), sim->getName());
+
   delete norm;
   for (hid_t g : group)
     if (g != -1)
@@ -783,3 +786,41 @@ bool HDF5Writer::writeLog (const std::string& data, const std::string& name)
   return false;
 #endif
 }
+
+
+bool HDF5Writer::writeAnaSol (const std::map<std::string,std::string>& aSols,
+                              const std::string& name)
+{
+  if (aSols.empty())
+    return true;
+
+#ifdef HAS_HDF5
+  hid_t group;
+  if (!checkGroupExistence(m_file,"/anasol")) {
+    group = H5Gcreate2(m_file,"/anasol",0,H5P_DEFAULT,H5P_DEFAULT);
+    H5Gclose(group);
+  }
+
+  std::stringstream str;
+  str << "/anasol/" << name;
+
+  if (checkGroupExistence(m_file,str.str().c_str()))
+    group = H5Gopen2(m_file,str.str().c_str(),H5P_DEFAULT);
+  else
+    group = H5Gcreate2(m_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT);
+
+  int rank = 0;
+#ifdef HAVE_MPI
+  MPI_Comm_rank(*m_adm.getCommunicator(), &rank);
+#endif
+
+  for (const auto& it : aSols)
+    this->writeArray(group,it. first, -1,
+                     rank == 0 ? it.second.size() : 0, it.second.data(), H5T_NATIVE_CHAR);
+
+  H5Gclose(group);
+  return true;
+#else
+  return false;
+#endif
+}
diff --git a/src/Utility/HDF5Writer.h b/src/Utility/HDF5Writer.h
index b65a69ad7..3e7f02807 100644
--- a/src/Utility/HDF5Writer.h
+++ b/src/Utility/HDF5Writer.h
@@ -97,6 +97,12 @@ class HDF5Writer : public DataWriter, public HDF5Base
   //! \param name Name of log
   virtual bool writeLog(const std::string& data, const std::string& name);
 
+  //! \brief Write analytical solutions to file.
+  //! \param aSol Map of functions
+  //! \param name Name of simulator
+  bool writeAnaSol(const std::map<std::string, std::string>& aSols,
+                   const std::string& name);
+
 #ifdef HAS_HDF5
 protected:
   //! \brief Internal helper function writing a data array to file.