Updaing NekRS 11k case

11145_pebbles/NEKMOOSE

+      parameter (lsurf_m=1700000)
+      common /point_cloudx/ pc_x(lsurf_m*4)
+      common /point_cloudy/ pc_y(lsurf_m*4)
+      common /point_cloudz/ pc_z(lsurf_m*4)
+      common /point_cloudt/ pc_t(lsurf_m*4)
+      common /point_cloudf/ pc_f(lsurf_m*4)
+      common /ref_element/  pc_flag(lsurf_m*4)
+      common /tot_surf/ nw_bdt       
+      common /data_reference/ plmask(lelt), pflmask(lelt), 
+     & pplist(lelt), offset, nw_bd1
+      common /test_passing/ flux_moose, temp_nek
+      real pc_x, pc_y, pc_z, pc_flag, pc_t, pc_f
+      integer plmask, pflmask, pllist, offset, nw_bd1, nw_bdt 
+

11145_pebbles/nek_moose.f

+      subroutine nek_pointscloud()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'
+      parameter(lp=50000)
+      integer e,f
+      integer tots_np(lp), iwk(lp)
+
+      nxyz=lx1*ly1*lz1
+
+      if (np.gt.lp) then
+        write(6,*)"> Increasing lp to: ",np
+        call exitt()
+      endif
+
+      ! Local element arrays containing surface information 
+      nw_bd1 = 0
+      do e=1,nelt
+      eg = lglel(e)
+      pplist(e) = 0
+      pflmask(e) = 0
+      plmask(e) = 0 
+      do f=1,6
+        if (cbc(f,e,2).eq.'f  ') then
+           nw_bd1=nw_bd1+1
+           plmask(e) =1
+           pflmask(e)=f
+           pplist(e)=nw_bd1
+        endif
+      enddo
+      enddo
+
+      ! Local TOTAL number of surfaces per processor  
+      do i=1,lp
+         tots_np(i)=0.0
+      enddo
+      tots_np(nid+1)=nw_bd1
+      call igop(tots_np,iwk,'+  ', np)
+
+      offset=0
+      if (nid.gt.0) then
+        do i=1,nid
+        offset=offset+tots_np(i)
+        enddo
+      endif
+
+      ! Adding local offset 
+      do e=1,nelt
+         pplist(e)=pplist(e)+offset
+      enddo
+
+      nw_bdt=0
+      do i=1,np
+        nw_bdt = nw_bdt + tots_np(i)
+      enddo
+
+      if (nw_bdt.gt.lsurf_m) then
+        write(6,*)"> Increasing lsurf_m to:",nw_bdt
+        call exitt()  
+      endif
+
+      call assign_p(pc_x,xm1)
+      call assign_p(pc_y,ym1)
+      call assign_p(pc_z,zm1)
+
+      do i=1,nw_bdt
+      do j=1,4
+         pc_flag((i-1)*4+j)=dble(i)
+      enddo
+      enddo
+
+c      Uncomment for testing
+c
+c      write(6,*)"(nid, nw_bdt):",nid,nw_bdt
+c
+c      if (nid.eq.0) then
+c      do i=1,nw_bdt
+c      do j=1,4
+c      k=(i-1)*4+j
+c      write(6,*)i,j," - (x,y,z):",pc_x(k),pc_y(k),pc_z(k)
+c      enddo
+c      enddo
+c      endif
+c      call exitt();
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine assign_p (vpf,fu)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'
+      integer e,eg
+      real*8 vpf(lsurf_m*4)
+      real*8 fu(lx1,lx1,lx1,lelt)
+
+      integer i_us, ip_stp, iv_istp
+
+      real*8 wtmp(lsurf_m*4)
+
+      do i=1,4*lsurf_m
+      vpf(i)=0.0
+      enddo
+
+      do e=1,nelt
+
+      if (plmask(e).eq.1)  then
+
+      i_us=(pplist(e)-1)*4
+
+      if (pflmask(e).eq.1) then
+      vpf(i_us+1)=fu(1,1,1,e)
+      vpf(i_us+2)=fu(lx1,1,1,e)
+      vpf(i_us+4)=fu(1,1,lx1,e)
+      vpf(i_us+3)=fu(lx1,1,lx1,e)
+      endif
+
+      if (pflmask(e).eq.2) then
+      vpf(i_us+1)=fu(lx1,1,1,e)
+      vpf(i_us+2)=fu(lx1,lx1,1,e)
+      vpf(i_us+4)=fu(lx1,1,lx1,e)
+      vpf(i_us+3)=fu(lx1,lx1,lx1,e)
+      endif
+
+      if (pflmask(e).eq.3)  then
+      vpf(i_us+1)=fu(1,lx1,1,e)
+      vpf(i_us+2)=fu(lx1,lx1,1,e)
+      vpf(i_us+4)=fu(1,lx1,lx1,e)
+      vpf(i_us+3)=fu(lx1,lx1,lx1,e)
+      endif
+
+      if (pflmask(e).eq.4) then
+      vpf(i_us+1)=fu(1,1,1,e)
+      vpf(i_us+2)=fu(1,lx1,1,e)
+      vpf(i_us+4)=fu(1,1,lx1,e)
+      vpf(i_us+3)=fu(1,lx1,lx1,e)
+      endif
+
+      if (pflmask(e).eq.5) then
+      vpf(i_us+1)=fu(1,1,1,e)
+      vpf(i_us+2)=fu(lx1,1,1,e)
+      vpf(i_us+4)=fu(1,lx1,1,e)
+      vpf(i_us+3)=fu(lx1,lx1,1,e)
+      endif
+
+      if (pflmask(e).eq.6) then
+      vpf(i_us+1)=fu(1,1,lx1,e)
+      vpf(i_us+2)=fu(lx1,1,lx1,e)
+      vpf(i_us+4)=fu(1,lx1,lx1,e)
+      vpf(i_us+3)=fu(lx1,lx1,lx1,e)
+      endif
+
+      endif
+      enddo
+
+      call  gop(vpf,wtmp,'+  ',4*nw_bdt)
+
+      return
+      end
+C=======================================================================
+      subroutine nek_interpolation()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'
+      integer e,f
+
+      call assign_p(pc_t,t)
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine rassign (fu,vpf)
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'
+      integer e,eg
+
+      real*8 vpf(lsurf_m*4)
+      real*8 fu(3,3,3,lelt)
+
+c     Filling a lx1=3 brick for each face from MOOSE. Interpolating 
+c     mid-size points at the moment. 
+      
+      do e=1,nelt
+
+      do i=1,27
+         fu(i,1,1,e)=0.0
+      enddo 
+
+      if (plmask(e).eq.1)  then
+
+      i_us=(pplist(e)-1)*4
+
+      if (pflmask(e).eq.1) then
+      fu(1,1,1,e) = vpf(i_us+1)
+      fu(3,1,1,e) = vpf(i_us+2)
+      fu(1,1,3,e) = vpf(i_us+4)
+      fu(3,1,3,e) = vpf(i_us+3)   
+      fu(2,1,1,e) = 0.5*(fu(1,1,1,e)+fu(3,1,1,e)) ! interpolate midpoints 
+      fu(3,1,2,e) = 0.5*(fu(3,1,1,e)+fu(3,1,3,e))  
+      fu(2,1,3,e) = 0.5*(fu(1,1,3,e)+fu(3,1,3,e))
+      fu(1,1,2,e) = 0.5*(fu(1,1,1,e)+fu(1,1,3,e))
+      fu(2,1,2,e) = 
+     &  0.25*(fu(2,1,1,e)+fu(3,1,2,e)+fu(2,1,3,e)+fu(1,1,2,e)) 
+      do i=1,3         ! fill entire brick with surface solution to enable spectral interpolation
+        do j=1,3
+           fu(i,2,j,e)=fu(i,1,j,e)
+           fu(i,3,j,e)=fu(i,1,j,e) 
+        enddo
+      enddo 
+      endif
+
+      if (pflmask(e).eq.2) then
+      fu(3,1,1,e) = vpf(i_us+1)
+      fu(3,3,1,e) = vpf(i_us+2)
+      fu(3,1,3,e) = vpf(i_us+4)
+      fu(3,3,3,e) = vpf(i_us+3)
+      fu(3,2,1,e) = 0.5*(fu(3,1,1,e)+fu(3,3,1,e)) ! interpolate midpoints 
+      fu(3,3,2,e) = 0.5*(fu(3,3,1,e)+fu(3,3,3,e))
+      fu(3,2,3,e) = 0.5*(fu(3,1,3,e)+fu(3,3,3,e))
+      fu(3,1,2,e) = 0.5*(fu(3,1,1,e)+fu(3,1,3,e))
+      fu(3,2,2,e) =
+     &  0.25*(fu(3,2,1,e)+fu(3,3,2,e)+fu(3,2,3,e)+fu(3,1,2,e))
+      do i=1,3
+        do j=1,3        ! fill entire brick with surface solution to enable spectral interpolation
+  
+           fu(2,i,j,e)=fu(3,i,j,e)
+           fu(1,i,j,e)=fu(3,i,j,e)
+        enddo
+      enddo
+      endif
+
+      if (pflmask(e).eq.3)  then
+      fu(1,3,1,e) = vpf(i_us+1)
+      fu(3,3,1,e) = vpf(i_us+2)
+      fu(1,3,3,e) = vpf(i_us+4)
+      fu(3,3,3,e) = vpf(i_us+3)
+      fu(2,3,1,e) = 0.5*(fu(1,3,1,e)+fu(3,3,1,e)) ! interpolate midpoints 
+      fu(3,3,2,e) = 0.5*(fu(3,3,1,e)+fu(3,3,3,e))
+      fu(2,3,3,e) = 0.5*(fu(1,3,3,e)+fu(3,3,3,e))
+      fu(1,3,2,e) = 0.5*(fu(1,3,1,e)+fu(1,3,3,e))
+      fu(2,3,2,e) =
+     &  0.25*(fu(2,3,1,e)+fu(3,3,2,e)+fu(2,3,3,e)+fu(1,3,2,e))
+      do i=1,3         ! fill entire brick with surface solution to enable spectral interpolation
+        do j=1,3
+           fu(i,1,j,e)=fu(i,3,j,e)
+           fu(i,2,j,e)=fu(i,3,j,e)
+        enddo
+      enddo 
+      endif
+
+      if (pflmask(e).eq.4) then
+      fu(1,1,1,e) = vpf(i_us+1)
+      fu(1,3,1,e) = vpf(i_us+2)
+      fu(1,1,3,e) = vpf(i_us+4)
+      fu(1,3,3,e) = vpf(i_us+3)
+      fu(1,2,1,e) = 0.5*(fu(1,1,1,e)+fu(1,3,1,e)) ! interpolate midpoints 
+      fu(1,3,2,e) = 0.5*(fu(1,3,1,e)+fu(1,3,3,e))
+      fu(1,2,3,e) = 0.5*(fu(1,1,3,e)+fu(1,3,3,e))
+      fu(1,1,2,e) = 0.5*(fu(1,1,1,e)+fu(1,1,3,e))
+      fu(1,2,2,e) =
+     &  0.25*(fu(1,2,1,e)+fu(1,3,2,e)+fu(1,2,3,e)+fu(1,1,2,e))
+      do i=1,3
+        do j=1,3        ! fill entire brick with surface solution to enable spectral interpolation
+           fu(2,i,j,e)=fu(1,i,j,e)
+           fu(3,i,j,e)=fu(1,i,j,e)
+        enddo
+      enddo
+      endif
+
+      if (pflmask(e).eq.5) then
+      fu(1,1,1,e) = vpf(i_us+1)
+      fu(3,1,1,e) = vpf(i_us+2)
+      fu(1,3,1,e) = vpf(i_us+4)
+      fu(3,3,1,e) = vpf(i_us+3)
+      fu(2,1,1,e) = 0.5*(fu(1,1,1,e)+fu(3,1,1,e)) ! interpolate midpoints 
+      fu(3,2,1,e) = 0.5*(fu(3,1,1,e)+fu(3,3,1,e))
+      fu(2,3,1,e) = 0.5*(fu(1,3,1,e)+fu(3,3,1,e))
+      fu(1,2,1,e) = 0.5*(fu(1,1,1,e)+fu(1,3,1,e))
+      fu(2,2,1,e) =
+     &  0.25*(fu(2,1,1,e)+fu(3,2,1,e)+fu(2,3,1,e)+fu(1,2,1,e))
+      do i=1,3
+        do j=1,3        ! fill entire brick with surface solution to enable spectral interpolation
+           fu(i,j,2,e)=fu(i,j,1,e)
+           fu(i,j,3,e)=fu(i,j,1,e)
+        enddo
+      enddo
+      endif
+
+      if (pflmask(e).eq.6) then
+      fu(1,1,3,e) = vpf(i_us+1)
+      fu(3,1,3,e) = vpf(i_us+2)
+      fu(1,3,3,e) = vpf(i_us+4)
+      fu(3,3,3,e) = vpf(i_us+3)
+      fu(2,1,3,e) = 0.5*(fu(1,1,3,e)+fu(3,1,3,e)) ! interpolate midpoints 
+      fu(3,2,3,e) = 0.5*(fu(3,1,3,e)+fu(3,3,3,e))
+      fu(2,3,3,e) = 0.5*(fu(1,3,3,e)+fu(3,3,3,e))
+      fu(1,2,3,e) = 0.5*(fu(1,1,3,e)+fu(1,3,3,e))
+      fu(2,2,3,e) =
+     &  0.25*(fu(2,1,3,e)+fu(3,2,3,e)+fu(2,3,3,e)+fu(1,2,3,e))
+      do i=1,3
+        do j=1,3        ! fill entire brick with surface solution to enable spectral interpolation
+           fu(i,j,1,e)=fu(i,j,3,e)
+           fu(i,j,2,e)=fu(i,j,3,e)
+        enddo
+      enddo
+      endif
+
+      endif
+      enddo  
+
+      return
+      end
+C=======================================================================
+      subroutine int_surf12(fo,fi,n2,n1,ne)
+      include 'SIZE'
+
+      parameter(lx=20)
+      real*8 z1(lx),z2(lx),w(lx)
+      real*8 fm12(lx*lx),fm12t(lx*lx)
+      real*8 work(lx*lx*lx)
+      real*8 fo(n2,n2,n2,lelt), fi(n1,n1,n1,lelt)
+      integer n1, n2, e, ne
+
+      if (n1.gt.lx.or.n2.gt.lx) then
+         write(6,*)'ERROR: increase lx in setmap to max:',n1,n2
+         call exitt
+      endif
+
+      call rzero(fo,n2*n2*n2*ne)
+
+      call zwgll(z1,w,n1)
+      call zwgll(z2,w,n2)
+      call igllm(fm12,fm12t,z1,z2,n1,n2,n1,n2)
+
+      do e=1,ne
+      call specmpf(fo(1,1,1,e),n2,fi(1,1,1,e),n1,fm12,fm12t,work)
+      enddo
+
+      return
+      end
+c-----------------------------------------------------------------------
+      subroutine specmpf(b,nb,a,na,ba,ab,w)
+C
+C     -  Spectral interpolation from A to B via tensor products
+C     -  scratch arrays: w(na*na*nb)
+C
+C
+      include 'SIZE'
+      include 'INPUT'
+      real b(nb,nb,nb),a(na,na,na)
+      real w(1)
+C
+
+         nab = na*nb
+         nbb = nb*nb
+         call mxm(ba,nb,a,na,b,na*na)
+         k=1
+         l=1
+         do iz=1,na
+            call mxm(b(k,1,1),nb,ab,na,w(l),nb)
+            k=k+nab
+            l=l+nbb
+         enddo
+         call mxm(w,nbb,ab,na,b,nb)
+      return
+      end
+C=======================================================================
+      subroutine flux_reconstruction()
+      include 'SIZE'
+      include 'TOTAL'
+      include 'NEKMOOSE'
+      common /flux_reconstructed/flux_recon(lx1,ly1,lz1,lelt)
+      integer e,f,lx_low
+      real*8  flux_low(3,3,3,lelt) 
+      real*8  flux_recon
+      ntot=nx1*ny1*nz1*nelt
+!     Call to rassign and spectral interpolation
+
+!     Testing on pc_f
+      ntot=nw_bdt*4
+      pc_f_max=glmax(pc_f,ntot)
+      pc_f_min=glmin(pc_f,ntot) 
+
+      call rassign(flux_low,pc_f)
+!     Testing on flux_low
+      ntot=nelt*3*3*3
+      fl_max=glmax(flux_low,ntot)
+      fl_min=glmin(flux_low,ntot)
+
+      lx_low=3
+      call int_surf12(flux_recon,flux_low,lx1,lx_low,lelt) 
+
+      ntot=nelt*lx1*lx1*lx1
+      fh_max=glmax(flux_recon,ntot)
+      fh_min=glmin(flux_recon,ntot)
+
+      if (nid.eq.0) then
+        write(6,*)"> Min/Max flux 1: ",pc_f_min,pc_f_max
+        write(6,*)"> Min/Max flux 2: ",fl_min,fl_max
+        write(6,*)"> Min/Max flux 3: ",fh_min,fh_max
+      endif
+
+
+      sint1=0.0
+      sarea1=0.0
+
+      do e=1,nelt
+        do f=1,6
+          call surface_int(sint,sarea,flux_recon,e,f)
+          if (cbc(f,e,1).eq.'W  ') then
+           sint1=sint1+sint
+           sarea1=sarea1+sarea
+          endif
+         enddo
+      enddo
+
+      call  gop(sint1,wtmp,'+  ',1)
+      call  gop(sarea1,wtmp,'+  ',1)
+
+      flux_moose1=sint1
+
+      if (nid.eq.0) then
+        write(6,*)"> Average reconstructed flux: ",flux_moose1
+        write(6,*)"> Imposing flux: ",flux_moose
+      endif
+
+      ntot  = lx1*ly1*lz1*nelt
+      scale=1.0
+      if (abs(flux_moose1).gt.0.0) scale = flux_moose/flux_moose1 
+      do i=1,ntot  
+         flux_recon(i,1,1,1)=flux_recon(i,1,1,1)*scale
+      enddo
+
+c      call exitt()  
+
+      return
+      end
+c-----------------------------------------------------------------------

11145_pebbles/peb11k.oudf

 // Boundary conditions
-/* wall 1, inflow 2, outflow 3, x-slip 4, y-slip 5, z-slip 6 */
-void insVelocityDirichletConditions3D(bcData *bc)
+void velocityDirichletConditions(bcData *bc)
 {                                                                        
-  bc->uP = 0.0;
-  bc->vP = 0.0;
-  bc->wP = 1.0;
+  bc->u = 0.0;
+  bc->v = 0.0;
+  bc->w = 0.17;
 }
 
-void insVelocityNeumannConditions3D(bcData *bc)
+void scalarDirichletConditions(bcData *bc)
 {
+  bc->s = 573.0;
 }
 
-void cdsDirichletConditions3D(bcData *bc)
+void scalarNeumannConditions(bcData *bc)
 {
-  bc->sP = 0.0;
-}
-void cdsNeumannConditions3D(bcData *bc)
-{
-  bc->sF = 1.0;
+  bc->flux = bc->wrk[bc->idM];
 }
 
 // Stabilized outflow (Dong et al)
-void insPressureDirichletConditions3D(bcData *bc)
+void pressureDirichletConditions(bcData *bc)
 {
   const dfloat iU0delta = 10.0;
-  const dfloat un = bc->uM*bc->nx + bc->vM*bc->ny + bc->wM*bc->nz;
-  const dfloat s0 = 0.5 * (1.0 - tanh(un*iU0delta)); 
-  bc->pP = -0.5 * (bc->uM*bc->uM + bc->vM*bc->vM + bc->wM*bc->wM) * s0;
+  const dfloat un = bc->u*bc->nx + bc->v*bc->ny + bc->w*bc->nz;
+  const dfloat s0 = 0.5 * (1.0 - tanh(un*iU0delta));
+  bc->p = -0.5 * (bc->u*bc->u + bc->v*bc->v + bc->w*bc->w) * s0;
 }
+

11145_pebbles/peb11k.udf

@@ -4,31 +4,35 @@
 
 #include <math.h>
 #include "udf.hpp"
-#include "cfl.hpp"
-
-int my_iostep=100; 
+#include <iostream>
 
 /* UDF Functions */                                                      
 
-void UDF_LoadKernels(ins_t *ins)
+void UDF_LoadKernels(nrs_t *nrs)
 {
 }
 
-void UDF_Setup(ins_t *ins)
+void UDF_Setup(nrs_t *nrs)
 {
   // get IC from nek
-  if (!ins->readRestartFile) nek_copyTo(ins->startTime);
+  nrs->o_usrwrk.free(); free(nrs->usrwrk);
+  nrs->usrwrk   = (dfloat*) calloc(nrs->mesh->Nelements*nrs->mesh->Np, sizeof(dfloat));
+  nrs->o_usrwrk = nrs->mesh->device.malloc(nrs->mesh->Nelements*nrs->mesh->Np*sizeof(dfloat), nrs->usrwrk);
+  // nrs->cds->o_usrwrk = nrs->o_usrwrk;
+  // printf ("> Test %d \n", nrs->mesh->Nelements*nrs->mesh->Np*sizeof(dfloat));
 }
 
-void UDF_ExecuteStep(ins_t *ins, dfloat time, int tstep)
+void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
 {
-  dfloat cfl=computeCFL(ins, time, tstep);
-  if (ins->isOutputStep || cfl>4.0) {
-    nek_ocopyFrom(time, tstep);
-    nek_userchk();
-  }
-  if (ins->isOutputStep) {
-    nek_ocopyFrom(time, tstep);
-    nek_userchk();
-  }
+  nek_userchk();
+  //std::cout << "Finished userchk" << std::endl;
+
+  int nwbdt=nekData.cbscnrs[0];
+  auto *wrk1 =  nekData.cbscnrs; 
+  wrk1 = wrk1 + nwbdt*5*4 + 2; 
+  nrs->o_usrwrk.copyFrom(wrk1,nrs->mesh->Nelements*nrs->mesh->Np*sizeof(dfloat));
+
+  //auto *ptr=nekData.cbscnrs;
+  //nrs->o_usrwrk.copyFrom(wrk, mesh->Nelements*mesh->Np*sizeof(dfloat));
+  //std::cout << "Finished UDF_ExecuteStep" << std::endl;
 }

11145_pebbles/peb11k.usr

+      include 'nek_moose.f'
 c-----------------------------------------------------------------------
       subroutine uservp(i,j,k,eg) ! set variable properties
       include 'SIZE'
@@ -84,7 +85,7 @@ c-----------------------------------------------------------------------
       ux   = 0.0
       uy   = 0.0
       uz   = 1.0
-      temp = 0.0
+      temp = 573.0
 
       return
       end
@@ -92,10 +93,21 @@ c-----------------------------------------------------------------------
       subroutine userchk
       include 'SIZE'
       include 'TOTAL'
+      include 'NEKMOOSE'
+c      implicit none
+      COMMON /SCNRS/ SC_NRS(LX1*LY1*LZ1*LELT*7*5)
+      real           SC_NRS
+      common /flux_reconstructed/flux_recon(lx1,ly1,lz1,lelt)
+
+      integer n_count
+      save n_count
 
       n=lx1*ly1*lz1*nelt
 
+      time=0.0
+ 
       if (istep.eq.0) then
+        time=0.0 
         xxmax = glmax(xm1,n)
         yymax = glmax(ym1,n)