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

merge dox

Browse Source
This commit is contained in:
JamesEMcclure 2021-09-20 09:22:38 -04:00
commit 3cf0bc67d8
14 changed files with 575 additions and 774 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
analysis/ElectroChemistry.cpp
View File

@ -1,5 +1,6 @@
#include "analysis/ElectroChemistry.h"
ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr <Domain> dm):
Dm(dm)
{

61
analysis/FlowAdaptor.h
View File

@ -12,15 +12,74 @@
#include "models/ColorModel.h"
/**
* \class FlowAdaptor
*
* @brief
* The FlowAdaptor class operates on a lattice Boltzmann model to alter the flow conditions
*
*/
class FlowAdaptor{
public:
/**
* \brief Create a flow adaptor to operate on the LB model
* @param M ScaLBL_ColorModel
*/
FlowAdaptor(ScaLBL_ColorModel &M);
/**
* \brief Destructor
*/
~FlowAdaptor();
/**
* \brief Fast-forward interface motion
* \details Accelerate the movement of interfaces based on the time derivative
* Optional keys to control behavior can be specified in the input database:
* move_interface_cutoff -- identifies the diffuse interface region
* move_interface_factor -- determines how much to ``fast forward"
* @param M ScaLBL_ColorModel
*/
double MoveInterface(ScaLBL_ColorModel &M);
/**
* \brief Image re-initialization
* \details Re-initialize LB simulation from image data
* @param M ScaLBL_ColorModel
* @param Filename name of input file to be used to read image
*/
double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
/**
* \details Update volume fraction based on morphological algorithm. Dilation / erosion algorithm will be applied to
* grow / shrink the phase regions
* @param M ScaLBL_ColorModel
* @param delta_volume target change in volume fraction
*/
double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
/**
* \details Update fractional flow condition. Mass will be preferentially added or removed from
* phase regions based on where flow is occurring
* @param M ScaLBL_ColorModel
*/ double UpdateFractionalFlow(ScaLBL_ColorModel &M);
/**
* \brief image re-initialization
* \details Re-initialize LB simulation from image data
* @param M ScaLBL_ColorModel
* @param seed_water_in_oil controls amount of mass to randomly seed into fluids
*/
double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
/**
* \brief Re-initialize LB simulation
* @param M ScaLBL_ColorModel
*/
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;
DoubleArray phi_t;

8
analysis/FreeEnergy.h
View File

@ -19,6 +19,14 @@
#include "IO/Writer.h"
#include "models/FreeLeeModel.h"
/**
* \class FreeEnergyAnalyzer
*
* @brief
* The FreeEnergyAnalyzer class is constructed to analyze the LBPM free energy model for liquid-gas systems
*
*/
class FreeEnergyAnalyzer{
public:
std::shared_ptr <Domain> Dm;

16
analysis/GreyPhase.h
View File

@ -15,6 +15,14 @@
#include "IO/Reader.h"
#include "IO/Writer.h"
/**
* \class GreyPhase
*
* @brief
* The GreyPhase class tracks pressure, mass and momentum within a grey phase
*
*/
class GreyPhase{
public:
double p;
@ -26,6 +34,14 @@ class GreyPhase{
private:
};
/**
* \class GreyPhaseAnalysis
*
* @brief
* The GreyPhaseAnalysis class is constructed to analyze the LBPM greyscale model
*
*/
class GreyPhaseAnalysis{
public:
std::shared_ptr <Domain> Dm;

2
analysis/Minkowski.cpp
View File

@ -66,7 +66,7 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue)
{
PROFILE_START("ComputeScalar");
Xi = Ji = Ai = 0.0;
DECL object;
DCEL object;
int e1,e2,e3;
double s,s1,s2,s3;
double a1,a2,a3;

50
analysis/Minkowski.h
View File

@ -34,6 +34,14 @@
#include "IO/Reader.h"
#include "IO/Writer.h"
/**
* \class Minkowski
*
* @brief
* The Minkowski class is constructed to analyze the geometric properties of structures based on the Minkowski functionals
*
*/
class Minkowski{
//...........................................................................
@ -61,6 +69,7 @@ public:
int n_connected_components;
//...........................................................................
int Nx,Ny,Nz;
double V(){
return Vi;
}
@ -75,15 +84,56 @@ public:
}
//..........................................................................
/**
* \brief Null constructor
*/
Minkowski(){};//NULL CONSTRUCTOR
/**
* \brief Constructor based on an existing Domain
* @param Dm - Domain structure
*/
Minkowski(std::shared_ptr <Domain> Dm);
~Minkowski();
/**
* \brief Compute scalar minkowski functionals
* step 1. compute the distance to an object
* step 2. construct dcel to represent the isosurface
* step 3. compute the scalar Minkowski functionals
* THIS ALGORITHM ASSUMES THAT id() is populated with phase id to distinguish objects
* 0 - labels the object
* 1 - labels everything else
*/
void MeasureObject();
void MeasureObject(double factor, const DoubleArray &Phi);
/**
* \details Compute scalar minkowski functionals for connected part of a structure
* step 1. compute connected components and extract largest region by volume
* step 2. compute the distance to the connected part of the structure
* step 3. construct dcel to represent the isosurface
* step 4. compute the scalar Minkowski functionals
* THIS ALGORITHM ASSUMES THAT id() is populated with phase id to distinguish objects
* 0 - labels the object
* 1 - labels everything else
*/
int MeasureConnectedPathway();
int MeasureConnectedPathway(double factor, const DoubleArray &Phi);
/**
* \brief Compute scalar minkowski functionals
* \details Construct an isosurface and return the geometric invariants based on the triangulated list
* @param isovalue - threshold value to use to determine iso-surface
* @param Field - DoubleArray containing the field to threshold
*/
void ComputeScalar(const DoubleArray& Field, const double isovalue);
/**
* \brief print the scalar invariants
*/
void PrintAll();
};

18
analysis/dcel.cpp
View File

@ -1,19 +1,19 @@
#include "analysis/dcel.h"
DECL::DECL(){
DCEL::DCEL(){
}
DECL::~DECL(){
DCEL::~DCEL(){
TriangleCount=0;
VertexCount=0;
}
int DECL::Face(int index){
int DCEL::Face(int index){
return FaceData[index];
}
void DECL::Write(){
void DCEL::Write(){
int e1,e2,e3;
FILE *TRIANGLES;
TRIANGLES = fopen("triangles.stl","w");
@ -32,7 +32,7 @@ void DECL::Write(){
fclose(TRIANGLES);
}
void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, const int j, const int k){
void DCEL::LocalIsosurface(const DoubleArray& A, double value, const int i, const int j, const int k){
Point P,Q;
Point PlaceHolder;
Point C0,C1,C2,C3,C4,C5,C6,C7;
@ -174,7 +174,7 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
}
int nTris = TriangleCount;
// Now add the local values to the DECL data structure
// Now add the local values to the DCEL data structure
if (nTris>0){
FaceData.resize(TriangleCount);
//printf("Construct halfedge structure... \n");
@ -250,7 +250,7 @@ void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, cons
}
}
Point DECL::TriNormal(int edge)
Point DCEL::TriNormal(int edge)
{
Point P,Q,R;
Point U,V,W;
@ -294,7 +294,7 @@ Point DECL::TriNormal(int edge)
return W;
}
double DECL::EdgeAngle(int edge)
double DCEL::EdgeAngle(int edge)
{
double angle;
double dotprod;
@ -369,7 +369,7 @@ double DECL::EdgeAngle(int edge)
void iso_surface(const Array<double>&Field, const double isovalue)
{
DECL object;
DCEL object;
int e1,e2,e3;
FILE *TRIANGLES;
TRIANGLES = fopen("isosurface.stl","w");

22
analysis/dcel.h
View File

@ -4,8 +4,9 @@
#include <vector>
#include "analysis/pmmc.h"
/*
Doubly-connected edge list (DECL)
/**
* \class Vertex
* @brief store vertex for DCEL data structure
*/
// Vertex structure
@ -34,8 +35,10 @@ private:
};
// Halfedge structure
// Face
/**
* \class Halfedge
* @brief store half edge for DCEL data structure
*/
class Halfedge{
public:
Halfedge() = default;
@ -60,11 +63,14 @@ private:
std::vector<std::array<int,6>> d_data;
};
// DECL
class DECL{
/**
* \class DCEL
* @details doubly connected edge list data structure
*/
class DCEL{
public:
DECL();
~DECL();
DCEL();
~DCEL();
int face();
Vertex vertex;

355
common/ScaLBL.h
View File

@ -14,8 +14,8 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/* ScaLBL.h
* Header file for Scalable Lattice Boltzmann Library
/** @file ScaLBL.h */
/* \detail Header file for Scalable Lattice Boltzmann Library
* Separate implementations for GPU and CPU must both follow the conventions defined in this header
* This libarry contains the essential components of the LBM
* - streaming implementations
@ -27,48 +27,194 @@
#define ScalLBL_H
#include "common/Domain.h"
/**
* \brief Set compute device
* @param rank rank of MPI process
*/
extern "C" int ScaLBL_SetDevice(int rank);
/**
* \brief Allocate memory
* @param address memory address
* @param size size in bytes
*/
extern "C" void ScaLBL_AllocateDeviceMemory(void** address, size_t size);
/**
* \brief Free memory
* @param pointer pointer to memory to free
*/
extern "C" void ScaLBL_FreeDeviceMemory(void* pointer);
/**
* \brief Copy memory from host to device
* \detail Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
* @param source memory region to copy from
* @param size size of the region to copy in bytes
*/
extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size);
/**
* \brief Copy memory from device to host
* \detail Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
* @param source memory region to copy from
* @param size size of the region to copy in bytes
*/
extern "C" void ScaLBL_CopyToHost(void* dest, const void* source, size_t size);
/**
=* \brief Allocate zero copy memory buffer (i.e. shared memory)
* @param address memory address
* @param size size in bytes
*/
extern "C" void ScaLBL_AllocateZeroCopy(void** address, size_t size);
/**
* \brief Copy memory from host to zero copy buffer
* \detail Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
* @param source memory region to copy from
* @param size size of the region to copy in bytes
*/
extern "C" void ScaLBL_CopyToZeroCopy(void* dest, const void* source, size_t size);
/**
* \brief Device barrier routine
*/
extern "C" void ScaLBL_DeviceBarrier();
/**
* \brief Pack D3Q19 distributions for communication
* @param q - index for distribution based on D3Q19 discrete velocity structure
* @param list - list of distributions to communicate
* @param start - index to start parsing the list
* @param count - number of values to pack
* @param sendbuf - memory buffer to hold values that will be sent
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Pack(int q, int *list, int start, int count, double *sendbuf, double *dist, int N);
/**
* \brief Unpack D3Q19 distributions after communication
* @param q - index for distribution based on D3Q19 discrete velocity structure
* @param list - list of distributions to communicate
* @param start - index to start parsing the list
* @param count - number of values to unppack
* @param revbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N);
/**
* \brief Unpack D3Q7 distributions after communication
* @param q - index for distribution based on D3Q19 discrete velocity structure
* @param list - list of distributions to communicate
* @param start - index to start parsing the list
* @param count - number of values to unppack
* @param revbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_Unpack(int q, int *list, int start, int count, double *recvbuf, double *dist, int N);
/**
* \brief Pack halo for scalar field to be prepare for communication
* @param list - list of distributions to communicate
* @param count - number of values to ppack
* @param sendbuf - memory buffer to pack values into
* @param Data - scalar field
* @param N - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_Scalar_Pack(int *list, int count, double *sendbuf, double *Data, int N);
/**
* \brief Pack halo for scalar field to be prepare for communication
* @param list - list of distributions to communicate
* @param count - number of values to unppack
* @param recvbuf - memory buffer where recieved values have been stored
* @param Data - scalar field
* @param N - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_Scalar_Unpack(int *list, int count, double *recvbuf, double *Data, int N);
/**
* \brief Unpack values and compute Shan-Chen type of gradient
* @param weight - weight value for gradient sum
* @param Cqx - contribution to x-part of gradient
* @param Cqy - contribution to y-part of gradient
* @param Cqz - contribution to z-part of gradient
* @param list - list of distributions to communicate
* @param start - location to start reading the list
* @param count - number of values to unppack
* @param recvbuf - memory buffer where recieved values have been stored
* @param phi - scalar field
* @param grad - gradient
* @param N - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_Gradient_Unpack(double weight, double Cqx, double Cqy, double Cqz,
int *list, int start, int count, double *recvbuf, double *phi, double *grad, int N);
extern "C" void ScaLBL_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N);
extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N);
/**
* \brief Initialize D3Q19 distributions
* @param Dist - D3Q19 distributions
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
/**
* \brief Compute momentum from D3Q19 distribution
* @param Dist - D3Q19 distributions
* @param vel - memory buffer to store the momentum that is computed
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
/**
* \brief Compute pressure from D3Q19 distribution
* @param Dist - D3Q19 distributions
* @param press - memory buffer to store the pressure field that is computed
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *press, int Np);
// BGK MODEL
/**
* \brief BGK collision based on AA even access pattern for D3Q19
* @param dist - D3Q19 distributions
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
* @param rlx - relaxation parameter
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
*/
extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz);
/**
* \brief BGK collision based on AA odd access pattern for D3Q19
* @param neighborList - neighbors based on D3Q19 lattice structure
* @param dist - D3Q19 distributions
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
* @param rlx - relaxation parameter
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
*/
extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz);
// GREYSCALE MODEL (Single-component)
@ -136,26 +282,73 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
// LBM Poisson solver
/**
* \brief Poisson-Boltzmann collision based on AA odd access pattern for D3Q19
* @param neighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Den_charge - charge density
* @param Psi -
* @param ElectricField - electric field
* @param tau - relaxation time
* @param epsilon_LB -
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
/**
* \brief Poisson-Boltzmann collision based on AA even access pattern for D3Q7
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Den_charge - charge density
* @param Psi -
* @param ElectricField - electric field
* @param tau - relaxation time
* @param epsilon_LB -
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau,
double epsilon_LB, int start, int finish, int Np);
/**
* \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7
* @param neighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Psi -
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np);
/**
* \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Psi -
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np);
/**
* \brief Initialize Poisson-Boltzmann solver
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Psi -
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish);
//maybe deprecated
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np);
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
@ -165,27 +358,138 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);
// MRT MODEL
/**
* \brief MRT collision based on AA even access pattern for D3Q19
* @param dist - D3Q19 distributions
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
* @param rlx_setA - relaxation parameter for viscous modes
* @param rlx_setB - relaxation parameter for non-viscous modes
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
*/
extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
double Fy, double Fz);
extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *d_neighborList, double *dist, int start, int finish, int Np,
/**
* \brief MRT collision based on AA even access pattern for D3Q19
* @param neighborList - index of neighbors based on D3Q19 lattice structure
* @param dist - D3Q19 distributions
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
* @param rlx_setA - relaxation parameter for viscous modes
* @param rlx_setB - relaxation parameter for non-viscous modes
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
*/
extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *neighborList, double *dist, int start, int finish, int Np,
double rlx_setA, double rlx_setB, double Fx, double Fy, double Fz);
// COLOR MODEL
/**
* \brief Color model collision based on AA even access pattern for D3Q19
* @param Map - mapping between sparse and dense data structures
* @param dist - D3Q19 distributions
* @param Aq - D3Q7 distribution for component A
* @param Bq - D3Q7 distribution for component B
* @param Den - density field
* @param Phi - phase indicator field
* @param Vel - velocity field
* @param rhoA - density of component A
* @param rhoB - density of component B
* @param tauA - relaxation time for component A
* @param tauB - relaxation time for component B
* @param alpha - parameter to control interfacial tension
* @param beta - parameter to control interface width
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
* @param strideY - stride in y-direction for gradient computation
* @param strideZ - stride in z-direction for gradient computation
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_Color(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
/**
* \brief Color model collision based on AA even access pattern for D3Q19
* @param NeighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense data structures
* @param dist - D3Q19 distributions
* @param Aq - D3Q7 distribution for component A
* @param Bq - D3Q7 distribution for component B
* @param Den - density field
* @param Phi - phase indicator field
* @param Vel - velocity field
* @param rhoA - density of component A
* @param rhoB - density of component B
* @param tauA - relaxation time for component A
* @param tauB - relaxation time for component B
* @param alpha - parameter to control interfacial tension
* @param beta - parameter to control interface width
* @param Fx - force in x direction
* @param Fy - force in y direction
* @param Fz - force in z direction
* @param strideY - stride in y-direction for gradient computation
* @param strideZ - stride in z-direction for gradient computation
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_AAodd_Color(int *NeighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
/**
* \brief Compute phase field based on AA odd access pattern for D3Q19
* @param NeighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense data structures
* @param Aq - D3Q7 distribution for component A
* @param Bq - D3Q7 distribution for component B
* @param Den - density field
* @param Phi - phase indicator field
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *NeighborList, int *Map, double *Aq, double *Bq,
double *Den, double *Phi, int start, int finish, int Np);
/**
* \brief Compute phase field based on AA even access pattern for D3Q19
* @param Map - mapping between sparse and dense data structures
* @param Aq - D3Q7 distribution for component A
* @param Bq - D3Q7 distribution for component B
* @param Den - density field
* @param Phi - phase indicator field
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi,
int start, int finish, int Np);
/**
* \brief Initialize phase field for color model
* @param Map - mapping between sparse and dense data structures
* @param Phi - phase indicator field
* @param Den - density field
* @param Aq - D3Q7 distribution for component A
* @param Bq - D3Q7 distribution for component B
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, int *Map, double *Aq, double *Bq, double *Den,
double *Phi, double *ColorGrad, double *Vel, double rhoA, double rhoB, double beta, int start, int finish, int Np);
@ -196,7 +500,6 @@ extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad,
extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz);
extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np);
// Density functional hydrodynamics LBM
extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np);
@ -347,12 +650,24 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion_Flux_DiffAdvcElec_BC_Z(int *list, double
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Flux_DiffAdvcElec_BC_z(int *d_neighborList, int *list, double *dist, double Cin, double tau, double *VelocityZ,double *ElectricField,double Di,double zi,double Vt, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Flux_DiffAdvcElec_BC_Z(int *d_neighborList, int *list, double *dist, double Cout, double tau, double *VelocityZ,double *ElectricField,double Di,double zi,double Vt, int count, int Np);
/**
* \class ScaLBL_Communicator
*
* @brief Generalized communication routines for lattice Boltzmann methods.
*
*/
class ScaLBL_Communicator{
public:
//......................................................................................
ScaLBL_Communicator(std::shared_ptr <Domain> Dm);
/**
*\brief Constructor
* @param Dm - Domain information
*/
ScaLBL_Communicator(std::shared_ptr <Domain> Dm);
//ScaLBL_Communicator(Domain &Dm, IntArray &Map);
/**
*\brief Destructor
*/
~ScaLBL_Communicator();
//......................................................................................
unsigned long int CommunicationCount,SendCount,RecvCount;

2
doxygen/Doxyfile.in
View File

@ -1173,7 +1173,7 @@ HIDE_UNDOC_RELATIONS = YES
# toolkit from AT&T and Lucent Bell Labs. The other options in this section
# have no effect if this option is set to NO (the default)
HAVE_DOT = YES
HAVE_DOT = NO
# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen
# will generate a graph for each documented class showing the direct and

12
doxygen/DoxygenMainpage.h
View File

@ -1,9 +1,17 @@
/** \mainpage LBPM
*
* This is the documentation for LBPM
* C/C++ routines
*
* - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)"
* - \ref models "Lattice Boltzmann models"
* - \ref ColorModel "Color model"
* - \ref analysis "Analysis routines"
* - \ref FlowAdaptor "FlowAdaptor"
* - \ref DCEL "Doubly connected edge list"
* - \ref Minkowski "Minkowski functionals"
* - \ref IO "IO routines"
* - \ref Utilities "Utility routines"
* - \ref tests "Unit tests"
*
* \author James McClure
* \author J.E. McClure, M. Berrill
*/

734
models/ColorModel.cpp
View File

@ -959,7 +959,11 @@ double ScaLBL_ColorModel::Run(int returntime){
void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
int analysis_interval = 1000; // number of timesteps in between in situ analysis
if (analysis_db->keyExists( "analysis_interval" )){
analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
}
/*
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
@ -978,7 +982,6 @@ void ScaLBL_ColorModel::Run(){
int CURRENT_MORPH_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
int CURRENT_STEADY_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
int morph_interval = 100000;
int analysis_interval = 1000; // number of timesteps in between in situ analysis
int morph_timesteps = 0;
double morph_delta = 0.0;
double seed_water = 0.0;
@ -989,7 +992,6 @@ void ScaLBL_ColorModel::Run(){
double delta_volume = 0.0;
double delta_volume_target = 0.0;
/* history for morphological algoirthm */
double KRA_MORPH_FACTOR=0.5;
double volA_prev = 0.0;
double log_krA_prev = 1.0;
@ -1041,9 +1043,7 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "tolerance" )){
tolerance = analysis_db->getScalar<double>( "tolerance" );
}
if (analysis_db->keyExists( "analysis_interval" )){
analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
}
if (analysis_db->keyExists( "min_steady_timesteps" )){
MIN_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "min_steady_timesteps" );
}
@ -1054,7 +1054,6 @@ void ScaLBL_ColorModel::Run(){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
/* defaults for simulation protocols */
auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){
// Get the list of images
@ -1100,8 +1099,7 @@ void ScaLBL_ColorModel::Run(){
double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB;
}
}
if (rank==0){
} if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
printf(" using protocol = image sequence \n");
@ -1156,7 +1154,7 @@ void ScaLBL_ColorModel::Run(){
printf("No. of timesteps: %i \n", timestepMax);
fflush(stdout);
}
*/
//************ MAIN ITERATION LOOP ***************************************/
comm.barrier();
PROFILE_START("Loop");
@ -1256,242 +1254,6 @@ void ScaLBL_ColorModel::Run(){
}
// Run the analysis
analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
// allow initial ramp-up to get closer to steady state
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
analysis.finish();
CURRENT_STEADY_TIMESTEPS += analysis_interval;
double volB = Averages->gwb.V;
double volA = Averages->gnb.V;
volA /= Dm->Volume;
volB /= Dm->Volume;
//initial_volume = volA*Dm->Volume;
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
double vB_x = Averages->gwb.Px/Averages->gwb.M;
double vB_y = Averages->gwb.Py/Averages->gwb.M;
double vB_z = Averages->gwb.Pz/Averages->gwb.M;
double muA = rhoA*(tauA-0.5)/3.f;
double muB = rhoB*(tauB-0.5)/3.f;
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
double current_saturation = volB/(volA+volB);
double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
if ( morph_timesteps > morph_interval ){
bool isSteady = false;
if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_STEADY_TIMESTEPS > MIN_STEADY_TIMESTEPS))
isSteady = true;
if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
isSteady = true;
if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_STEADY_TIMESTEPS > RESCALE_FORCE_AFTER_TIMESTEP){
RESCALE_FORCE = false;
double RESCALE_FORCE_FACTOR = capillary_number / Ca;
if (RESCALE_FORCE_FACTOR > 2.0) RESCALE_FORCE_FACTOR = 2.0;
if (RESCALE_FORCE_FACTOR < 0.5) RESCALE_FORCE_FACTOR = 0.5;
Fx *= RESCALE_FORCE_FACTOR;
Fy *= RESCALE_FORCE_FACTOR;
Fz *= RESCALE_FORCE_FACTOR;
force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
color_db->putVector<double>("F",{Fx,Fy,Fz});
}
if ( isSteady ){
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
//****** ENDPOINT ADAPTATION ********/
double krA_TMP= fabs(muA*flow_rate_A / force_mag);
double krB_TMP= fabs(muB*flow_rate_B / force_mag);
log_krA = log(krA_TMP);
if (krA_TMP < 0.0){
// cannot do endpoint adaptation if kr is negative
log_krA = log_krA_prev;
}
else if (krA_TMP < krB_TMP && morph_delta > 0.0){
/** morphological target based on relative permeability for A **/
log_krA_target = log(KRA_MORPH_FACTOR*(krA_TMP));
slope_krA_volume = (log_krA - log_krA_prev)/(Dm->Volume*(volA - volA_prev));
delta_volume_target=min(delta_volume_target,Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume));
if (rank==0){
printf(" Enabling endpoint adaptation: krA = %f, krB = %f \n",krA_TMP,krB_TMP);
printf(" log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
}
}
log_krA_prev = log_krA;
volA_prev = volA;
//******************************** **/
/** compute averages & write data **/
Averages->Full();
Averages->Write(timestep);
analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.finish();
if (rank==0){
printf("** WRITE STEADY POINT *** ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
double h = Dm->voxel_length;
// pressures
double pA = Averages->gnb.p;
double pB = Averages->gwb.p;
double pAc = Averages->gnc.p;
double pBc = Averages->gwc.p;
double pAB = (pA-pB)/(h*6.0*alpha);
double pAB_connected = (pAc-pBc)/(h*6.0*alpha);
// connected contribution
double Vol_nc = Averages->gnc.V/Dm->Volume;
double Vol_wc = Averages->gwc.V/Dm->Volume;
double Vol_nd = Averages->gnd.V/Dm->Volume;
double Vol_wd = Averages->gwd.V/Dm->Volume;
double Mass_n = Averages->gnc.M + Averages->gnd.M;
double Mass_w = Averages->gwc.M + Averages->gwd.M;
double vAc_x = Averages->gnc.Px/Mass_n;
double vAc_y = Averages->gnc.Py/Mass_n;
double vAc_z = Averages->gnc.Pz/Mass_n;
double vBc_x = Averages->gwc.Px/Mass_w;
double vBc_y = Averages->gwc.Py/Mass_w;
double vBc_z = Averages->gwc.Pz/Mass_w;
// disconnected contribution
double vAd_x = Averages->gnd.Px/Mass_n;
double vAd_y = Averages->gnd.Py/Mass_n;
double vAd_z = Averages->gnd.Pz/Mass_n;
double vBd_x = Averages->gwd.Px/Mass_w;
double vBd_y = Averages->gwd.Py/Mass_w;
double vBd_z = Averages->gwd.Pz/Mass_w;
double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z);
double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z);
double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z);
double flow_rate_B_disconnected = (Vol_wd)*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z);
double kAeff_connected = h*h*muA*flow_rate_A_connected/(force_mag);
double kBeff_connected = h*h*muB*flow_rate_B_connected/(force_mag);
double kAeff_disconnected = h*h*muA*flow_rate_A_disconnected/(force_mag);
double kBeff_disconnected = h*h*muB*flow_rate_B_disconnected/(force_mag);
double kAeff = h*h*muA*(flow_rate_A)/(force_mag);
double kBeff = h*h*muB*(flow_rate_B)/(force_mag);
double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag;
double Mobility = muA/muB;
bool WriteHeader=false;
FILE * kr_log_file = fopen("relperm.csv","r");
if (kr_log_file != NULL)
fclose(kr_log_file);
else
WriteHeader=true;
kr_log_file = fopen("relperm.csv","a");
if (WriteHeader)
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
fclose(kr_log_file);
printf(" Measured capillary number %f \n ",Ca);
}
if (SET_CAPILLARY_NUMBER ){
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
color_db->putVector<double>("F",{Fx,Fy,Fz});
}
CURRENT_STEADY_TIMESTEPS = 0;
}
else{
if (rank==0){
printf("** Continue to simulate steady *** \n ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
}
}
morph_timesteps=0;
Ca_previous = Ca;
}
if (MORPH_ADAPT ){
CURRENT_MORPH_TIMESTEPS += analysis_interval;
if (USE_DIRECT){
// Use image sequence
IMAGE_INDEX++;
MORPH_ADAPT = false;
if (IMAGE_INDEX < IMAGE_COUNT){
std::string next_image = ImageList[IMAGE_INDEX];
if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
color_db->putScalar<int>("image_index",IMAGE_INDEX);
ImageInit(next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
timestep = timestepMax;
}
}
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
double massChange = SeedPhaseField(seed_water);
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target);
}
else if (USE_MORPHOPEN_OIL){
delta_volume = volA*Dm->Volume - initial_volume;
if (rank==0) printf("***Morphological opening of connected oil, target volume change %f ***\n", delta_volume_target);
MorphOpenConnected(delta_volume_target);
}
else {
if (rank==0) printf("***Shell aggregation, target volume change %f ***\n", delta_volume_target);
//double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A
delta_volume += MorphInit(beta,delta_volume_target-delta_volume);
}
if ( (delta_volume - delta_volume_target)/delta_volume_target > 0.0 ){
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
else if (!(USE_DIRECT) && CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
RESCALE_FORCE = true;
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
}
morph_timesteps += analysis_interval;
}
comm.barrier();
}
analysis.finish();
PROFILE_STOP("Loop");
@ -1515,486 +1277,6 @@ void ScaLBL_ColorModel::Run(){
// ************************************************************************
}
double ScaLBL_ColorModel::ImageInit(std::string Filename){
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; // save what was read
double *PhaseLabel;
PhaseLabel = new double[Nx*Ny*Nz];
AssignComponentLabels(PhaseLabel);
double Count = 0.0;
double PoreCount = 0.0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (id[Nx*Ny*k+Nx*j+i] == 2){
PoreCount++;
Count++;
}
else if (id[Nx*Ny*k+Nx*j+i] == 1){
PoreCount++;
}
}
}
}
Count=Dm->Comm.sumReduce( Count);
PoreCount=Dm->Comm.sumReduce( PoreCount);
if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
comm.barrier();
ScaLBL_D3Q19_Init(fq, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
comm.barrier();
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
double saturation = Count/PoreCount;
return saturation;
}
double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
int nx = Nx;
int ny = Ny;
int nz = Nz;
int n;
int N = nx*ny*nz;
double volume_change=0.0;
if (target_volume_change < 0.0){
Array<char> id_solid(nx,ny,nz);
Array<int> phase_label(nx,ny,nz);
DoubleArray distance(Nx,Ny,Nz);
DoubleArray phase(nx,ny,nz);
signed char *id_connected;
id_connected = new signed char [nx*ny*nz];
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
// Extract only the connected part of NWP
double vF=0.0; double vS=0.0;
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
comm.barrier();
long long count_connected=0;
long long count_porespace=0;
long long count_water=0;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( phase_label(i,j,k) == 0){
count_connected++;
}
if (id[n] > 0){
count_porespace++;
}
if (id[n] == 2){
count_water++;
}
}
}
}
count_connected=Dm->Comm.sumReduce( count_connected);
count_porespace=Dm->Comm.sumReduce( count_porespace);
count_water=Dm->Comm.sumReduce( count_water);
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( phase_label(i,j,k) == 0){
id_solid(i,j,k) = 1;
id_connected[n] = 2;
id[n] = 2;
/* delete the connected component */
phase(i,j,k) = -1.0;
}
else{
id_solid(i,j,k) = 0;
id_connected[n] = 0;
}
}
}
}
CalcDist(distance,id_solid,*Dm);
signed char water=2;
signed char notwater=1;
double SW=-(target_volume_change)/count_connected;
MorphOpen(distance, id_connected, Dm, SW, water, notwater);
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( id_connected[n] == 1){
id_solid(i,j,k) = 0;
}
else{
id_solid(i,j,k) = 1;
}
}
}
}
CalcDist(distance,id_solid,*Dm);
// re-initialize
double beta = 0.95;
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
double d = distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){
if (d < 3.f){
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
int count_morphopen=0.0;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n=k*nx*ny+j*nx+i;
// only apply opening to connected component
if ( id_connected[n] == 1){
count_morphopen++;
}
}
}
}
count_morphopen=Dm->Comm.sumReduce( count_morphopen);
volume_change = double(count_morphopen - count_connected);
if (rank==0) printf(" opening of connected oil %f \n",volume_change/count_connected);
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
}
return(volume_change);
}
double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
double count =0.f;
double *Aq_tmp, *Bq_tmp;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
count= Dm->Comm.sumReduce( count);
mass_loss= Dm->Comm.sumReduce( mass_loss);
if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
ScaLBL_CopyToDevice(Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(Bq, Bq_tmp, 7*Np*sizeof(double));
return(mass_loss);
}
double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta_volume){
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
double vF = 0.f;
double vS = 0.f;
double delta_volume;
double WallFactor = 1.0;
bool USE_CONNECTED_NWP = false;
DoubleArray phase(Nx,Ny,Nz);
IntArray phase_label(Nx,Ny,Nz);;
DoubleArray phase_distance(Nx,Ny,Nz);
Array<char> phase_id(Nx,Ny,Nz);
fillHalo<double> fillDouble(Dm->Comm,Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
// Basic algorithm to
// 1. Copy phase field to CPU
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
double count = 0.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
double volume_initial = Dm->Comm.sumReduce( count);
double PoreVolume = Dm->Volume*Dm->Porosity();
/*ensure target isn't an absurdly small fraction of pore volume */
if (volume_initial < target_delta_volume*PoreVolume){
volume_initial = target_delta_volume*PoreVolume;
}
/*
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb");
fwrite(phase.data(),8,N,INPUT);
fclose(INPUT);
*/
// 2. Identify connected components of phase field -> phase_label
double volume_connected = 0.0;
double second_biggest = 0.0;
if (USE_CONNECTED_NWP){
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
comm.barrier();
// only operate on component "0"
count = 0.0;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int label = phase_label(i,j,k);
if (label == 0 ){
phase_id(i,j,k) = 0;
count += 1.0;
}
else
phase_id(i,j,k) = 1;
if (label == 1 ){
second_biggest += 1.0;
}
}
}
}
volume_connected = Dm->Comm.sumReduce( count);
second_biggest = Dm->Comm.sumReduce( second_biggest);
}
else {
// use the whole NWP
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (Averages->SDs(i,j,k) > 0.f){
if (phase(i,j,k) > 0.f ){
phase_id(i,j,k) = 0;
}
else {
phase_id(i,j,k) = 1;
}
}
else {
phase_id(i,j,k) = 1;
}
}
}
}
}
/*int reach_x, reach_y, reach_z;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
}
}
}*/
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*Dm);
double temp,value;
double factor=0.5/beta;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 3.f ){
value = phase(i,j,k);
if (value > 1.f) value=1.f;
if (value < -1.f) value=-1.f;
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
temp = -factor*log((1.0+value)/(1.0-value));
/// use this approximation close to the object
if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
// erase the original object
phase(i,j,k) = -1.0;
}
}
}
}
if (USE_CONNECTED_NWP){
if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
// if connected volume is less than 2% just delete the whole thing
if (rank==0) printf("Connected region has shrunk! \n");
REVERSE_FLOW_DIRECTION = true;
}
/* else{*/
if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
}
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
double target_delta_volume_incremental = target_delta_volume;
if (fabs(target_delta_volume) > 0.01*volume_initial)
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
else phase_id(i,j,k) = 1;
//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
}
}
}
CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
// 5. Update phase indicator field based on new distnace
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double d = phase_distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){
if (d < 3.f){
//phase(i,j,k) = -1.0;
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
fillDouble.fill(phase);
//}
count = 0.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f){
count+=1.f;
}
}
}
}
double volume_final= Dm->Comm.sumReduce( count);
delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
if (rank == 0) printf(" new saturation = %f \n", volume_final/(Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs)));
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
/*
sprintf(LocalRankFilename,"dist_final.%05i.raw",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(phase_distance.data(),8,N,DIST);
fclose(DIST);
sprintf(LocalRankFilename,"phi_final.%05i.raw",rank);
FILE *PHI = fopen(LocalRankFilename,"wb");
fwrite(phase.data(),8,N,PHI);
fclose(PHI);
*/
// 7. Re-initialize phase field and density
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
return delta_volume;
}
void ScaLBL_ColorModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);

66
models/ColorModel.h
View File

@ -37,21 +37,78 @@ Implementation of color lattice boltzmann model
#ifndef ScaLBL_ColorModel_INC
#define ScaLBL_ColorModel_INC
/**
* \class ScaLBL_ColorModel
*
* @details
* The ScaLBL_ColorModel class contains routines to initialize and run a two-component color lattice Boltzmann model
* Momentum transport equations are described by a D3Q19 scheme
* Mass transport equations are described by D3Q7 scheme
*/
class ScaLBL_ColorModel{
public:
/**
* \brief Constructor
* @param RANK processor rank
* @param NP number of processors
* @param COMM MPI communicator
*/
ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM);
~ScaLBL_ColorModel();
// functions in they should be run
/**
* \brief Read simulation parameters
* @param filename input database file that includes "Color" section
*/
void ReadParams(string filename);
/**
* \brief Read simulation parameters
* @param db0 input database that includes "Color" section
*/
void ReadParams(std::shared_ptr<Database> db0);
/**
* \brief Create domain data structures
*/
void SetDomain();
/**
* \brief Read image data
*/
void ReadInput();
/**
* \brief Create color model data structures
*/
void Create();
/**
* \brief Initialize the simulation
*/
void Initialize();
/**
* \brief Run the simulation
*/
void Run();
/**
* \brief Run the simulation
* @param returntime - timestep at which the routine will return
*/
double Run(int returntime);
/**
* \brief Debugging function to dump simulation state to disk
*/
void WriteDebug();
/**
* \brief Copy the phase field for use by external methods
* @param f - DoubleArray to hold the phase field
*/
void getPhaseField(DoubleArray &f);
bool Restart,pBC;
@ -90,6 +147,9 @@ public:
double *Velocity;
double *Pressure;
/**
* \brief Assign wetting affinity values
*/
void AssignComponentLabels(double *phase);
private:
@ -104,10 +164,6 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
double ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);
double MorphOpenConnected(double target_volume_change);
};
#endif

2
tests/test_dcel_tri_normal.cpp
View File

@ -54,7 +54,7 @@ int main(int argc, char **argv)
}
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
DECL object;
DCEL object;
Point P1,P2,P3;
Point U,V,W;
int e1,e2,e3;