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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2018-09-22 16:00:44 -04:00
commit c2ec00941e
6 changed files with 228 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
analysis/TwoPhase.cpp
View File

@ -548,7 +548,6 @@ void TwoPhase::ComputeLocal()
} }
} }
} }
Array <char> phase_label(Nx,Ny,Nz); Array <char> phase_label(Nx,Ny,Nz);
Array <double> phase_distance(Nx,Ny,Nz); Array <double> phase_distance(Nx,Ny,Nz);
// Analyze the wetting fluid // Analyze the wetting fluid
@ -598,7 +597,6 @@ void TwoPhase::ComputeLocal()
} }
CalcDist(phase_distance,phase_label,*Dm); CalcDist(phase_distance,phase_label,*Dm);
nonwet_morph.ComputeScalar(phase_distance,0.f); nonwet_morph.ComputeScalar(phase_distance,0.f);
} }

229
analysis/dcel.cpp
View File

@ -1,5 +1,78 @@
#include "analysis/dcel.h" #include "analysis/dcel.h"
/*
Double connected edge list (DECL)
*/
Vertex::Vertex(){
size_ = 0;
vertex_data.resize(36);
}
Vertex::~Vertex(){
}
void Vertex::add(Point P){
vertex_data.push_back(P.x);
vertex_data.push_back(P.y);
vertex_data.push_back(P.z);
size_++;
}
void Vertex::assign(int idx, Point P){
vertex_data[3*idx] = P.x;
vertex_data[3*idx+1] = P.y;
vertex_data[3*idx+2] = P.z;
}
int Vertex::size(){
return size_;
}
Point Vertex::coords(int idx){
Point P;
P.x = vertex_data[3*idx];
P.y = vertex_data[3*idx+1];
P.z = vertex_data[3*idx+2];
return P;
}
Halfedge::Halfedge(){
size_=0;
}
Halfedge::~Halfedge(){
}
int Halfedge::v1(int edge){
return data(0,edge);
}
int Halfedge::v2(int edge){
return data(1,edge);
}
int Halfedge::face(int edge){
return data(2,edge);
}
int Halfedge::twin(int edge){
return data(3,edge);
}
int Halfedge::prev(int edge){
return data(4,edge);
}
int Halfedge::next(int edge){
return data(5,edge);
}
int Halfedge::size(){
return size_;
}
DECL::DECL(){ DECL::DECL(){
@ -12,20 +85,24 @@ DECL::~DECL(){
} }
int DECL::Face(int index){ int DECL::Face(int index){
return FaceData[index]; return FaceData(index);
} }
void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, const int j, const int k){ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const int j, const int k){
Point P,Q; Point P,Q;
Point PlaceHolder; Point PlaceHolder;
Point C0,C1,C2,C3,C4,C5,C6,C7; Point C0,C1,C2,C3,C4,C5,C6,C7;
int CubeIndex;
int nTris = 0;
int nVert =0;
Point VertexList[12]; Point VertexList[12];
Point NewVertexList[12]; Point NewVertexList[12];
int LocalRemap[12]; int LocalRemap[12];
Point cellvertices[20]; DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
std::array<std::array<int,3>,20> Triangles; IntArray Triangles = IntArray(3,20);
// Values from array 'A' at the cube corners // Values from array 'A' at the cube corners
double CubeValues[8]; double CubeValues[8];
@ -52,7 +129,7 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
//Determine the index into the edge table which //Determine the index into the edge table which
//tells us which vertices are inside of the surface //tells us which vertices are inside of the surface
int CubeIndex = 0; CubeIndex = 0;
if (CubeValues[0] < 0.0f) CubeIndex |= 1; if (CubeValues[0] < 0.0f) CubeIndex |= 1;
if (CubeValues[1] < 0.0f) CubeIndex |= 2; if (CubeValues[1] < 0.0f) CubeIndex |= 2;
if (CubeValues[2] < 0.0f) CubeIndex |= 4; if (CubeValues[2] < 0.0f) CubeIndex |= 4;
@ -145,30 +222,31 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
//P.x += i; //P.x += i;
//P.y += j; //P.y += j;
//P.z += k; //P.z += k;
cellvertices[idx] = P; cellvertices(idx) = P;
} }
nVert = VertexCount;
TriangleCount = 0; TriangleCount = 0;
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) { for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
Triangles[TriangleCount][0] = LocalRemap[triTable[CubeIndex][idx+0]]; Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
Triangles[TriangleCount][1] = LocalRemap[triTable[CubeIndex][idx+1]]; Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
Triangles[TriangleCount][2] = LocalRemap[triTable[CubeIndex][idx+2]]; Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
TriangleCount++; TriangleCount++;
} }
int nTris = TriangleCount; nTris = TriangleCount;
// Now add the local values to the DECL data structure // Now add the local values to the DECL data structure
if (nTris>0){ if (nTris>0){
FaceData.resize(TriangleCount); FaceData.resize(TriangleCount);
//printf("Construct halfedge structure... \n"); //printf("Construct halfedge structure... \n");
//printf(" Construct %i triangles \n",nTris); //printf(" Construct %i triangles \n",nTris);
halfedge.resize(nTris*3); halfedge.data.resize(6,nTris*3);
int idx_edge=0; int idx_edge=0;
for (int idx=0; idx<TriangleCount; idx++){ for (int idx=0; idx<TriangleCount; idx++){
int V1 = Triangles[idx][0]; int V1 = Triangles(0,idx);
int V2 = Triangles[idx][1]; int V2 = Triangles(1,idx);
int V3 = Triangles[idx][2]; int V3 = Triangles(2,idx);
FaceData[idx] = idx_edge; FaceData(idx) = idx_edge;
// first edge: V1->V2 // first edge: V1->V2
halfedge.data(0,idx_edge) = V1; // first vertex halfedge.data(0,idx_edge) = V1; // first vertex
halfedge.data(1,idx_edge) = V2; // second vertex halfedge.data(1,idx_edge) = V2; // second vertex
@ -212,8 +290,8 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
} }
} }
// Use "ghost" twins if edge is on a cube face // Use "ghost" twins if edge is on a cube face
P = cellvertices[V1]; P = cellvertices(V1);
Q = cellvertices[V2]; Q = cellvertices(V2);
if (P.x == 0.0 && Q.x == 0.0) halfedge.data(3,idx) = -1; // ghost twin for x=0 face if (P.x == 0.0 && Q.x == 0.0) halfedge.data(3,idx) = -1; // ghost twin for x=0 face
if (P.x == 1.0 && Q.x == 1.0) halfedge.data(3,idx) = -4; // ghost twin for x=1 face if (P.x == 1.0 && Q.x == 1.0) halfedge.data(3,idx) = -4; // ghost twin for x=1 face
if (P.y == 0.0 && Q.y == 0.0) halfedge.data(3,idx) = -2; // ghost twin for y=0 face if (P.y == 0.0 && Q.y == 0.0) halfedge.data(3,idx) = -2; // ghost twin for y=0 face
@ -225,7 +303,7 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
// Map vertices to global coordinates // Map vertices to global coordinates
for (int idx=0;idx<VertexCount;idx++) { for (int idx=0;idx<VertexCount;idx++) {
P = cellvertices[idx]; P = cellvertices(idx);
P.x += i; P.x += i;
P.y += j; P.y += j;
P.z += k; P.z += k;
@ -280,7 +358,8 @@ Point DECL::TriNormal(int edge)
double DECL::EdgeAngle(int edge) double DECL::EdgeAngle(int edge)
{ {
double angle; double angle;
double dotprod; double nx,ny,nz;
double dotprod,length,hypotenuse;
Point P,Q,R; // triangle vertices Point P,Q,R; // triangle vertices
Point U,V,W; // normal vectors Point U,V,W; // normal vectors
int e2 = halfedge.next(edge); int e2 = halfedge.next(edge);
@ -293,14 +372,14 @@ double DECL::EdgeAngle(int edge)
if (halfedge.twin(edge) < 0 ){ if (halfedge.twin(edge) < 0 ){
// compute edge normal in plane of cube face // compute edge normal in plane of cube face
W = P - Q; // edge tangent vector W = P - Q; // edge tangent vector
double length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z); length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z);
W.x /= length; W.x /= length;
W.y /= length; W.y /= length;
W.z /= length; W.z /= length;
// edge normal within the plane of the cube face // edge normal within the plane of the cube face
double nx = W.y*V.z - W.z*V.y; nx = W.y*V.z - W.z*V.y;
double ny = W.z*V.x - W.x*V.z; ny = W.z*V.x - W.x*V.z;
double nz = W.x*V.y - W.y*V.x; nz = W.x*V.y - W.y*V.x;
length = sqrt(nx*nx+ny*ny+nz*nz); length = sqrt(nx*nx+ny*ny+nz*nz);
// new value for V is this normal vector // new value for V is this normal vector
V.x = nx/length; V.y = ny/length; V.z = nz/length; V.x = nx/length; V.y = ny/length; V.z = nz/length;
@ -352,19 +431,20 @@ void Isosurface(DoubleArray &A, const double &v)
{ {
Point P,Q; Point P,Q;
Point PlaceHolder; Point PlaceHolder;
double temp;
Point C0,C1,C2,C3,C4,C5,C6,C7; Point C0,C1,C2,C3,C4,C5,C6,C7;
int TriangleCount; int TriangleCount;
int VertexCount; int VertexCount;
int CubeIndex; int CubeIndex;
int nTris, nVert;
Point VertexList[12]; Point VertexList[12];
Point NewVertexList[12]; Point NewVertexList[12];
int LocalRemap[12]; int LocalRemap[12];
Point cellvertices[20]; DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
std::array<std::array<int,3>,20> Triangles; IntArray Triangles = IntArray(3,20);
Triangles.fill( { 0 } );
// Values from array 'A' at the cube corners // Values from array 'A' at the cube corners
double CubeValues[8]; double CubeValues[8];
@ -383,7 +463,6 @@ void Isosurface(DoubleArray &A, const double &v)
C6.x = 1.0; C6.y = 1.0; C6.z = 1.0; C6.x = 1.0; C6.y = 1.0; C6.z = 1.0;
C7.x = 0.0; C7.y = 1.0; C7.z = 1.0; C7.x = 0.0; C7.y = 1.0; C7.z = 1.0;
std::vector<std::array<int,6>> HalfEdge;
for (int k=1; k<Nz-1; k++){ for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){ for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){ for (int i=1; i<Nx-1; i++){
@ -490,81 +569,83 @@ void Isosurface(DoubleArray &A, const double &v)
//P.x += i; //P.x += i;
//P.y += j; //P.y += j;
//P.z += k; //P.z += k;
cellvertices[idx] = P; cellvertices(idx) = P;
} }
nVert = VertexCount;
TriangleCount = 0; TriangleCount = 0;
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) { for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
Triangles[TriangleCount][0] = LocalRemap[triTable[CubeIndex][idx+0]]; Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
Triangles[TriangleCount][1] = LocalRemap[triTable[CubeIndex][idx+1]]; Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
Triangles[TriangleCount][2] = LocalRemap[triTable[CubeIndex][idx+2]]; Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
TriangleCount++; TriangleCount++;
} }
int nTris = TriangleCount; nTris = TriangleCount;
// Now add the local values to the DECL data structure // Now add the local values to the DECL data structure
HalfEdge.resize(nTris*3); IntArray HalfEdge(6,nTris*3);
DoubleArray EdgeAngles(nTris*3);
int idx_edge=0; int idx_edge=0;
for (int idx=0; idx<TriangleCount; idx++){ for (int idx=0; idx<TriangleCount; idx++){
int V1 = Triangles[idx][0]; int V1 = Triangles(0,idx);
int V2 = Triangles[idx][1]; int V2 = Triangles(1,idx);
int V3 = Triangles[idx][2]; int V3 = Triangles(2,idx);
// first edge: V1->V2 // first edge: V1->V2
HalfEdge[idx_edge][0] = V1; // first vertex HalfEdge(0,idx_edge) = V1; // first vertex
HalfEdge[idx_edge][1] = V2; // second vertex HalfEdge(1,idx_edge) = V2; // second vertex
HalfEdge[idx_edge][2] = idx; // triangle HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge[idx_edge][3] = -1; // twin HalfEdge(3,idx_edge) = -1; // twin
HalfEdge[idx_edge][4] = idx_edge+2; // previous edge HalfEdge(4,idx_edge) = idx_edge+2; // previous edge
HalfEdge[idx_edge][5] = idx_edge+1; // next edge HalfEdge(5,idx_edge) = idx_edge+1; // next edge
idx_edge++; idx_edge++;
// second edge: V2->V3 // second edge: V2->V3
HalfEdge[idx_edge][0] = V2; // first vertex HalfEdge(0,idx_edge) = V2; // first vertex
HalfEdge[idx_edge][1] = V3; // second vertex HalfEdge(1,idx_edge) = V3; // second vertex
HalfEdge[idx_edge][2] = idx; // triangle HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge[idx_edge][3] = -1; // twin HalfEdge(3,idx_edge) = -1; // twin
HalfEdge[idx_edge][4] = idx_edge-1; // previous edge HalfEdge(4,idx_edge) = idx_edge-1; // previous edge
HalfEdge[idx_edge][5] = idx_edge+1; // next edge HalfEdge(5,idx_edge) = idx_edge+1; // next edge
idx_edge++; idx_edge++;
// third edge: V3->V1 // third edge: V3->V1
HalfEdge[idx_edge][0] = V3; // first vertex HalfEdge(0,idx_edge) = V3; // first vertex
HalfEdge[idx_edge][1] = V1; // second vertex HalfEdge(1,idx_edge) = V1; // second vertex
HalfEdge[idx_edge][2] = idx; // triangle HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge[idx_edge][3] = -1; // twin HalfEdge(3,idx_edge) = -1; // twin
HalfEdge[idx_edge][4] = idx_edge-1; // previous edge HalfEdge(4,idx_edge) = idx_edge-1; // previous edge
HalfEdge[idx_edge][5] = idx_edge-2; // next edge HalfEdge(5,idx_edge) = idx_edge-2; // next edge
idx_edge++; idx_edge++;
} }
int EdgeCount=idx_edge; int EdgeCount=idx_edge;
for (int idx=0; idx<EdgeCount; idx++){ for (int idx=0; idx<EdgeCount; idx++){
int V1=HalfEdge[idx][0]; int V1=HalfEdge(0,idx);
int V2=HalfEdge[idx][1]; int V2=HalfEdge(1,idx);
// Find all the twins within the cube // Find all the twins within the cube
for (int jdx=0; idx<EdgeCount; jdx++){ for (int jdx=0; idx<EdgeCount; jdx++){
if (HalfEdge[jdx][1] == V1 && HalfEdge[jdx][0] == V2){ if (HalfEdge(1,jdx) == V1 && HalfEdge(0,jdx) == V2){
// this is the pair // this is the pair
HalfEdge[idx][3] = jdx; HalfEdge(3,idx) = jdx;
HalfEdge[jdx][3] = idx; HalfEdge(3,jdx) = idx;
} }
if (HalfEdge[jdx][1] == V2 && HalfEdge[jdx][0] == V1 && !(idx==jdx)){ if (HalfEdge(1,jdx) == V2 && HalfEdge(0,jdx) == V1 && !(idx==jdx)){
std::printf("WARNING: half edges with identical orientation! \n"); std::printf("WARNING: half edges with identical orientation! \n");
} }
} }
// Use "ghost" twins if edge is on a cube face // Use "ghost" twins if edge is on a cube face
P = cellvertices[V1]; P = cellvertices(V1);
Q = cellvertices[V2]; Q = cellvertices(V2);
if (P.x == 0.0 && Q.x == 0.0) HalfEdge[idx_edge][3] = -1; // ghost twin for x=0 face if (P.x == 0.0 && Q.x == 0.0) HalfEdge(3,idx_edge) = -1; // ghost twin for x=0 face
if (P.x == 1.0 && Q.x == 1.0) HalfEdge[idx_edge][3] = -2; // ghost twin for x=1 face if (P.x == 1.0 && Q.x == 1.0) HalfEdge(3,idx_edge) = -2; // ghost twin for x=1 face
if (P.y == 0.0 && Q.y == 0.0) HalfEdge[idx_edge][3] = -3; // ghost twin for y=0 face if (P.y == 0.0 && Q.y == 0.0) HalfEdge(3,idx_edge) = -3; // ghost twin for y=0 face
if (P.y == 1.0 && Q.y == 1.0) HalfEdge[idx_edge][3] = -4; // ghost twin for y=1 face if (P.y == 1.0 && Q.y == 1.0) HalfEdge(3,idx_edge) = -4; // ghost twin for y=1 face
if (P.z == 0.0 && Q.z == 0.0) HalfEdge[idx_edge][3] = -5; // ghost twin for z=0 face if (P.z == 0.0 && Q.z == 0.0) HalfEdge(3,idx_edge) = -5; // ghost twin for z=0 face
if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6; // ghost twin for z=1 face if (P.z == 1.0 && Q.z == 1.0) HalfEdge(3,idx_edge) = -6; // ghost twin for z=1 face
} }
// Find all the angles // Find all the angles
for (int idx=0; idx<EdgeCount; idx++){ for (int idx=0; idx<EdgeCount; idx++){
int V1=HalfEdge[idx][0]; int V1=HalfEdge(0,idx);
int V2=HalfEdge[idx][1]; int V2=HalfEdge(1,idx);
int T1= HalfEdge[idx_edge][2]; int T1= HalfEdge(2,idx_edge);
int twin=HalfEdge[idx_edge][3]; int twin=HalfEdge(3,idx_edge);
if (twin == -1){ if (twin == -1){
} }
@ -572,11 +653,11 @@ void Isosurface(DoubleArray &A, const double &v)
// Map vertices to global coordinates // Map vertices to global coordinates
for (int idx=0;idx<VertexCount;idx++) { for (int idx=0;idx<VertexCount;idx++) {
P = cellvertices[idx]; P = cellvertices(idx);
P.x += i; P.x += i;
P.y += j; P.y += j;
P.z += k; P.z += k;
cellvertices[idx] = P; cellvertices(idx) = P;
} }
} }
} }

59
analysis/dcel.h
View File

@ -8,53 +8,36 @@ Doubly-connected edge list (DECL)
// Vertex structure // Vertex structure
class Vertex{ class Vertex{
public: public:
Vertex() { d_data.resize(12); } Vertex();
~Vertex() = default; ~Vertex();
Vertex( const Vertex& ) = delete; void add(Point P);
Vertex operator=( const Vertex& ) = delete; void assign( int idx, Point P);
int size();
// Add/assign a point Point coords(int idx);
inline void add( const Point& P ) { d_data.push_back( P ); }
inline void assign( int idx, const Point& P ) { d_data[idx] = P; }
// Get a point
inline Point& coords( int idx ) { return d_data[idx]; }
inline const Point& coords( int idx ) const { return d_data[idx]; }
int IncidentEdge(); int IncidentEdge();
// Return the number of points
inline int size() const { return d_data.size(); }
private: private:
std::vector<Point> d_data; std::vector<double> vertex_data;
int size_;
}; };
// Halfedge structure // Halfedge structure
// Face // Face
class Halfedge{ class Halfedge{
public: public:
Halfedge() = default; Halfedge();
~Halfedge() = default; ~Halfedge();
Halfedge( const Halfedge& ) = delete;
Halfedge operator=( const Halfedge& ) = delete;
inline int v1(int edge) const { return d_data[edge][0]; } int v1(int edge);
inline int v2(int edge) const { return d_data[edge][1]; } int v2(int edge);
inline int face(int edge) const { return d_data[edge][2]; } int twin(int edge);
inline int twin(int edge) const { return d_data[edge][3]; } int face(int edge);
inline int prev(int edge) const { return d_data[edge][4]; } int next(int edge);
inline int next(int edge) const { return d_data[edge][5]; } int prev(int edge);
int size();
inline int size() const { return d_data.size(); }
inline void resize( int N ) { d_data.resize( N ); }
inline int& data( int i, int j ) { return d_data[j][i]; }
inline const int& data( int i, int j ) const { return d_data[j][i]; }
Array<int> data;
private: private:
std::vector<std::array<int,6>> d_data; int size_;
}; };
// DECL // DECL
@ -66,7 +49,7 @@ public:
int face(); int face();
Vertex vertex; Vertex vertex;
Halfedge halfedge; Halfedge halfedge;
void LocalIsosurface(const DoubleArray& A, double value, int i, int j, int k); void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k);
int Face(int index); int Face(int index);
double origin(int edge); double origin(int edge);
@ -76,5 +59,5 @@ public:
int VertexCount; int VertexCount;
private: private:
std::vector<int> FaceData; Array <int> FaceData;
}; };

29
example/Piston/Piston.py Normal file
View File

@ -0,0 +1,29 @@
import numpy
nx=96
ny=24
nz=24
N=nx*ny*nz
mesh=(nx,ny,nz)
data=numpy.ones(mesh,dtype=numpy.int8)
#print(data)
print("Writing piston")
print("Mesh size: "+repr(mesh))
radius = 8
# assign a bubble in the middle
for x in range(0,nx):
for y in range(0,ny):
for z in range(0,nz):
Y = y - ny/2
Z = z - nz/2
if Y*Y+Z*Z > radius*radius:
data[x,y,z]=0
elif x < 12:
data[x,y,z]=1
else:
data[x,y,z]=2
data.tofile("Piston.raw")

38
example/Piston/Piston.sh
View File

@ -1,37 +1,7 @@
#!/bin/bash #!/bin/bash
# Lines assigning various pressure BC for Color.in LBPM_DIR=../../tests
echo "0 1 1.01 0.99" > Color.in.pressures
echo "0 1 1.0125 0.9875" >> Color.in.pressures
echo "0 1 1.015 0.985" >> Color.in.pressures
echo "0 1 1.02 0.98" >> Color.in.pressures
echo "0 1 1.025 0.975" >> Color.in.pressures
echo "0 1 1.03 0.97" >> Color.in.pressures
for i in `seq 1 6`; do python Piston.py
# Set up cases for each boundary pressure pair mpirun -np 1 $LBPM_DIR/lbpm_serial_decomp input.db
dir="Case"$i mpirun -np 4 $LBPM_DIR/lbpm_color_simulator input.db
echo $dir
mkdir -p $dir
# copy the domain file
cp Domain.in $dir
# set up each case -- parameters are fixed in Color.in, with multiple cases to set the boundary pressure
sed -n '1p' Color.in > $dir/Color.in
sed -n '2p' Color.in >> $dir/Color.in
sed -n '3p' Color.in >> $dir/Color.in
sed -n '4p' Color.in >> $dir/Color.in
# sed -n '5p' Color.in >> $dir/Color.in
# print the pressure values into the input file
sed -n "${i}p" Color.in.pressures >> $dir/Color.in
sed -n '6p' Color.in >> $dir/Color.in
done
# simulations should be run using the following syntax
# PRE-PROCESSOR - set the radius to 18 voxel lengths
#mpirun -np 10 ~/install-LBPM-WIA/bin/lbpm_captube_pp 18 1
# RUN THE SIMULAUTION
#mpirun -np 10 ~/install-LBPM-WIA/bin/lbpm_color_simulator
exit;

28
example/Piston/input.db
View File

@ -1,34 +1,42 @@
Color { Color {
tau = 1.0; tauA = 0.7;
alpha = 1e-2; tauB = 0.7;
rhoA = 1.0;
rhoB = 1.0;
alpha = 1e-3;
beta = 0.95; beta = 0.95;
phi_s = 0.8;
wp_saturation = 0.7
F = 0, 0, 0 F = 0, 0, 0
Restart = false Restart = false
pBC = 0 pBC = 0
din = 1.0 din = 1.0
dout = 1.0 dout = 1.0
timestepMax = 200 timestepMax = 3000
interval = 1000 interval = 1000
tol = 1e-5; tol = 1e-5;
das = 0.1 das = 0.1
dbs = 0.9 dbs = 0.9
flux = 2.0
ComponentLabels = 0
ComponentAffinity = -1.0;
} }
Domain { Domain {
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz) Filename = "Piston.raw"
n = 16, 16, 16 // Size of local domain (Nx,Ny,Nz) nproc = 1, 1, 4 // Number of processors (Npx,Npy,Npz)
n = 24, 24, 24 // Size of local domain (Nx,Ny,Nz)
N = 24, 24, 96 // size of the input image
n_spheres = 1 // Number of spheres n_spheres = 1 // Number of spheres
L = 1, 1, 1 // Length of domain (x,y,z) L = 1, 1, 1 // Length of domain (x,y,z)
BC = 0 // Boundary condition type BC = 4 // Boundary condition type
ReadValues = 0, 1, 2
WriteValues = 0, 1, 2
} }
Analysis { Analysis {
blobid_interval = 1000 // Frequency to perform blob identification blobid_interval = 1000 // Frequency to perform blob identification
analysis_interval = 1000 // Frequency to perform analysis analysis_interval = 1000 // Frequency to perform analysis
restart_interval = 20000 // Frequency to write restart data restart_interval = 1000 // Frequency to write restart data
vis_interval = 20000 // Frequency to write visualization data visualization_interval = 1000 // Frequency to write visualization data
restart_file = "Restart" // Filename to use for restart file (will append rank) restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent" load_balance = "independent" // Load balance method to use: "none", "default", "independent"