OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixing some compile warnings

Browse Source
This commit is contained in:
Mark Berrill
2017-09-14 08:40:35 -04:00
parent adbd8d2937
commit 182a661fdf
3 changed files with 153 additions and 152 deletions
Download Patch File Download Diff File Expand all files Collapse all files

266
analysis/eikonal.hpp
View File

@@ -9,12 +9,12 @@
inline float minmod(float &a, float &b)
{
float value = a;
if ( a*b < 0.0)
value=0.0;
else if (fabs(a) > fabs(b))
float value = a;
if ( a*b < 0.0)
value=0.0;
else if (fabs(a) > fabs(b))
value = b;
return value;
return value;
}
@@ -23,79 +23,79 @@ inline float minmod(float &a, float &b)
******************************************************************/
inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps)
{
PROFILE_START("Eikonal3D");
PROFILE_START("Eikonal3D");
/*
* 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
*/
/*
* 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;
float dt=0.1;
float Dx,Dy,Dz;
float Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
float Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
float sign,norm;
float LocalVar,GlobalVar,LocalMax,GlobalMax;
int i,j,k;
float dt=0.1;
float Dx,Dy,Dz;
float Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
float Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
float sign,norm;
float LocalVar,GlobalVar,LocalMax,GlobalMax;
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<float> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<float> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
// Arrays to store the second derivatives
Array<float> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
// Arrays to store the second derivatives
Array<float> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
Array<float> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
int count = 0;
while (count < timesteps){
int count = 0;
while (count < timesteps){
// Communicate the halo of values
fillData.fill(Distance);
// 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);
// 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
LocalMax=LocalVar=0.0;
// Execute the next timestep
// f(n+1) = f(n) + dt*sign(1-|grad f|)
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
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;
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
sign = -1;
if (ID(i,j,k) == 1) sign = 1;
sign = -1;
if (ID(i,j,k) == 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));
// 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 ...................
/* //............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;
@@ -103,64 +103,64 @@ inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Dom
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);
Dyp = Distance(i,j+1,k);
Dzp = Distance(i,j,k+1);
*/
Dxp = Distance(i+1,j,k);
Dyp = Distance(i,j+1,k);
Dzp = Distance(i,j,k+1);
Dxm = Distance(i-1,j,k);
Dym = Distance(i,j-1,k);
Dzm = Distance(i,j,k-1);
Dxm = Distance(i-1,j,k);
Dym = Distance(i,j-1,k);
Dzm = Distance(i,j,k-1);
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){
if (Dxp > Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){
if (Dxp > Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
else{
if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
else{
if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
if (norm > 1.0) norm=1.0;
norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
if (norm > 1.0) norm=1.0;
Distance(i,j,k) += dt*sign*(1.0 - norm);
LocalVar += dt*sign*(1.0 - norm);
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));
}
}
}
if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
LocalMax = fabs(dt*sign*(1.0 - norm));
}
}
}
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_FLOAT,MPI_SUM,Dm.Comm);
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_FLOAT,MPI_MAX,Dm.Comm);
GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz;
count++;
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_FLOAT,MPI_SUM,Dm.Comm);
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_FLOAT,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 (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;
}
}
PROFILE_STOP("Eikonal3D");
return GlobalVar;
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf(" Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
}
PROFILE_STOP("Eikonal3D");
return GlobalVar;
}
@@ -210,7 +210,7 @@ inline bool CalcDist3DIteration( Array<float> &Distance, const Domain &Dm )
}
inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
{
PROFILE_START("Calc Distance");
PROFILE_START("Calc Distance");
// Initialize the distance to be 0 fore the cells adjacent to the interface
Distance.fill(1e100);
for (size_t k=1; k<ID.size(2)-1; k++) {
@@ -223,10 +223,10 @@ inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Dom
}
}
// Compute the distance everywhere
fillHalo<float> fillData(Dm.Comm, Dm.rank_info,Dm.Nx,Dm.Ny,Dm.Nz,1,1,1,0,1);
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
fillHalo<float> fillData(Dm.Comm, Dm.rank_info,Dm.Nx,Dm.Ny,Dm.Nz,1,1,1,0,1);
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
// The distance of the cell is the minimum of the distance of the neighbors plus the distance to that node
bool changed = CalcDist3DIteration( Distance, Dm );
changed = sumReduce(Dm.Comm,changed);
@@ -236,7 +236,7 @@ inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Dom
// Update the sign of the distance
for (size_t i=0; i<ID.length(); i++)
Distance(i) *= ID(i)>0 ? 1:-1;
PROFILE_STOP("Calc Distance");
PROFILE_STOP("Calc Distance");
}
@@ -268,29 +268,29 @@ inline void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
ASSERT(int(Distance.size(0))==Nx+2&&int(Distance.size(1))==Ny+2&&int(Distance.size(2))==Nz+2);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
if ( Nx%ratio==0 && Nx>8 && Ny%ratio==0 && Ny>8 && Nz%ratio==0 && Nz>8 ) {
// Use recursive version
int Nr = std::max(std::max(ratio,ratio),ratio);
// Run Nr iterations, communicate, run Nr iterations
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
/*fillData.fill(Distance);
/*fillData.fill(Distance);
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );*/
// Coarsen
Array<float> dist(Nx,Ny,Nz);
fillData.copy(Distance,dist);
fillData.copy(Distance,dist);
Domain Dm2(Nx/ratio,Ny/ratio,Nz/ratio,Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz,Dm.Lx,Dm.Ly,Dm.Lz,0);
Dm2.CommInit(Dm.Comm);
fillHalo<float> fillData2(Dm2.Comm,Dm2.rank_info,Nx/ratio,Ny/ratio,Nz/ratio,1,1,1,0,1);
Dm2.CommInit(Dm.Comm);
fillHalo<float> fillData2(Dm2.Comm,Dm2.rank_info,Nx/ratio,Ny/ratio,Nz/ratio,1,1,1,0,1);
auto dist2 = dist.coarsen( {ratio,ratio,ratio}, coarsen );
Array<float> Distance2(Nx/ratio+2,Ny/ratio+2,Nz/ratio+2);
fillData2.copy(dist2,Distance2);
fillData2.copy(dist2,Distance2);
// Solve
CalcDistMultiLevelHelper( Distance2, Dm2 );
// Interpolate the coarse grid to the fine grid
fillData2.copy(Distance2,dist2);
fillData2.copy(Distance2,dist2);
for (int k=0; k<Nz; k++) {
int k2 = k/ratio;
float z = (k-k2*ratio)-0.5*(ratio-1);
@@ -304,18 +304,18 @@ inline void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
}
}
}
fillData.copy(dist,Distance);
fillData.copy(dist,Distance);
// Run Nr iterations, communicate, run Nr iterations
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
fillData.fill(Distance);
fillData.fill(Distance);
for (int i=0; i<Nr; i++)
CalcDist3DIteration( Distance, Dm );
} else {
// Use coarse-grid version
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
while ( true ) {
// Communicate the halo of values
fillData.fill(Distance);
// The distance of the cell is the minimum of the distance of the neighbors plus the distance to that node
bool changed = CalcDist3DIteration( Distance, Dm );
changed = sumReduce(Dm.Comm,changed);
@@ -326,12 +326,12 @@ inline void CalcDistMultiLevelHelper( Array<float> &Distance, const Domain &Dm )
}
inline void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm )
{
PROFILE_START("Calc Distance Multilevel");
PROFILE_START("Calc Distance Multilevel");
int Nx = Dm.Nx-2;
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
ASSERT(int(Distance.size(0))==Nx+2&&int(Distance.size(1))==Ny+2&&int(Distance.size(2))==Nz+2);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
fillHalo<float> fillData(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
// Initialize the distance to be 0 fore the cells adjacent to the interface
Distance.fill(1e100);
for (size_t k=1; k<ID.size(2)-1; k++) {
@@ -348,13 +348,13 @@ inline void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, c
// Update the sign of the distance
for (size_t i=0; i<ID.length(); i++)
Distance(i) *= ID(i)>0 ? 1:-1;
fillData.fill(Distance);
fillData.fill(Distance);
// Run a quick filter to smooth the data
float sigma = 0.6;
Array<float> H = imfilter::create_filter<float>( { 1 }, "gaussian", &sigma );
std::vector<imfilter::BC> BC(3,imfilter::BC::replicate);
Distance = imfilter::imfilter_separable<float>( Distance, {H,H,H}, BC );
PROFILE_STOP("Calc Distance Multilevel");
PROFILE_STOP("Calc Distance Multilevel");
}
#endif