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

---
 analysis/distance.cpp | 11 ++++++-----
 analysis/distance.h   | 13 ++++++++++++-
 2 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/analysis/distance.cpp b/analysis/distance.cpp
index 3b7641d6..be4cbac6 100644
--- a/analysis/distance.cpp
+++ b/analysis/distance.cpp
@@ -182,7 +182,7 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
     }
 }
 
-double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
+double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps, const std::array<bool,3>& periodic){
 
   /*
    * This routine converts the data in the Distance array to a signed distance
@@ -207,7 +207,8 @@ double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
   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);
+  //fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
+  fillHalo<double> fillData(  Dm.Comm, Dm.rank_info, {xdim, ydim, zdim}, {1,1,1}, 50, 1, {true,true,true}, periodic );
 
   // Arrays to store the second derivatives
   DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
@@ -310,14 +311,14 @@ double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
 
     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;
+    GlobalVar /= Dm.Volume;
     count++;
 
-    if (count%50 == 0 && Dm.rank==0 )
+    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");
+      if (Dm.rank()==0) printf("Exiting with max tolerance of 1e-5 \n");
       count=timesteps;
     }
   }
diff --git a/analysis/distance.h b/analysis/distance.h
index 5ca00c83..b382a49f 100644
--- a/analysis/distance.h
+++ b/analysis/distance.h
@@ -16,6 +16,16 @@ struct Vec {
 };
 inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.z < r.x*r.x+r.y*r.y+r.z*r.z; }
 
+inline double minmod(double &a, double &b){
+
+  double value;
+
+  value = a;
+  if ( a*b < 0.0)    value=0.0;
+  else if (fabs(a) > fabs(b)) value = b;
+
+  return value;
+}
 
 /*!
  * @brief  Calculate the distance using a simple method
@@ -48,7 +58,8 @@ void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
  * @param[in] ID            Domain id
  * @param[in] Dm            Domain information
  * @param[in] timesteps      number of timesteps to run for Eikonal solver
+ * @param[in] periodic      Directions that are periodic 
  */
-double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps);
+double Eikonal(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps, const std::array<bool,3>& periodic);
 
 #endif