From d89dcf2648d0b9e23d786bb6d3f8cc0d380111f3 Mon Sep 17 00:00:00 2001
From: JamesEMcclure <jemcclur@gmail.com>
Date: Sat, 27 Jun 2020 07:28:48 -0400
Subject: [PATCH] Eikonal solver in distance

---
 analysis/distance.cpp | 141 ++++++++++++++++++++++++++++++++++++++++++
 analysis/distance.h   |  11 ++++
 2 files changed, 152 insertions(+)

diff --git a/analysis/distance.cpp b/analysis/distance.cpp
index e297b435..78c344d2 100644
--- a/analysis/distance.cpp
+++ b/analysis/distance.cpp
@@ -182,6 +182,147 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
     }
 }
 
+double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
+
+  /*
+   * This routine converts the data in the Distance array to a signed distance
+   * by solving the equation df/dt = sign(1-|grad f|), where Distance provides
+   * the values of f on the mesh associated with domain Dm
+   * It has been tested with segmented data initialized to values [-1,1]
+   * and will converge toward the signed distance to the surface bounding the associated phases
+   *
+   * Reference:
+   * Min C (2010) On reinitializing level set functions, Journal of Computational Physics229
+   */
+
+  int i,j,k;
+  double dt=0.1;
+  double Dx,Dy,Dz;
+  double Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
+  double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
+  double sign,norm;
+  double LocalVar,GlobalVar,LocalMax,GlobalMax;
+
+  int xdim,ydim,zdim;
+  xdim=Dm.Nx-2;
+  ydim=Dm.Ny-2;
+  zdim=Dm.Nz-2;
+  fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
+
+  // Arrays to store the second derivatives
+  DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
+  DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
+  DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
+
+  int count = 0;
+  while (count < timesteps){
+
+    // Communicate the halo of values
+    fillData.fill(Distance);
+
+    // Compute second order derivatives
+    for (k=1;k<Dm.Nz-1;k++){
+      for (j=1;j<Dm.Ny-1;j++){
+	for (i=1;i<Dm.Nx-1;i++){
+	  Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
+	  Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
+	  Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
+	}
+      }
+    }
+    fillData.fill(Dxx);
+    fillData.fill(Dyy);
+    fillData.fill(Dzz);
+
+    LocalMax=LocalVar=0.0;
+    // Execute the next timestep
+    for (k=1;k<Dm.Nz-1;k++){
+      for (j=1;j<Dm.Ny-1;j++){
+	for (i=1;i<Dm.Nx-1;i++){
+
+	  int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
+
+	  sign = -1;
+	  if (ID[n] == 1) sign = 1;
+
+	  // local second derivative terms
+	  Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
+	  Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
+	  Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
+	  Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
+	  Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
+	  Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
+
+	  /* //............Compute upwind derivatives ...................
+                    Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
+                    Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
+                    Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
+                    Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
+                    Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
+                    Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
+	  */
+	  Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
+	  Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
+	  Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
+
+	  Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
+	  Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
+	  Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
+
+	  // Compute upwind derivatives for Godunov Hamiltonian
+	  if (sign < 0.0){
+	    if (Dxp + Dxm > 0.f)  Dx = Dxp*Dxp;
+	    elseDx = Dxm*Dxm;
+
+	    if (Dyp + Dym > 0.f)  Dy = Dyp*Dyp;
+	    elseDy = Dym*Dym;
+
+	    if (Dzp + Dzm > 0.f)  Dz = Dzp*Dzp;
+	    elseDz = Dzm*Dzm;
+	  }
+	  else{
+
+	    if (Dxp + Dxm < 0.f)  Dx = Dxp*Dxp;
+	    elseDx = Dxm*Dxm;
+
+	    if (Dyp + Dym < 0.f)  Dy = Dyp*Dyp;
+	    elseDy = Dym*Dym;
+
+	    if (Dzp + Dzm < 0.f)  Dz = Dzp*Dzp;
+	    elseDz = Dzm*Dzm;
+	  }
+
+	  //Dx = max(Dxp*Dxp,Dxm*Dxm);
+	  //Dy = max(Dyp*Dyp,Dym*Dym);
+	  //Dz = max(Dzp*Dzp,Dzm*Dzm);
+
+	  norm=sqrt(Dx + Dy + Dz);
+	  if (norm > 1.0) norm=1.0;
+
+	  Distance(i,j,k) += dt*sign*(1.0 - norm);
+	  LocalVar +=  dt*sign*(1.0 - norm);
+
+	  if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
+	    LocalMax = fabs(dt*sign*(1.0 - norm));
+	}
+      }
+    }
+
+    MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
+    MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm);
+    GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz;
+    count++;
+
+    if (count%50 == 0 && Dm.rank==0 )
+      printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
+
+    if (fabs(GlobalMax) < 1e-5){
+      if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
+      count=timesteps;
+    }
+  }
+  return GlobalVar;
+}
 
 // Explicit instantiations
 template void CalcDist<float>( Array<float>&, const Array<char>&, const Domain&, const std::array<bool,3>&, const std::array<double,3>& );
diff --git a/analysis/distance.h b/analysis/distance.h
index b3fc870e..5ca00c83 100644
--- a/analysis/distance.h
+++ b/analysis/distance.h
@@ -40,4 +40,15 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
 void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
     const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );
 
+
+/*!
+ * @brief  Calculate the distance based on solution of Eikonal equation
+ * @details  This routine calculates the signed distance to the nearest domain surface.
+ * @param[out] Distance     Distance function
+ * @param[in] ID            Domain id
+ * @param[in] Dm            Domain information
+ * @param[in] timesteps      number of timesteps to run for Eikonal solver
+ */
+double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps);
+
 #endif