(WIP) add neknek example with meshes coupling through solid

This is the minimium case to make nekrs support the one pebble case.
Need to extend nekrs to make this running
The cooresponding Nek5000 case runs and converge to the solution of the
single mesh

examples/chtNekNek/README.md

@@ -0,0 +1,65 @@
+## Conjugate Heat Transfer
+
+Run to steady state (measured by Dillin's limit.f)
+
+- Middle: y in [0,1]: 2D channel with vx=parabla. gap=1.0         
+  flow: v-O along x direction, double Walls for y        
+  heat: t-O       
+
+- Top: y in [1, 1.5]: solid      
+  heat: All I
+
+- Bottom: y in [-0.5, 0]: solid     
+  heat: All I
+
+
+### Cases
+
+- Baseline: re-learn neknek in Nek5000 2D
+
+  - run\_nek5k\_ref2d:
+    Single mesh case, cht
+
+  - run\_nek5k\_nn2d:    
+    neknek 2 meshes
+    - cht mesh: y in [-0.4, 1.4], E = 8 x (2+4+2)
+    - solid mesh: y in [-0.5, -0.25] U [1.25, 1.5],  E = 10x2 x 2 (top/bottom)
+
+
+  - run\_nek5k\_nn2d\_par:
+    neknek 2 meshes, use parfile
+    ```
+    ./neknekb cht_merge solid 1 1
+    ```
+
+- 3D Case: 
+
+  - run\_nek5k\_ref3d:
+
+  - run\_nek5k\_nn3d\_par:
+    ```
+    ./neknekb cht_mergeb solidb 1 1
+    ```
+
+  - run\_nekrs\_ref3d:
+
+  - run\_nekrs\_nn3d:
+    ```
+    ~/.local/nekrs_v23/bin/nrsbmpi conj_ht.sess 2
+    ```
+ 
+
+
+### Notes 
+- Nek5000
+  1. [Nek5000] neknek cannot have cht with ifflow=T in one session and ifflow=F in the other session
+  2. [Nek5000] Nek5000 cannot run ifflow=T + nelv=nelgv=0
+  3. [Nek5000] Because of the 1 and 2, we have to run ifflow=T with a fluid mesh that has zero solution (Walls)    
+     Note that, NekRS might have issues with residual norm = 0. maybe a sqrt(0)
+  4. [NekRS] neknek + CHT not supported !?
+     ```
+     Error in nrsSetup: Conjugate heat transfer + neknek not supported!
+      ```
+  5. [NekRS] idsess is not supported
+  6. [NekRS] cht number of boundary conditions different in t-meah and v-mesh causing error     
+     https://github.com/Nek5000/nekRS/issues/504

examples/chtNekNek/cht_merge/cht_merge.oudf

@@ -0,0 +1,40 @@
+// Boundary conditions
+void velocityDirichletConditions(bcData *bc)
+{
+  if (bc->id == 1) {
+    bc->u = 4.0 * bc->y * (1.0 - bc->y);
+    bc->v = 0.0;
+    bc->w = 0.0;
+  } else { 
+    bc->u = 0.0;
+    bc->v = 0.0;
+    bc->w = 0.0;
+  }
+}
+
+void scalarDirichletConditions(bcData *bc)
+{
+  if (bc->id == 1) {
+    bc->s = bc->sinterp;
+  } else {
+    bc->s = 0.0;
+  }
+}
+
+@kernel void cFill(const dlong Nelements,
+                   const dfloat CONST1,
+                   const dfloat CONST2,
+                   @ restrict const dlong *eInfo,
+                   @ restrict dfloat *QVOL)
+{
+  for (dlong e = 0; e < Nelements; ++e; @outer(0)) {
+    const dlong solid = eInfo[e];
+    for (int n = 0; n < p_Np; ++n; @inner(0)) {
+      const int id = e * p_Np + n;
+      QVOL[id] = CONST1;
+      if (solid) {
+        QVOL[id] = CONST2;
+      }
+    }
+  }
+}

examples/chtNekNek/cht_merge/cht_merge.par

@@ -0,0 +1,42 @@
+[GENERAL] 
+verbose = yes
+polynomialOrder = 7
+#startFrom = "restart.fld"
+numSteps = 10000
+dt = 5e-02
+timeStepper = tombo2
+
+writeControl = steps
+writeInterval = 1000
+
+[NEKNEK]
+boundaryEXTOrder = 1
+
+[PROBLEMTYPE]
+equation = navierStokes+variableViscosity
+
+[PRESSURE]
+residualTol = 1e-06
+
+[VELOCITY]
+boundaryTypeMap = v, O, W
+residualTol = 1e-08
+density = 1
+viscosity = 0.5
+initialGuess = projectionAconj+nVector=10
+
+[TEMPERATURE] # temperature with Hmholtz 
+boundaryTypeMap = int, t, O, I
+rhoCp = 1.0
+conductivity = 0.5
+residualTol = 1e-08
+initialGuess = projectionAconj+nVector=10
+
+#[SCALAR01] # temperature with CVODE
+#conjugateHeatTransfer = yes 
+#solver = cvode
+#absoluteTol = 1e-10
+#
+#[CVODE]
+#relativeTol = 1e-06
+#dtmax = 5e-02 

examples/chtNekNek/cht_merge/cht_merge.udf

@@ -0,0 +1,106 @@
+
+//
+// nekRS User Defined File
+//
+#include "udf.hpp"
+
+static int updateProperties = 1;
+static dfloat vtrans, vdiff, ttrans, tdiff;
+
+#ifdef __okl__
+
+#include "cht_merge.oudf"
+
+#endif
+
+void userq(nrs_t *nrs, dfloat time, occa::memory o_S, occa::memory o_FS)
+{
+  cds_t *cds = nrs->cds;
+  mesh_t *mesh = cds->mesh[0];
+  const dfloat qvolFluid = 0.0;
+  const dfloat qvolSolid = 1.0;
+  cFill(mesh->Nelements, qvolFluid, qvolSolid, mesh->o_elementInfo, o_FS);
+}
+
+void uservp(nrs_t *nrs,
+            dfloat time,
+            occa::memory o_U,
+            occa::memory o_S,
+            occa::memory o_UProp,
+            occa::memory o_SProp)
+{
+  cds_t *cds = nrs->cds;
+
+  if (updateProperties) {
+    const dfloat rho = vtrans;
+    const dfloat mue = vdiff;
+    const dfloat rhoCpFluid = ttrans;
+    const dfloat conFluid = tdiff;
+    const dfloat rhoCpSolid = 1.0;
+    const dfloat conSolid = 10.0 * tdiff;
+
+    if (platform->comm.mpiRank == 0) {
+      printf("Updating properties: v (%g,%g) tfluid (%g,%g) tsolid (%g,%g)\n",
+             rho,mue,rhoCpFluid,conFluid,rhoCpSolid,conSolid);
+    }
+
+    // velocity
+    const occa::memory o_mue = o_UProp.slice(0 * nrs->fieldOffset * sizeof(dfloat));
+    const occa::memory o_rho = o_UProp.slice(1 * nrs->fieldOffset * sizeof(dfloat));
+    cFill(nrs->meshV->Nelements, mue, 0, nrs->meshV->o_elementInfo, o_mue);
+    cFill(nrs->meshV->Nelements, rho, 0, nrs->meshV->o_elementInfo, o_rho);
+    // temperature
+    const occa::memory o_con = o_SProp.slice(0 * cds->fieldOffset[0] * sizeof(dfloat));
+    const occa::memory o_rhoCp = o_SProp.slice(1 * cds->fieldOffset[0] * sizeof(dfloat));
+    cFill(cds->mesh[0]->Nelements, conFluid, conSolid, cds->mesh[0]->o_elementInfo, o_con);
+    cFill(cds->mesh[0]->Nelements, rhoCpFluid, rhoCpSolid, cds->mesh[0]->o_elementInfo, o_rhoCp);
+    updateProperties = 0;
+  }
+}
+
+void UDF_Setup0(MPI_Comm comm, setupAide &options)
+{
+  options.getArgs("DENSITY", vtrans);
+  options.getArgs("VISCOSITY", vdiff);
+  options.getArgs("RHOCP", ttrans);
+  options.getArgs("CONDUCTIVITY", tdiff);
+}
+
+int timeStepConverged(nrs_t *nrs, int stage)
+{
+  if (nrs->neknek->nEXT == 1) {
+    return 1;
+  }
+
+  // do one corrector step
+  return stage > 1;
+}
+
+void UDF_Setup(nrs_t *nrs)
+{
+  udf.sEqnSource = &userq;
+  udf.properties = &uservp;
+  udf.timeStepConverged = timeStepConverged;
+}
+
+void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
+{
+
+/*
+  auto *mesh = nrs->meshV;
+
+  const auto LinfUx = platform->linAlg->max(mesh->Nlocal, o_UerrX, platform->comm.mpiCommParent);
+  const auto LinfUy = platform->linAlg->max(mesh->Nlocal, o_UerrY, platform->comm.mpiCommParent);
+  const auto LinfUz = platform->linAlg->max(mesh->Nlocal, o_UerrZ, platform->comm.mpiCommParent);
+
+  if (platform->comm.mpiRank == 0) {
+    printf("LinfUx = %g, LinfUy = %g, LinfUz = %g\n", LinfUx, LinfUy, LinfUz);
+  }
+*/
+  // TODO: add Dillon's print_limits
+
+  if (nrs->isOutputStep) {
+    nek::ocopyToNek(time, tstep);
+    nek::userchk();
+  }
+}