From 38ceb73cfe74c3ad382ea16351fee2938894618a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 24 Jan 2018 16:19:10 -0500 Subject: [PATCH] Moved double precision Eikonal() function into eikonal.h / .hpp --- analysis/eikonal.h | 2 + analysis/eikonal.hpp | 157 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 159 insertions(+) diff --git a/analysis/eikonal.h b/analysis/eikonal.h index 8f4545d8..6cc0896c 100644 --- a/analysis/eikonal.h +++ b/analysis/eikonal.h @@ -21,6 +21,8 @@ * @param[in] timesteps Maximum number of timesteps to process * @return Returns the global variation */ +inline double Eikonal(DoubleArray &Distance, const char *ID, const Domain &Dm, int timesteps); + float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps); diff --git a/analysis/eikonal.hpp b/analysis/eikonal.hpp index 8e960ae2..119677d3 100644 --- a/analysis/eikonal.hpp +++ b/analysis/eikonal.hpp @@ -18,9 +18,166 @@ inline float minmod(float &a, float &b) } +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; +} + + /****************************************************************** * Solve the eikonal equation * ******************************************************************/ + + +inline 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 Physics 229 + */ + + 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 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 0.f) Dx = Dxp*Dxp; + else Dx = Dxm*Dxm; + + if (Dyp + Dym > 0.f) Dy = Dyp*Dyp; + else Dy = Dym*Dym; + + if (Dzp + Dzm > 0.f) Dz = Dzp*Dzp; + else Dz = Dzm*Dzm; + } + else{ + + if (Dxp + Dxm < 0.f) Dx = Dxp*Dxp; + else Dx = Dxm*Dxm; + + if (Dyp + Dym < 0.f) Dy = Dyp*Dyp; + else Dy = Dym*Dym; + + if (Dzp + Dzm < 0.f) Dz = Dzp*Dzp; + else Dz = 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; +} + inline float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps) { PROFILE_START("Eikonal3D");