Provide mechanism to obtain scaling factor used in constant flow rate. (

#189)

examples/turbPipePeriodic/turbPipe.udf

@@ -54,6 +54,21 @@ void UDF_Setup(nrs_t *nrs)
 
 void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
 {
+  if(platform->options.compareArgs("CONSTANT FLOW RATE", "TRUE")){
+    mesh_t * mesh = nrs->meshV;
+    occa::memory o_Uz = nrs->o_U + 2 * nrs->fieldOffset * sizeof(dfloat);
+    const dfloat ubar = platform->linAlg->innerProd(mesh->Nlocal, o_Uz, mesh->o_LMM, platform->comm.mpiComm)/mesh->volume;
+
+    dfloat expectedUbar;
+    platform->options.getArgs("FLOW RATE", expectedUbar);
+
+    if(platform->comm.mpiRank == 0){
+      printf("constantFlowScale = %g, flowRateError = %g\n",
+        ConstantFlowRate::scaleFactor(),
+        ubar - expectedUbar);
+    }
+    fflush(stdout);
+  }
   if (nrs->lastStep)
     if (ciMode) ciTestErrors(nrs, time, tstep);
 }

src/lns/constantFlowRate.cpp

@@ -5,6 +5,7 @@
 #include <limits>
 
 namespace {
+static dfloat constantFlowScale = 0.0;
 inline dfloat distance(
     dfloat x1, dfloat x2, dfloat y1, dfloat y2, dfloat z1, dfloat z2) {
   const dfloat dist_x = x1 - x2;
@@ -330,21 +331,25 @@ bool apply(nrs_t *nrs, int tstep, dfloat time) {
 
   const dfloat deltaFlowRate = volumetricFlowRate - currentFlowRate;
 
-  const dfloat scale = deltaFlowRate / baseFlowRate;
+  constantFlowScale = deltaFlowRate / baseFlowRate;
 
   // vx += scale * vxc
   platform->linAlg->axpbyMany(mesh->Nlocal,
       nrs->NVfields,
       nrs->fieldOffset,
-      scale,
+      constantFlowScale,
       nrs->o_Uc,
       1.0,
       nrs->o_U);
-  platform->linAlg->axpby(mesh->Nlocal, scale, nrs->o_Pc, 1.0, nrs->o_P);
+  platform->linAlg->axpby(mesh->Nlocal, constantFlowScale, nrs->o_Pc, 1.0, nrs->o_P);
 
   return recomputeBaseFlowRate;
 }
 
+dfloat scaleFactor(){
+  return constantFlowScale;
+}
+
 void compute(nrs_t *nrs, double lengthScale, dfloat time) {
 
   constexpr int ndim = 3;

src/lns/constantFlowRate.hpp

@@ -7,6 +7,7 @@ namespace ConstantFlowRate{
 bool apply(nrs_t *nrs, int tstep, dfloat time);
 void compute(nrs_t *nrs, dfloat lengthScale, dfloat time);
 bool checkIfRecompute(nrs_t* nrs, int tstep);
+dfloat scaleFactor();
 }
 
 #endif

src/udf/udf.hpp

@@ -6,6 +6,7 @@
 #include "nrs.hpp"
 #include "nekInterfaceAdapter.hpp"
 #include "parReader.hpp"
+#include "constantFlowRate.hpp"
 
 extern "C" {
 void UDF_Setup0(MPI_Comm comm, setupAide &options);