From 393fdaf5467ba71c3c40b2582ed37582f58b692e Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 18 Sep 2021 12:16:42 -0400 Subject: [PATCH 1/7] update doxygen --- analysis/FlowAdaptor.h | 36 ++++++++++++++++++++++++++++++++++++ doxygen/DoxygenMainpage.h | 7 ++++++- 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/analysis/FlowAdaptor.h b/analysis/FlowAdaptor.h index 331c7d7a..d5382339 100644 --- a/analysis/FlowAdaptor.h +++ b/analysis/FlowAdaptor.h @@ -12,12 +12,48 @@ #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); + double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume); double UpdateFractionalFlow(ScaLBL_ColorModel &M); double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil); diff --git a/doxygen/DoxygenMainpage.h b/doxygen/DoxygenMainpage.h index 0d42f215..1e369ff8 100644 --- a/doxygen/DoxygenMainpage.h +++ b/doxygen/DoxygenMainpage.h @@ -3,7 +3,12 @@ * This is the documentation for LBPM * * - \ref IO "IO routines" + * - \ref analysis "Analysis routines" + * - \ref FlowAdaptor "Flow adaptor" + * - \ref models "Lattice Boltzmann models" + * - \ref ScaLBL "Lattice Boltzmann building blocks" + * - \ref tests "Unit tests" * - \ref Utilities "Utility routines" * - * \author James McClure + * \author J.E. McClure, M. Berrill */ From 6eb19ad7644f4b7c6538e3400363d854e13099fc Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 18 Sep 2021 16:32:32 -0400 Subject: [PATCH 2/7] document dcel class --- analysis/FlowAdaptor.h | 29 ++++++++++++++++++++++++++--- analysis/dcel.cpp | 18 +++++++++--------- analysis/dcel.h | 23 ++++++++++++++++------- 3 files changed, 51 insertions(+), 19 deletions(-) diff --git a/analysis/FlowAdaptor.h b/analysis/FlowAdaptor.h index d5382339..89bb03d4 100644 --- a/analysis/FlowAdaptor.h +++ b/analysis/FlowAdaptor.h @@ -47,16 +47,39 @@ public: double MoveInterface(ScaLBL_ColorModel &M); /** - * \brief image re-initialization - * \details Re-initialize LB simulation from image data + * \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; diff --git a/analysis/dcel.cpp b/analysis/dcel.cpp index 91102acf..9e41673d 100644 --- a/analysis/dcel.cpp +++ b/analysis/dcel.cpp @@ -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&Field, const double isovalue) { - DECL object; + DCEL object; int e1,e2,e3; FILE *TRIANGLES; TRIANGLES = fopen("isosurface.stl","w"); diff --git a/analysis/dcel.h b/analysis/dcel.h index 0e6b2b1d..14897fef 100644 --- a/analysis/dcel.h +++ b/analysis/dcel.h @@ -5,7 +5,11 @@ #include "analysis/pmmc.h" /* -Doubly-connected edge list (DECL) +*/ + +/** + * \class Vertex + * @brief store vertex for DCEL data structure */ // Vertex structure @@ -34,8 +38,10 @@ private: }; -// Halfedge structure -// Face +/** + * \class Halfedge + * @brief store half edge for DCEL data structure +*/ class Halfedge{ public: Halfedge() = default; @@ -60,11 +66,14 @@ private: std::vector> 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; From daba5c7521eb9e6d6194865b1f25e1c3aece6b4a Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 18 Sep 2021 16:33:26 -0400 Subject: [PATCH 3/7] extend doxygen support --- analysis/ElectroChemistry.cpp | 1 + analysis/FreeEnergy.h | 8 ++++++ analysis/GreyPhase.h | 16 +++++++++++ analysis/Minkowski.cpp | 2 +- analysis/Minkowski.h | 50 ++++++++++++++++++++++++++++++++++ doxygen/DoxygenMainpage.h | 11 ++++++-- tests/test_dcel_tri_normal.cpp | 2 +- 7 files changed, 85 insertions(+), 5 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index dc932f91..a762fd82 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -1,5 +1,6 @@ #include "analysis/ElectroChemistry.h" + ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm): Dm(dm) { diff --git a/analysis/FreeEnergy.h b/analysis/FreeEnergy.h index fbb1ba31..75e6aef6 100644 --- a/analysis/FreeEnergy.h +++ b/analysis/FreeEnergy.h @@ -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 Dm; diff --git a/analysis/GreyPhase.h b/analysis/GreyPhase.h index 3ab46752..b2632a70 100644 --- a/analysis/GreyPhase.h +++ b/analysis/GreyPhase.h @@ -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 Dm; diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index 9dfff477..5e2e28da 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -50,7 +50,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; diff --git a/analysis/Minkowski.h b/analysis/Minkowski.h index 6cf2673a..e408fceb 100644 --- a/analysis/Minkowski.h +++ b/analysis/Minkowski.h @@ -18,6 +18,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{ //........................................................................... @@ -45,6 +53,7 @@ public: int n_connected_components; //........................................................................... int Nx,Ny,Nz; + double V(){ return Vi; } @@ -59,15 +68,56 @@ public: } //.......................................................................... + /** + * \brief Null constructor + */ Minkowski(){};//NULL CONSTRUCTOR + + /** + * \brief Constructor based on an existing Domain + * @param Dm - Domain structure + */ Minkowski(std::shared_ptr 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(); }; diff --git a/doxygen/DoxygenMainpage.h b/doxygen/DoxygenMainpage.h index 1e369ff8..4f1687b5 100644 --- a/doxygen/DoxygenMainpage.h +++ b/doxygen/DoxygenMainpage.h @@ -1,12 +1,17 @@ /** \mainpage LBPM * - * This is the documentation for LBPM + * C/C++ routines * * - \ref IO "IO routines" * - \ref analysis "Analysis routines" - * - \ref FlowAdaptor "Flow adaptor" + * - \ref FlowAdaptor "FlowAdaptor" + * - \ref DCEL "Doubly connected edge list" + * - \ref Minkowski "Minkowski functionals" * - \ref models "Lattice Boltzmann models" - * - \ref ScaLBL "Lattice Boltzmann building blocks" + * - \ref ScaLBL "Scalable Lattice Boltzmann Library (ScaLBL)" + * - \ref cpu "cpu implementation for library routines" + * - \ref cuda "CUDA implementation for library routines" + * - \ref hip "HIP implementation for library routines" * - \ref tests "Unit tests" * - \ref Utilities "Utility routines" * diff --git a/tests/test_dcel_tri_normal.cpp b/tests/test_dcel_tri_normal.cpp index 622c9fa8..53a42977 100644 --- a/tests/test_dcel_tri_normal.cpp +++ b/tests/test_dcel_tri_normal.cpp @@ -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; From 229feb93fb9d224780b7ba14f4042698eace37e6 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 18 Sep 2021 16:35:51 -0400 Subject: [PATCH 4/7] extend doxygen support --- doxygen/DoxygenMainpage.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/doxygen/DoxygenMainpage.h b/doxygen/DoxygenMainpage.h index 4f1687b5..9b0f1596 100644 --- a/doxygen/DoxygenMainpage.h +++ b/doxygen/DoxygenMainpage.h @@ -9,9 +9,6 @@ * - \ref Minkowski "Minkowski functionals" * - \ref models "Lattice Boltzmann models" * - \ref ScaLBL "Scalable Lattice Boltzmann Library (ScaLBL)" - * - \ref cpu "cpu implementation for library routines" - * - \ref cuda "CUDA implementation for library routines" - * - \ref hip "HIP implementation for library routines" * - \ref tests "Unit tests" * - \ref Utilities "Utility routines" * From 7575d5c116d4c97140db29205bb29b5855301ee1 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 19 Sep 2021 09:52:06 -0400 Subject: [PATCH 5/7] depcrecate internal color flow adaptor --- analysis/dcel.h | 3 - models/ColorModel.cpp | 734 +----------------------------------------- models/ColorModel.h | 66 +++- 3 files changed, 69 insertions(+), 734 deletions(-) diff --git a/analysis/dcel.h b/analysis/dcel.h index 14897fef..e18f6930 100644 --- a/analysis/dcel.h +++ b/analysis/dcel.h @@ -4,9 +4,6 @@ #include #include "analysis/pmmc.h" -/* -*/ - /** * \class Vertex * @brief store vertex for DCEL data structure diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 2f15d320..b94b4b1f 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -921,7 +921,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( "analysis_interval" ); + } + /* int IMAGE_INDEX = 0; int IMAGE_COUNT = 0; std::vector ImageList; @@ -940,7 +944,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; @@ -951,7 +954,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; @@ -1003,9 +1005,7 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "tolerance" )){ tolerance = analysis_db->getScalar( "tolerance" ); } - if (analysis_db->keyExists( "analysis_interval" )){ - analysis_interval = analysis_db->getScalar( "analysis_interval" ); - } + if (analysis_db->keyExists( "min_steady_timesteps" )){ MIN_STEADY_TIMESTEPS = analysis_db->getScalar( "min_steady_timesteps" ); } @@ -1016,7 +1016,6 @@ void ScaLBL_ColorModel::Run(){ MAX_MORPH_TIMESTEPS = analysis_db->getScalar( "max_morph_timesteps" ); } - /* defaults for simulation protocols */ auto protocol = color_db->getWithDefault( "protocol", "none" ); if (protocol == "image sequence"){ // Get the list of images @@ -1062,8 +1061,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"); @@ -1118,7 +1116,7 @@ void ScaLBL_ColorModel::Run(){ printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); } - + */ //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); PROFILE_START("Loop"); @@ -1218,242 +1216,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("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("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("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"); @@ -1477,486 +1239,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; iid[i]; // save what was read - for (int i=0; iid[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; kComm.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 id_solid(nx,ny,nz); - Array 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 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; kSDs(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; kComm.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 phase_id(Nx,Ny,Nz); - fillHalo 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 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; kComm.sumReduce( count); - second_biggest = Dm->Comm.sumReduce( second_biggest); - } - else { - // use the whole NWP - for (int k=0; kSDs(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 phase_distance - CalcDist(phase_distance,phase_id,*Dm); - - double temp,value; - double factor=0.5/beta; - for (int k=0; k 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; kSDs(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 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); diff --git a/models/ColorModel.h b/models/ColorModel.h index a3eddb29..6be45fe1 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -21,21 +21,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 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; @@ -74,6 +131,9 @@ public: double *Velocity; double *Pressure; + /** + * \brief Assign wetting affinity values + */ void AssignComponentLabels(double *phase); private: @@ -88,10 +148,6 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr 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 From 75f31678a02aab587e82fa5aa2a37568d60b0516 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 19 Sep 2021 20:24:40 -0400 Subject: [PATCH 6/7] adding dox --- common/ScaLBL.h | 133 ++++++++++++++++++++++++++++++++++++-- doxygen/Doxyfile.in | 2 +- doxygen/DoxygenMainpage.h | 9 +-- 3 files changed, 134 insertions(+), 10 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 07fb72bf..916bf4bd 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -1,5 +1,5 @@ -/* 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 @@ -11,43 +11,166 @@ #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 diff --git a/doxygen/Doxyfile.in b/doxygen/Doxyfile.in index b654532f..851b36e2 100755 --- a/doxygen/Doxyfile.in +++ b/doxygen/Doxyfile.in @@ -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 diff --git a/doxygen/DoxygenMainpage.h b/doxygen/DoxygenMainpage.h index 9b0f1596..e46ec3c9 100644 --- a/doxygen/DoxygenMainpage.h +++ b/doxygen/DoxygenMainpage.h @@ -2,15 +2,16 @@ * * C/C++ routines * - * - \ref IO "IO 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 models "Lattice Boltzmann models" - * - \ref ScaLBL "Scalable Lattice Boltzmann Library (ScaLBL)" - * - \ref tests "Unit tests" + * - \ref IO "IO routines" * - \ref Utilities "Utility routines" + * - \ref tests "Unit tests" * * \author J.E. McClure, M. Berrill */ From a3a2bf24ddb71ae078f6166eb8e9f823c57a7e21 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 20 Sep 2021 06:29:05 -0400 Subject: [PATCH 7/7] expanding dox for ScaLBL --- common/ScaLBL.h | 222 ++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 207 insertions(+), 15 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 916bf4bd..0ec447a5 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -174,8 +174,31 @@ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np); 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) @@ -243,26 +266,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); @@ -272,27 +342,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); @@ -303,7 +484,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); @@ -454,12 +634,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 Dm); + /** + *\brief Constructor + * @param Dm - Domain information + */ + ScaLBL_Communicator(std::shared_ptr Dm); - //ScaLBL_Communicator(Domain &Dm, IntArray &Map); + /** + *\brief Destructor + */ ~ScaLBL_Communicator(); //...................................................................................... unsigned long int CommunicationCount,SendCount,RecvCount;