Eikonal solver in distance
This commit is contained in:
JamesEMcclure
parent
4fd74a78a7
commit
d89dcf2648
@ -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>& );
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user