diff --git a/CMakeLists.txt b/CMakeLists.txt index 31518d86..c9b7e66a 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,8 +18,7 @@ SET( TEST_MAX_PROCS 16 ) # Initialize the project -PROJECT( ${PROJ} LANGUAGES CXX ) - +PROJECT( ${PROJ} LANGUAGES CXX) # Prevent users from building in place IF ("${CMAKE_CURRENT_SOURCE_DIR}" STREQUAL "${CMAKE_CURRENT_BINARY_DIR}" ) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 2a0b7169..dc932f91 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -53,19 +53,32 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss Poisson.getElectricPotential(ElectricalPotential); /* local sub-domain averages */ - double rho_avg_local[Ion.number_ion_species]; - double rho_mu_avg_local[Ion.number_ion_species]; - double rho_mu_fluctuation_local[Ion.number_ion_species]; - double rho_psi_avg_local[Ion.number_ion_species]; - double rho_psi_fluctuation_local[Ion.number_ion_species]; + double *rho_avg_local; + double *rho_mu_avg_local; + double *rho_mu_fluctuation_local; + double *rho_psi_avg_local; + double *rho_psi_fluctuation_local; /* global averages */ - double rho_avg_global[Ion.number_ion_species]; - double rho_mu_avg_global[Ion.number_ion_species]; - double rho_mu_fluctuation_global[Ion.number_ion_species]; - double rho_psi_avg_global[Ion.number_ion_species]; - double rho_psi_fluctuation_global[Ion.number_ion_species]; + double *rho_avg_global; + double *rho_mu_avg_global; + double *rho_mu_fluctuation_global; + double *rho_psi_avg_global; + double *rho_psi_fluctuation_global; - for (int ion=0; ionrank()==0){ fprintf(TIMELOG,"%i ",timestep); - for (int ion=0; ion( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); auto ElectricPotential = std::make_shared(); std::vector> IonConcentration; - for (int ion=0; ion()); } auto VxVar = std::make_shared(); @@ -163,8 +176,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P } if (vis_db->getWithDefault( "save_concentration", true )){ - for (int ion=0; ionname = VisName; IonConcentration[ion]->type = IO::VariableType::VolumeVariable; IonConcentration[ion]->dim = 1; @@ -199,8 +212,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P } if (vis_db->getWithDefault( "save_concentration", true )){ - for (int ion=0; ionname = VisName; ASSERT(visData[0].vars[1+ion]->name==VisName); Array& IonConcentrationData = visData[0].vars[1+ion]->data; diff --git a/analysis/FreeEnergy.cpp b/analysis/FreeEnergy.cpp index 6a641a95..b49f2ec5 100644 --- a/analysis/FreeEnergy.cpp +++ b/analysis/FreeEnergy.cpp @@ -47,8 +47,6 @@ void FreeEnergyAnalyzer::SetParams(){ void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){ - int i,j,k; - if (Dm->rank()==0){ fprintf(TIMELOG,"%i ",timestep); /*for (int ion=0; ion input_db, int timestep){ auto vis_db = input_db->getDatabase( "Visualization" ); - char VisName[40]; std::vector visData; fillHalo fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1); diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 324c45cb..4b023db6 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -120,6 +120,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr sendtag = recvtag = 7; int ii,jj,kk; + int imin,jmin,kmin,imax,jmax,kmax; int Nx = nx; int Ny = ny; int Nz = nz; @@ -128,17 +129,13 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr double void_fraction_new=1.0; double void_fraction_diff_old = 1.0; double void_fraction_diff_new = 1.0; - - // Increase the critical radius until the target saturation is met - double deltaR=0.05; // amount to change the radius in voxel units - double Rcrit_old; - - int imin,jmin,kmin,imax,jmax,kmax; - if (ErodeLabel == 1){ VoidFraction = 1.0 - VoidFraction; } + // Increase the critical radius until the target saturation is met + double deltaR=0.05; // amount to change the radius in voxel units + double Rcrit_old = maxdistGlobal; double Rcrit_new = maxdistGlobal; while (void_fraction_new > VoidFraction) @@ -406,6 +403,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr2){ // Rcrit_new = strtod(argv[2],NULL); @@ -457,7 +452,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptrComm.dup(); d_Np = ColorModel.Np; - bool Regular = false; auto input_db = ColorModel.db; auto db = input_db->getDatabase( "Analysis" ); diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index ec13efc3..4c956e7f 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1046,29 +1046,28 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in } } } - + int *bb_dist_tmp = new int [local_count]; int *bb_interactions_tmp = new int [local_count]; ScaLBL_AllocateDeviceMemory((void **) &bb_dist, sizeof(int)*local_count); ScaLBL_AllocateDeviceMemory((void **) &bb_interactions, sizeof(int)*local_count); - int *fluid_boundary_tmp; - double *lattice_weight_tmp; - float *lattice_cx_tmp; - float *lattice_cy_tmp; - float *lattice_cz_tmp; - if(SlippingVelBC==true){ - fluid_boundary_tmp = new int [local_count]; - lattice_weight_tmp = new double [local_count]; - lattice_cx_tmp = new float [local_count]; - lattice_cy_tmp = new float [local_count]; - lattice_cz_tmp = new float [local_count]; - ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count); - ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count); - ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count); - ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count); - ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count); - } - + int *fluid_boundary_tmp; + double *lattice_weight_tmp; + float *lattice_cx_tmp; + float *lattice_cy_tmp; + float *lattice_cz_tmp; + /* allocate memory for bounce-back sites */ + fluid_boundary_tmp = new int [local_count]; + lattice_weight_tmp = new double [local_count]; + lattice_cx_tmp = new float [local_count]; + lattice_cy_tmp = new float [local_count]; + lattice_cz_tmp = new float [local_count]; + ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count); + local_count=0; for (k=1;k= 1.3 + else: + html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] + html_theme = 'sphinx_rtd_theme' + +#def setup(app): +# app.add_stylesheet("fix_rtd.css") diff --git a/docs/source/developerGuide/buildingModels/overview.rst b/docs/source/developerGuide/buildingModels/overview.rst new file mode 100644 index 00000000..223b7c99 --- /dev/null +++ b/docs/source/developerGuide/buildingModels/overview.rst @@ -0,0 +1,20 @@ +=========================== +Implementing a new LB model +=========================== + +While LBPM includes a range of fully-functioning lattice Boltzmann models, the commonly used +Bhatnager-Gross-Krook (BGK) model has been deliberately excluded. While the physical limitations +of this model are well-known, implementing the BGK model is an excellent way to understand +how to implement new LB models within the more general framework of LBPM. In this excercise +you will + +* learn "what goes where" + + + +* don't modify core data structures (unless you have a really good reason) + +* re-use existing components whenever possible + + + diff --git a/docs/source/developerGuide/designOverview/dataStructures.rst b/docs/source/developerGuide/designOverview/dataStructures.rst new file mode 100644 index 00000000..bea4f91b --- /dev/null +++ b/docs/source/developerGuide/designOverview/dataStructures.rst @@ -0,0 +1,6 @@ +=============== +Data Structures +=============== + +LBPM includes a variety of generalized data structures to facilitate the implementation +of different lattice Boltzmann models. diff --git a/docs/source/developerGuide/index.rst b/docs/source/developerGuide/index.rst new file mode 100644 index 00000000..d8448fa1 --- /dev/null +++ b/docs/source/developerGuide/index.rst @@ -0,0 +1,18 @@ +############################################################################### +Developer guide +############################################################################### + +The LBPM developer guide provides essential information on how to add new physics +into the framework. + +.. toctree:: + :glob: + :maxdepth: 2 + + designOverview/* + + buildingModels/* + + testingModels/* + + diff --git a/docs/source/developerGuide/testingModels/unitTests.rst b/docs/source/developerGuide/testingModels/unitTests.rst new file mode 100644 index 00000000..4a4e84bf --- /dev/null +++ b/docs/source/developerGuide/testingModels/unitTests.rst @@ -0,0 +1,9 @@ +================= +Adding unit tests +================= + +Unit tests in LBPM are implemented using ctest + +* general overview + +* launching unit tests for GPU (MPI flags etc.) diff --git a/docs/source/examples/running.rst b/docs/source/examples/running.rst new file mode 100644 index 00000000..83993f59 --- /dev/null +++ b/docs/source/examples/running.rst @@ -0,0 +1,9 @@ +============ +Running LBPM +============ + +There are two main components to running LBPM simulators. +First is understanding how to launch MPI tasks on your system, +which depends on the particular implementation of MPI that you are using, +as well as other details of the local configuration. The second component is +understanding the LBPM input file structure. diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 00000000..fc5225d1 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,25 @@ +.. LBPM documentation master file, created by + sphinx-quickstart on Thu May 20 12:19:14 2021. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +LBPM -- Documentation +=================================================== + +.. toctree:: + :glob: + :maxdepth: 2 + :caption: Contents: + + install + examples/* + userGuide/* + developerGuide/* + publications/* + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/source/install.rst b/docs/source/install.rst new file mode 100644 index 00000000..232da8d0 --- /dev/null +++ b/docs/source/install.rst @@ -0,0 +1,10 @@ +============ +Installation +============ + +LBPM requires CMake, MPI and HDF5 as required dependencies. + +.. code-block:: bash + + ls $LBPM_SOURCE/sample_scripts + diff --git a/docs/source/publications/publications.rst b/docs/source/publications/publications.rst new file mode 100644 index 00000000..aea8e5b6 --- /dev/null +++ b/docs/source/publications/publications.rst @@ -0,0 +1,13 @@ +============ +Publications +============ + +* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9 + + +* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670 + + +* Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966 + + diff --git a/docs/source/userGuide/IO/fileformat.rst b/docs/source/userGuide/IO/fileformat.rst new file mode 100644 index 00000000..faec04a6 --- /dev/null +++ b/docs/source/userGuide/IO/fileformat.rst @@ -0,0 +1,13 @@ +======================== +I/O conventions for LBPM +======================== + +There are three main kinds of output file that are supported by LBPM. + + +* CSV files -- + +* formatted binary files -- + +* unformatted binary files -- + diff --git a/docs/source/userGuide/analysis/analysis.rst b/docs/source/userGuide/analysis/analysis.rst new file mode 100644 index 00000000..10d5282a --- /dev/null +++ b/docs/source/userGuide/analysis/analysis.rst @@ -0,0 +1,5 @@ +=========================== +Internal Analysis Framework +=========================== + +placeholder for analysis diff --git a/docs/source/userGuide/index.rst b/docs/source/userGuide/index.rst new file mode 100644 index 00000000..49025202 --- /dev/null +++ b/docs/source/userGuide/index.rst @@ -0,0 +1,17 @@ +############################################################################### +User Guide +############################################################################### + +Welcome to the LBPM user guide. + +.. toctree:: + :glob: + :maxdepth: 2 + + models/* + + analysis/* + + visualization/* + + IO/* diff --git a/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst b/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst new file mode 100644 index 00000000..93dbaff0 --- /dev/null +++ b/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst @@ -0,0 +1,6 @@ +============================================= +Poisson-Boltzmann model +============================================= + +The LBPM Poisson-Boltzmann solver is designed to solve the Poisson-Boltzmann equation +to solve for the electric field in an ionic fluid. diff --git a/docs/source/userGuide/models/color/index.rst b/docs/source/userGuide/models/color/index.rst new file mode 100644 index 00000000..23bb5949 --- /dev/null +++ b/docs/source/userGuide/models/color/index.rst @@ -0,0 +1,91 @@ +############################################################################### +Color model +############################################################################### + +The LBPM color model is implemented by combining a multi-relaxation time D3Q19 +lattice Boltzmann equation (LBE) to solve for the momentum transport with two D3Q7 +LBEs for the mass transport. The color model will obey strict mass and momentum +conservation while minimizing diffusive fluxes across the interface between fluids. +The color model is a good choice for modeling dense fluids that are strongly immiscible +(e.g. water-oil systems). Due to the strong anti-diffusion in the interface region, +the color model is not suitable for modeling processes such as Ostwald ripening that +depend on diffusive fluxes between fluid phases. + +A typical command to launch the LBPM color simulator is as follows + +``` +mpirun -np $NUMPROCS lbpm_color_simulator input.db +``` + +where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is +the name of the input database that provides the simulation parameters. +Note that the specific syntax to launch MPI tasks may vary depending on your system. +For additional details please refer to your local system documentation. + +**************************** +Simulation protocols +**************************** + +Simulation protocols are designed to make it simpler to design and execute common +computational experiments. Protocols will automatically determine boundary conditions +needed to perform a particular simulation. LBPM will internall set default simulation paramaters +that can be over-ridden to develop customized simulations. + +.. toctree:: + :glob: + :maxdepth: 2 + + protocols/* + + +*************************** +Model parameters +*************************** + +The essential model parameters for the color model are + +- :math:`\alpha` -- control the interfacial tension between fluids with key ``alpha`` +- :math:`\beta` -- control the width of the interface with key ``beta`` +- :math:`\tau_A` -- control the viscosity of fluid A with key ``tauA`` +- :math:`\tau_B` -- control the viscosity of fluid B with key ``tauB`` + +**************************** +Model Formulation +**************************** + + + +**************************** +Boundary Conditions +**************************** + +The following external boundary conditions are supported by ``lbpm_color_simulator`` +and can be set by setting the ``BC`` key values in the ``Domain`` section of the +input file database + +- ``BC = 0`` -- fully periodic boundary conditions +- ``BC = 3`` -- constant pressure boundary condition +- ``BC = 4`` -- constant volumetric flux boundary condition + +For ``BC = 0`` any mass that exits on one side of the domain will re-enter at the other +side. If the pore-structure for the image is tight, the mismatch between the inlet and +outlet can artificially reduce the permeability of the sample due to the blockage of +flow pathways at the boundary. LBPM includes an internal utility that will reduce the impact +of the boundary mismatch by eroding the solid labels within the inlet and outlet layers +(https://doi.org/10.1007/s10596-020-10028-9) to create a mixing layer. +The number mixing layers to use can be set using the key values in the ``Domain`` section +of the input database + +- ``InletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the inlet +- ``OUtletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the outlet + +For the other boundary conditions a thin reservoir of fluid (default ``3`` voxels) +is established at either side of the domain. The inlet is defined as the boundary face +where ``z = 0`` and the outlet is the boundary face where ``z = nprocz*nz``. By default a +reservoir of fluid A is established at the inlet and a reservoir of fluid B is established at +the outlet, each with a default thickness of three voxels. To over-ride the default label at +the inlet or outlet, the ``Domain`` section of the database may specify the following key values + +- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet +- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet + diff --git a/docs/source/userGuide/models/color/protocols/centrifuge.rst b/docs/source/userGuide/models/color/protocols/centrifuge.rst new file mode 100644 index 00000000..21cf4eee --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/centrifuge.rst @@ -0,0 +1,12 @@ +====================================== +Color model -- Centrifuge Protocol +====================================== + +The centrifuge protocol is designed to mimic SCAL centrifuge experiments that +are used to infer the capillary pressure. The centrifuge protocol + +.. code-block:: bash + + image_sequence = "image1.raw", "image2.raw" + + diff --git a/docs/source/userGuide/models/color/protocols/coreFlooding.rst b/docs/source/userGuide/models/color/protocols/coreFlooding.rst new file mode 100644 index 00000000..12d1e799 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/coreFlooding.rst @@ -0,0 +1,63 @@ +====================================== +Color model -- Core Flooding +====================================== + +The core flooding protocol is designed to mimic SCAL experiments where one +immiscible fluid is injected into the sample at a constant rate, displacing the +other fluid. The core flooding protocol relies on a flux boundary condition +to ensure that fluid is injected into the sample at a constant rate. The flux +boundary condition implements a time-varying pressure boundary condition that +adapts to ensure a constant volumetric flux. Details for the flux boundary +condition are available +(see: https://doi.org/10.1016/j.compfluid.2020.104670) + +.. code-block:: bash + + protocol = "core flooding" + + +To match experimental conditions, it is usually important to match the capillary +number, which is + +.. math:: + \mbox{Ca} = \frac{\mu u_z}{\gamma} + + +where :math:`\mu` is the dynamic viscosity, :math:`u_z` is the fluid +(usually water) velocity and :math:`\gamma` is the interfacial tension. +The volumetric flow rate is related to the fluid velocity based on + +.. math:: + Q_z = \epsilon C_{xy} u_z + +where :math:`C_{xy}` is the cross-sectional area and :math:`\epsilon` +is the porosity. Given a particular experimental system +self-similar conditions can be determined for the lattice Boltzmann +simulation by matching the non-dimensional :math:`mbox{Ca}`. It is nearly +awlays advantageous for the timestep to be as large as possible so +that time-to-solution will be more favorable. This is accomplished by + +* use a high value for the numerical surface tension (e.g. ``alpha=1.0e-2``) +* use a small value for the fluid viscosity (e.g. ``tau_w = 0.7`` and ``tau_n = 0.7`` ) +* determine the volumetric flow rate needed to match :math:`\mbox{Ca}` + +For the color LBM the interfacial tension is +:math:`\gamma = 6 \alpha` and the dynamic viscosity is :math:`\mu = \rho(\tau-1/2)/3`, +where the units are relative to the lattice spacing, timestep and mass +density. Agreemetn between the experimental and simulated values for +:math:`\mbox{Ca}` is ensured by setting the volumetric flux + +.. math:: + Q_z = \frac{\epsilon C_{xy} \gamma }{\mu} \mbox{Ca} + +where the LB units of the volumetric flux will be voxels per timestep. + +In some situations it may also be important to match other non-dimensional numbers, +such as the viscosity ratio, density ratio, and/or Ohnesorge/Laplace number. This +can be accomplished with an analogous procedure. Enforcing additional constraints +will necessarily restrict the LB parameters that can be used, which are ultimately +manifested as a constraint on the size of a timestep. + + + + diff --git a/docs/source/userGuide/models/color/protocols/fractionalFlow.rst b/docs/source/userGuide/models/color/protocols/fractionalFlow.rst new file mode 100644 index 00000000..7f5a26f0 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/fractionalFlow.rst @@ -0,0 +1,17 @@ +========================================== +Color model -- Fractional Flow Protocol +========================================== + +The fractional flow protocol is designed to perform steady-state relative +permeability simulations by using an internal routine to change the fluid +saturation by adding and subtracting mass to the fluid phases. The +mass density is updated for each fluid based on the locations where +the local rate of flow is high. + +.. code-block:: bash + + protocol = "fractional flow" + + + + diff --git a/docs/source/userGuide/models/color/protocols/imageSequence.rst b/docs/source/userGuide/models/color/protocols/imageSequence.rst new file mode 100644 index 00000000..25c91e6a --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/imageSequence.rst @@ -0,0 +1,27 @@ +====================================== +Color model -- Image Sequence Protocol +====================================== + +The image sequence protocol is designed to perform a set steady-state +simulations based on a sequence of 3D (8-bit) images provided by the user. +The images might be the output of a previous LBPM simulation, a sequence of +(segmented) experimental data, or data generated from a custom routine. +The image sequence protocol will apply the same set of flow conditions +to all images in the sequence. This means + +* the image labels and any associated properties are the same +* the external boundary conditions are the same +* the physical simulation parameters are the same + +The image sequence protocol does not set boundary conditions by default. +It is up to the user to determine the flow condition, with the understanding +that the same set of will be applied to each image in the sequence. + +To enable the image sequence protocol, the following keys should be set +within the ``Color`` section of the input database + +.. code-block:: bash + + protocol = "image sequence" + image_sequence = "image1.raw", "image2.raw" + diff --git a/docs/source/userGuide/models/color/protocols/shellAggregation.rst b/docs/source/userGuide/models/color/protocols/shellAggregation.rst new file mode 100644 index 00000000..1ab54131 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/shellAggregation.rst @@ -0,0 +1,18 @@ +========================================== +Color model -- Shell Aggregation Protocol +========================================== + +The shell aggregation protocol is designed to perform steady-state relative +permeability simulations by using an internal routine to change the fluid +saturation by moving the interface. The basic design for the shell aggregation +protocol is described by Wang et al. ( https://doi.org/10.1016/j.jcp.2019.108966 ). + + +.. code-block:: bash + + + protocol = "shell aggregation" + + + + diff --git a/docs/source/userGuide/models/freeEnergy/freeEnergy.rst b/docs/source/userGuide/models/freeEnergy/freeEnergy.rst new file mode 100644 index 00000000..f53cfaad --- /dev/null +++ b/docs/source/userGuide/models/freeEnergy/freeEnergy.rst @@ -0,0 +1,6 @@ +============================================= +Free energy model +============================================= + +The LBPM free energy model is constructed to solve the Allen-Cahn equations, +which are typically used to describe liquid-gas systems. diff --git a/docs/source/userGuide/models/greyscale/greyscale.rst b/docs/source/userGuide/models/greyscale/greyscale.rst new file mode 100644 index 00000000..47c0621f --- /dev/null +++ b/docs/source/userGuide/models/greyscale/greyscale.rst @@ -0,0 +1,7 @@ +============================================= +Greyscale model model +============================================= + +The LBPM greyscale lattice Boltzmann model is constructed to approximate the +solution of the Darcy-Brinkman equations in grey regions, coupled to a Navier-Stokes +solution in open regions. diff --git a/docs/source/userGuide/models/index.rst b/docs/source/userGuide/models/index.rst new file mode 100644 index 00000000..6a158aa6 --- /dev/null +++ b/docs/source/userGuide/models/index.rst @@ -0,0 +1,23 @@ +############################################################################### +LBPM model summary +############################################################################### + +Currently supported lattice Boltzmann models + +.. toctree:: + :glob: + :maxdepth: 2 + + color/* + + mrt/* + + nernstPlanck/* + + PoissonBoltzmann/* + + greyscale/* + + freeEnergy/* + + diff --git a/docs/source/userGuide/models/mrt/mrt.rst b/docs/source/userGuide/models/mrt/mrt.rst new file mode 100644 index 00000000..1e232ea6 --- /dev/null +++ b/docs/source/userGuide/models/mrt/mrt.rst @@ -0,0 +1,6 @@ +============================================= +MRT model +============================================= + +The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the +solution of the Navier-Stokes equations. diff --git a/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst b/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst new file mode 100644 index 00000000..3842e2e6 --- /dev/null +++ b/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst @@ -0,0 +1,6 @@ +============================================= +Nernst-Planck model +============================================= + +The Nernst-Planck model is designed to model ion transport based on the +Nernst-Planck equation. diff --git a/docs/source/userGuide/visualization/visit.rst b/docs/source/userGuide/visualization/visit.rst new file mode 100644 index 00000000..f9ff2c8a --- /dev/null +++ b/docs/source/userGuide/visualization/visit.rst @@ -0,0 +1,5 @@ +====================================== +Visualizing simulation data with visit +====================================== + +placeholder for visit diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index cf4b1e13..9d96a667 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -145,6 +145,20 @@ void ScaLBL_ColorModel::ReadParams(string filename){ else if (domain_db->keyExists( "BC" )){ BoundaryCondition = domain_db->getScalar( "BC" ); } + if (domain_db->keyExists( "InletLayersPhase" )){ + int inlet_layers_phase = domain_db->getScalar( "InletLayersPhase" ); + if (inlet_layers_phase == 2 ) { + inletA = 0.0; + inletB = 1.0; + } + } + if (domain_db->keyExists( "OutletLayersPhase" )){ + int outlet_layers_phase = domain_db->getScalar( "OutletLayersPhase" ); + if (outlet_layers_phase == 1 ) { + inletA = 1.0; + inletB = 0.0; + } + } // Override user-specified boundary condition for specific protocols auto protocol = color_db->getWithDefault( "protocol", "none" ); @@ -176,6 +190,29 @@ void ScaLBL_ColorModel::ReadParams(string filename){ } domain_db->putScalar( "BC", BoundaryCondition ); } + else if (protocol == "centrifuge"){ + if (BoundaryCondition != 3 ){ + BoundaryCondition = 3; + if (rank==0) printf("WARNING: protocol (centrifuge) supports only constant pressure boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } + else if (protocol == "core flooding"){ + if (BoundaryCondition != 4){ + BoundaryCondition = 4; + if (rank==0) printf("WARNING: protocol (core flooding) supports only volumetric flux boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + if (color_db->keyExists( "capillary_number" )){ + double capillary_number = color_db->getScalar( "capillary_number" ); + double MuB = rhoB*(tauB - 0.5)/3.0; + double IFT = 6.0*alpha; + double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2)); + flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB; + if (rank==0) printf(" protocol (core flooding): set flux=%f to achieve Ca=%f \n",flux, capillary_number); + } + color_db->putScalar( "flux", flux ); + } } void ScaLBL_ColorModel::SetDomain(){ @@ -306,11 +343,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) if (WettingConvention == "SCAL"){ for (size_t idx=0; idxkeyExists( "tolerance" )){ tolerance = analysis_db->getScalar( "tolerance" ); } - + runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); auto t1 = std::chrono::system_clock::now(); int CURRENT_TIMESTEP = 0; @@ -878,32 +916,7 @@ void ScaLBL_ColorModel::Run(){ if (color_db->keyExists( "krA_morph_factor" )){ KRA_MORPH_FACTOR = color_db->getScalar( "krA_morph_factor" ); } - /* defaults for simulation protocols */ - auto protocol = color_db->getWithDefault( "protocol", "none" ); - if (protocol == "image sequence"){ - // Get the list of images - USE_DIRECT = true; - ImageList = color_db->getVector( "image_sequence"); - IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); - IMAGE_COUNT = ImageList.size(); - morph_interval = 10000; - USE_MORPH = true; - } - else if (protocol == "seed water"){ - morph_delta = -0.05; - seed_water = 0.01; - USE_SEED = true; - USE_MORPH = true; - } - else if (protocol == "open connected oil"){ - morph_delta = -0.05; - USE_MORPH = true; - USE_MORPHOPEN_OIL = true; - } - else if (protocol == "shell aggregation"){ - morph_delta = -0.05; - USE_MORPH = true; - } + if (color_db->keyExists( "capillary_number" )){ capillary_number = color_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; @@ -922,7 +935,6 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "seed_water" )){ seed_water = analysis_db->getScalar( "seed_water" ); if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water); - ASSERT(protocol == "seed water"); } if (analysis_db->keyExists( "morph_delta" )){ morph_delta = analysis_db->getScalar( "morph_delta" ); @@ -952,7 +964,54 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "max_morph_timesteps" )){ 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 + USE_DIRECT = true; + ImageList = color_db->getVector( "image_sequence"); + IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); + IMAGE_COUNT = ImageList.size(); + morph_interval = 10000; + USE_MORPH = true; + USE_SEED = false; + } + else if (protocol == "seed water"){ + morph_delta = -0.05; + seed_water = 0.01; + USE_SEED = true; + USE_MORPH = true; + } + else if (protocol == "open connected oil"){ + morph_delta = -0.05; + USE_SEED = false; + USE_MORPH = true; + USE_MORPHOPEN_OIL = true; + } + else if (protocol == "shell aggregation"){ + morph_delta = -0.05; + USE_MORPH = true; + USE_SEED = false; + } + else if (protocol == "fractional flow"){ + USE_MORPH = false; + USE_SEED = false; + } + else if (protocol == "centrifuge"){ + USE_MORPH = false; + USE_SEED = false; + } + else if (protocol == "core flooding"){ + USE_MORPH = false; + USE_SEED = false; + if (SET_CAPILLARY_NUMBER){ + double MuB = rhoB*(tauB - 0.5)/3.0; + double IFT = 6.0*alpha; + double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2)); + flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB; + } + } if (rank==0){ printf("********************************************************\n"); if (protocol == "image sequence"){ @@ -990,7 +1049,20 @@ void ScaLBL_ColorModel::Run(){ printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); printf(" tolerance = %f \n",tolerance); - printf(" morph_delta = %f \n",morph_delta); + } + else if (protocol == "centrifuge"){ + printf(" using protocol = centrifuge \n"); + printf(" driving force = %f \n",Fz); + if (Fz < 0){ + printf(" Component B displacing component A \n"); + } + else if (Fz > 0){ + printf(" Component A displacing component B \n"); + } + } + else if (protocol == "core flooding"){ + printf(" using protocol = core flooding \n"); + printf(" capillary number = %f \n", capillary_number); } printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); @@ -1104,7 +1176,7 @@ void ScaLBL_ColorModel::Run(){ double volB = Averages->gwb.V; double volA = Averages->gnb.V; volA /= Dm->Volume; - volB /= 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; @@ -1917,291 +1989,565 @@ FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M){ Nz = M.Dm->Nz; timestep=-1; timestep_previous=-1; - + phi.resize(Nx,Ny,Nz); phi.fill(0); // phase indicator field phi_t.resize(Nx,Ny,Nz); phi_t.fill(0); // time derivative for the phase indicator field } FlowAdaptor::~FlowAdaptor(){ - + } -double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ +double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename){ + int rank = M.rank; + int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz; + if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str()); + M.Mask->Decomp(Filename); + for (int i=0; iid[i]; // save what was read + for (int i=0; iid[i] = M.Mask->id[i]; // save what was read + + double *PhaseLabel; + PhaseLabel = new double[Nx*Ny*Nz]; + M.AssignComponentLabels(PhaseLabel); + + double Count = 0.0; + double PoreCount = 0.0; + for (int k=1; kComm.sumReduce( Count); + PoreCount=M.Dm->Comm.sumReduce( PoreCount); - double MASS_FRACTION_CHANGE = 0.05; - if (M.db->keyExists( "FlowAdaptor" )){ - auto flow_db = M.db->getDatabase( "FlowAdaptor" ); - MASS_FRACTION_CHANGE = flow_db->getWithDefault( "fractional_flow_increment", 0.05); - } + if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount); + ScaLBL_CopyToDevice(M.Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double)); + M.Dm->Comm.barrier(); + + ScaLBL_D3Q19_Init(M.fq, M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np); + M.Dm->Comm.barrier(); + + ScaLBL_CopyToHost(M.Averages->Phi.data(),M.Phi,Nx*Ny*Nz*sizeof(double)); - int Np = M.Np; - double dA, dB, phi; - double vx,vy,vz; - double vax,vay,vaz; - double vbx,vby,vbz; - double vax_global,vay_global,vaz_global; - double vbx_global,vby_global,vbz_global; + double saturation = Count/PoreCount; + return saturation; +} - double mass_a, mass_b, mass_a_global, mass_b_global; - - double *Aq_tmp, *Bq_tmp; - double *Vel_x, *Vel_y, *Vel_z, *Phase; - - Aq_tmp = new double [7*Np]; - Bq_tmp = new double [7*Np]; - Phase = new double [Np]; - Vel_x = new double [Np]; - Vel_y = new double [Np]; - Vel_z = new double [Np]; - ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); - ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); - ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double)); - ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double)); - ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double)); +double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ - /* compute the total momentum */ - vax = vay = vaz = 0.0; - vbx = vby = vbz = 0.0; - mass_a = mass_b = 0.0; - for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ - 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]; - 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]; - phi = (dA - dB) / (dA + dB); - Phase[n] = phi; - mass_a += dA; - mass_b += dB; - if (phi > 0.0){ - vax += Vel_x[n]; - vay += Vel_y[n]; - vaz += Vel_z[n]; - } - else { - vbx += Vel_x[n]; - vby += Vel_y[n]; - vbz += Vel_z[n]; - } - } - for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ - 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]; - 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]; - phi = (dA - dB) / (dA + dB); - Phase[n] = phi; - mass_a += dA; - mass_b += dB; - if (phi > 0.0){ - vax += Vel_x[n]; - vay += Vel_y[n]; - vaz += Vel_z[n]; - } - else { - vbx += Vel_x[n]; - vby += Vel_y[n]; - vbz += Vel_z[n]; - } - } - mass_a_global = M.Dm->Comm.sumReduce(mass_a); - mass_b_global = M.Dm->Comm.sumReduce(mass_b); - vax_global = M.Dm->Comm.sumReduce(vax); - vay_global = M.Dm->Comm.sumReduce(vay); - vaz_global = M.Dm->Comm.sumReduce(vaz); - vbx_global = M.Dm->Comm.sumReduce(vbx); - vby_global = M.Dm->Comm.sumReduce(vby); - vbz_global = M.Dm->Comm.sumReduce(vbz); - - double total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global); - double total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global); - - /* compute the total mass change */ - double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global); - if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global ) - TOTAL_MASS_CHANGE = 0.1*mass_a_global; - if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global ) - TOTAL_MASS_CHANGE = 0.1*mass_b_global; - - double LOCAL_MASS_CHANGE = 0.0; - for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ - phi = Phase[n]; - vx = Vel_x[n]; - vy = Vel_y[n]; - vz = Vel_z[n]; - double local_momentum = sqrt(vx*vx+vy*vy+vz*vz); - if (phi > 0.0){ - LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A; - Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE; - Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - } - else{ - LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B; - Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE; - Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - } - } + double MASS_FRACTION_CHANGE = 0.0005; + if (M.db->keyExists( "FlowAdaptor" )){ + auto flow_db = M.db->getDatabase( "FlowAdaptor" ); + MASS_FRACTION_CHANGE = flow_db->getWithDefault( "mass_fraction_factor", 0.0005); + } + int Np = M.Np; + double dA, dB, phi; + double vx,vy,vz; + double vax,vay,vaz; + double vbx,vby,vbz; + double vax_global,vay_global,vaz_global; + double vbx_global,vby_global,vbz_global; - for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ - phi = Phase[n]; - vx = Vel_x[n]; - vy = Vel_y[n]; - vz = Vel_z[n]; - double local_momentum = sqrt(vx*vx+vy*vy+vz*vz); - if (phi > 0.0){ - LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A; - Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE; - Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; - } - else{ - LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B; - Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE; - Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; - } - } + double mass_a, mass_b, mass_a_global, mass_b_global; - if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global); + double *Aq_tmp, *Bq_tmp; + double *Vel_x, *Vel_y, *Vel_z, *Phase; - // Need to initialize Aq, Bq, Den, Phi directly - //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double)); - ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); - ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; + Phase = new double [Np]; + Vel_x = new double [Np]; + Vel_y = new double [Np]; + Vel_z = new double [Np]; - return(TOTAL_MASS_CHANGE); + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double)); + + int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz; + + /* compute the total momentum */ + vax = vay = vaz = 0.0; + vbx = vby = vbz = 0.0; + mass_a = mass_b = 0.0; + double maxSpeed = 0.0; + double localMaxSpeed = 0.0; + for (int k=1; kSDs(i,j,k); + if (!(n<0) ){ + 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]; + 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]; + phi = (dA - dB) / (dA + dB); + Phase[n] = phi; + mass_a += dA; + mass_b += dB; + if (phi > 0.0){ + vax += Vel_x[n]; + vay += Vel_y[n]; + vaz += Vel_z[n]; + } + else { + vbx += Vel_x[n]; + vby += Vel_y[n]; + vbz += Vel_z[n]; + } + double speed = sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz); + if (distance > 3.0 && speed > localMaxSpeed){ + localMaxSpeed = speed; + } + } + } + } + } + maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed); + mass_a_global = M.Dm->Comm.sumReduce(mass_a); + mass_b_global = M.Dm->Comm.sumReduce(mass_b); + vax_global = M.Dm->Comm.sumReduce(vax); + vay_global = M.Dm->Comm.sumReduce(vay); + vaz_global = M.Dm->Comm.sumReduce(vaz); + vbx_global = M.Dm->Comm.sumReduce(vbx); + vby_global = M.Dm->Comm.sumReduce(vby); + vbz_global = M.Dm->Comm.sumReduce(vbz); + + double total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global); + double total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global); + /* compute the total mass change */ + double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global); + if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global ) + TOTAL_MASS_CHANGE = 0.1*mass_a_global; + if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global ) + TOTAL_MASS_CHANGE = 0.1*mass_b_global; + + double LOCAL_MASS_CHANGE = 0.0; + for (int k=1; k maxSpeed) local_momentum = maxSpeed; + if (phi > 0.0){ + LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A; + Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE; + Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + //DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE; + } + else{ + LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B; + Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE; + Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + //DebugMassB[n] = LOCAL_MASS_CHANGE; + } + } + } + } + } + + if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global); + + // Need to initialize Aq, Bq, Den, Phi directly + //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + + return(TOTAL_MASS_CHANGE); } void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){ - - int Np = M.Np; - double dA, dB, phi; - double mass_a, mass_b, mass_a_global, mass_b_global; - - double *Aq_tmp, *Bq_tmp; - double *Vel_x, *Vel_y, *Vel_z, *Phase; - - Aq_tmp = new double [7*Np]; - Bq_tmp = new double [7*Np]; + int Np = M.Np; + double dA, dB; - ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); - ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + double *Aq_tmp, *Bq_tmp; + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; - for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ - 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]; - 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]; - if (dA > 1.0){ - double mass_change = dA - 1.0; - Aq_tmp[n] -= 0.333333333333333*mass_change; - Aq_tmp[n+Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - if (dB > 1.0){ - double mass_change = dB - 1.0; - Bq_tmp[n] -= 0.333333333333333*mass_change; - Bq_tmp[n+Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - } - for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ - 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]; - 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]; - if (dA > 1.0){ - double mass_change = dA - 1.0; - Aq_tmp[n] -= 0.333333333333333*mass_change; - Aq_tmp[n+Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - if (dB > 1.0){ - double mass_change = dB - 1.0; - Bq_tmp[n] -= 0.333333333333333*mass_change; - Bq_tmp[n+Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - } + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); - ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); - ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ + 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]; + 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]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ + 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]; + 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]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); } double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){ - - double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.1 ); - double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); - ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); - /* compute the local derivative of phase indicator field */ - double beta = M.beta; - double factor = 0.5/beta; - double total_interface_displacement = 0.0; - double total_interface_sites = 0.0; - for (int n=0; nPhi(n); - double dist1 = factor*log((1.0+value1)/(1.0-value1)); - double value2 = phi(n); - double dist2 = factor*log((1.0+value2)/(1.0-value2)); - phi_t(n) = value2; - if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ - /* time derivative of distance */ - double dxdt = 0.125*(dist2-dist1); - /* extrapolate to move the distance further */ - double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; - /* compute the new phase interface */ - phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); - total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); - total_interface_sites += 1.0; - } - } - ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* 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); + double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.1 ); + double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); + + ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); + /* compute the local derivative of phase indicator field */ + double beta = M.beta; + double factor = 0.5/beta; + double total_interface_displacement = 0.0; + double total_interface_sites = 0.0; + for (int n=0; nPhi(n); + double dist1 = factor*log((1.0+value1)/(1.0-value1)); + double value2 = phi(n); + double dist2 = factor*log((1.0+value2)/(1.0-value2)); + phi_t(n) = value2; + if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ + /* time derivative of distance */ + double dxdt = 0.125*(dist2-dist1); + /* extrapolate to move the distance further */ + double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; + /* compute the new phase interface */ + phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); + total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); + total_interface_sites += 1.0; } } - */ + ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) ); + return total_interface_sites; } +double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){ + + const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz); + auto rank = M.rank; + auto Nx = M.Nx; auto Ny = M.Ny; auto Nz = M.Nz; + auto N = Nx*Ny*Nz; + 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(M.Dm->Comm,M.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(), M.Phi, N*sizeof(double)); + + double count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f) count+=1.f; + } + } + } + double volume_initial = M.Dm->Comm.sumReduce( count); + double PoreVolume = M.Dm->Volume*M.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; + } + + // 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,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm); + M.Dm->Comm.barrier(); + + // only operate on component "0"ScaLBL_ColorModel &M, + count = 0.0; + + for (int k=0; kComm.sumReduce( count); + second_biggest = M.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; + } + } + } + } + } + + // 3. Generate a distance map to the largest object -> phase_distance + CalcDist(phase_distance,phase_id,*M.Dm); + + double temp,value; + double factor=0.5/M.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 && M.Averages->SDs(i,j,k) > 1.f ){ + phase_distance(i,j,k) = temp; + } + // erase the original object + phase(i,j,k) = -1.0; + } + } + } + } + + + 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(M.Averages->SDs,phase_distance,phase_id,M.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*M.beta*d))/(1.f+exp(-2.f*M.beta*d))-1.f); + } + } + } + } + } + fillDouble.fill(phase); + + count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f){ + count+=1.f; + } + } + } + } + double volume_final= M.Dm->Comm.sumReduce( count); + + delta_volume = (volume_final-volume_initial); + if (rank == 0) printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial); + if (rank == 0) printf(" new saturation = %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs))); + + // 6. copy back to the device + //if (rank==0) printf("MorphInit: copy data back to device\n"); + ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double)); + + // 7. Re-initialize phase field and density + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np); + auto BoundaryCondition = M.BoundaryCondition; + if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ + if (M.Dm->kproc()==0){ + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2); + } + if (M.Dm->kproc() == M.nprocz-1){ + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3); + } + } + return delta_volume; +} + + +double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){ + srand(time(NULL)); + auto rank = M.rank; + auto Np = M.Np; + 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, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + + + for (int n=0; n < M.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=M.ScaLBL_Comm->FirstInterior(); n < M.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= M.Dm->Comm.sumReduce( count); + mass_loss= M.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(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + + return(mass_loss); +} diff --git a/models/ColorModel.h b/models/ColorModel.h index f942e26c..d69e3980 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -72,6 +72,8 @@ public: double *ColorGrad; double *Velocity; double *Pressure; + + void AssignComponentLabels(double *phase); private: Utilities::MPI comm; @@ -85,7 +87,6 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); - void AssignComponentLabels(double *phase); double ImageInit(std::string filename); double MorphInit(const double beta, const double morph_delta); double SeedPhaseField(const double seed_water_in_oil); @@ -97,7 +98,10 @@ public: FlowAdaptor(ScaLBL_ColorModel &M); ~FlowAdaptor(); double MoveInterface(ScaLBL_ColorModel &M); + 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); void Flatten(ScaLBL_ColorModel &M); DoubleArray phi; DoubleArray phi_t; diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index b1c0e99f..19a2afa8 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -63,7 +63,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray // Gets data (in optimized layout) from the HOST and stores in regular layout // Primarly for debugging int i,j,k,idx; - int n; // initialize the array regdata.fill(0.f); @@ -72,7 +71,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray for (k=0; kLastExterior(); n++){ -// va = cDen[n]; -// vb = cDen[Np + n]; -// value = (va-vb)/(va+vb); -// idx = TmpMap[n]; -// if (!(idx < 0) && idxFirstInterior(); nLastInterior(); n++){ -// va = cDen[n]; -// vb = cDen[Np + n]; -// value = (va-vb)/(va+vb); -// idx = TmpMap[n]; -// if (!(idx < 0) && idxBarrier(); -// comm.barrier(); -// -// if (rank==0) printf ("Initializing phase and density fields on device from Restart\n"); -// //TODO the following function is to be updated. -// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, 0, ScaLBL_Comm->LastExterior(), Np); -// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); } } double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){ - int nprocs=nprocx*nprocy*nprocz; int START_TIME = timestep; int EXIT_TIME = min(returntime, timestepMax); diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp index 308cc1e6..11f8ed63 100644 --- a/models/GreyscaleModel.cpp +++ b/models/GreyscaleModel.cpp @@ -150,7 +150,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){ // Generate the signed distance map // Initialize the domain and communication Array id_solid(Nx,Ny,Nz); - int count = 0; // Solve for the position of the solid phase for (int k=0;kid[i] = Mask->id[i]; - for (int idx=0; idxComm.sumReduce( label_count[idx]); + for (size_t idx=0; idxComm.sumReduce( label_count[idx]); //Initialize a weighted porosity after considering grey voxels GreyPorosity=0.0; for (unsigned int idx=0; idx0 @@ -691,23 +691,22 @@ void ScaLBL_GreyscaleModel::Run(){ //double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag; double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag; - if (rank==0){ + if (rank==0){ printf(" AbsPerm = %.5g [micron^2]\n",absperm); - bool WriteHeader=false; - FILE * log_file = fopen("Permeability.csv","r"); - if (log_file != NULL) - fclose(log_file); - else - WriteHeader=true; - log_file = fopen("Permeability.csv","a"); - if (WriteHeader) - fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n", - timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm); + bool WriteHeader=false; + FILE * log_file = fopen("Permeability.csv","r"); + if (log_file != NULL) + fclose(log_file); + else + WriteHeader=true; + log_file = fopen("Permeability.csv","a"); + if (WriteHeader) + fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n"); - fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, - h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm); - fclose(log_file); - } + fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, + h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm); + fclose(log_file); + } } if (timestep%visualization_interval==0){ diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 136af02f..f86bc6d2 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -97,7 +97,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector &num_iter){ } else{ time_conv.clear(); - for (int i=0; i &num_iter){ ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n"); } else{ - for (int i=0; i &num_iter){ ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n"); } else{ - for (int i=0; i &num_iter){ else { ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n"); } - for (unsigned int i=0;i &num_iter){ else { ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n"); } - for (unsigned int i=0;iid[n]; AFFINITY=0.f; // Assign the affinity from the paired list - for (unsigned int idx=0; idx < NLABELS; idx++){ + for (size_t idx=0; idx < NLABELS; idx++){ //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]); if (VALUE == LabelList[idx]){ AFFINITY=AffinityList[idx]; @@ -701,23 +703,23 @@ void ScaLBL_IonModel::Initialize(){ auto File_ion = ion_db->getVector( "IonConcentrationFile" ); double *Ci_host; Ci_host = new double[number_ion_species*Np]; - for (int ic=0; icFirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np); } @@ -734,35 +736,35 @@ void ScaLBL_IonModel::Initialize(){ break; } - for (int i=0; i rlx; - for (unsigned int ic=0;icBarrier(); comm.barrier(); //auto t1 = std::chrono::system_clock::now(); - for (int ic=0; icFirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np); } @@ -901,7 +903,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){ //if (rank==0) printf("********************************************************\n"); } -void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){ +void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic){ //This function wirte out the data in a normal layout (by aggregating all decomposed domains) ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration); @@ -913,13 +915,13 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const i void ScaLBL_IonModel::getIonConcentration_debug(int timestep){ //This function write out decomposed data DoubleArray PhaseField(Nx,Ny,Nz); - for (int ic=0; icRegularLayout(Map,&Ci[ic*Np],PhaseField); ScaLBL_Comm->Barrier(); comm.barrier(); IonConcentration_LB_to_Phys(PhaseField); FILE *OUTFILE; - sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank); + sprintf(LocalRankFilename,"Ion%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); OUTFILE = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,OUTFILE); fclose(OUTFILE); @@ -986,7 +988,7 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector &ci_avg_previous){ Ci_host = new double[Np]; vector error(number_ion_species,0.0); - for (int ic=0; ic BoundaryConditionInlet; vector BoundaryConditionOutlet; vector IonDiffusivity;//User input unit [m^2/sec] diff --git a/tests/DataAggregator.cpp b/tests/DataAggregator.cpp index 7fcbfb0f..137ebd61 100644 --- a/tests/DataAggregator.cpp +++ b/tests/DataAggregator.cpp @@ -9,8 +9,8 @@ using namespace std; int main(int argc, char **argv){ - - printf("Aggregating block data into single file \n"); + + printf("Aggregating block data into single file \n"); unsigned int Bx,By,Bz; uint64_t Nx,Ny,Nz; uint64_t N; @@ -22,25 +22,25 @@ int main(int argc, char **argv){ //Read the block size if (argc>9){ - printf("Input arguments accepted \n"); - Nx = atoi(argv[1]); - Ny = atoi(argv[2]); - Nz = atoi(argv[3]); - x0 = atoi(argv[4]); - y0 = atoi(argv[5]); - z0 = atoi(argv[6]); - NX = atol(argv[7]); - NY = atol(argv[8]); - NZ = atol(argv[9]); - printf("Size %i X %i X %i \n",NX,NY,NZ); - fflush(stdout); + printf("Input arguments accepted \n"); + Nx = atoi(argv[1]); + Ny = atoi(argv[2]); + Nz = atoi(argv[3]); + x0 = atoi(argv[4]); + y0 = atoi(argv[5]); + z0 = atoi(argv[6]); + NX = atol(argv[7]); + NY = atol(argv[8]); + NZ = atol(argv[9]); + printf("Size %llu X %llu X %llu \n",(unsigned long long) NX, (unsigned long long) NY, (unsigned long long) NZ); + fflush(stdout); } else{ - printf("setting defaults \n"); - Nx=Ny=Nz=1024; - x0=y0=z0=0; - NX=NY=8640; - NZ=6480; + printf("setting defaults \n"); + Nx=Ny=Nz=1024; + x0=y0=z0=0; + NX=NY=8640; + NZ=6480; } //Bx = By = Bz = 9; //Nx = Ny = Nz = 1024; @@ -53,29 +53,28 @@ int main(int argc, char **argv){ if (By>8) By=8; if (Bz>8) Bz=8; - printf("System size (output) is: %i x %i x %i \n",NX,NY,NZ); - printf("Block size (read) is: %i x %i x %i \n",Nx,Ny,Nz); - printf("Starting location (read) is: %i, %i, %i \n", x0,y0,z0); - printf("Block number (read): %i x %i x %i \n",Bx,By,Bz); + printf("System size (output) is: %llu x %llu x %llu \n",(unsigned long long) NX,(unsigned long long) NY, (unsigned long long) NZ); + printf("Block size (read) is: %llu x %llu x %llu \n",(unsigned long long) Nx,(unsigned long long) Ny,(unsigned long long) Nz); + printf("Starting location (read) is: %llu, %llu, %llu \n", (unsigned long long) x0,(unsigned long long) y0,(unsigned long long) z0); + printf("Block number (read): %llu x %llu x %llu \n",(unsigned long long) Bx,(unsigned long long) By,(unsigned long long) Bz); fflush(stdout); - + // Filenames used //char LocalRankString[8]; char LocalRankFilename[40]; char sx[2]; char sy[2]; char sz[2]; - char tmpstr[10]; //sprintf(LocalRankString,"%05d",rank); N = Nx*Ny*Nz; N_full=NX*NY*NZ; - + char *id; id = new char [N]; char *ID; ID = new char [N_full]; - + for (unsigned int bz=0; bz SW) { sw_diff_old = sw_diff_new; diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 51ffe96e..3dade4d5 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -9,8 +9,6 @@ #include "common/Utilities.h" #include "models/ColorModel.h" -//#define WRE_SURFACES - /* * Simulator for two-phase flow in porous media * James E. McClure 2013-2014 @@ -24,100 +22,170 @@ int main( int argc, char **argv ) { - // Initialize - Utilities::startup( argc, argv ); + // Initialize + Utilities::startup( argc, argv ); - { // Limit scope so variables that contain communicators will free before MPI_Finialize + { // Limit scope so variables that contain communicators will free before MPI_Finialize - Utilities::MPI comm( MPI_COMM_WORLD ); - int rank = comm.getRank(); - int nprocs = comm.getSize(); - std::string SimulationMode = "production"; - // Load the input database - auto db = std::make_shared( argv[1] ); - if (argc > 2) { - SimulationMode = "development"; - } + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + std::string SimulationMode = "production"; + // Load the input database + auto db = std::make_shared( argv[1] ); + if (argc > 2) { + SimulationMode = "development"; + } - if ( rank == 0 ) { - printf( "********************************************************\n" ); - printf( "Running Color LBM \n" ); - printf( "********************************************************\n" ); - if (SimulationMode == "development") - printf("**** DEVELOPMENT MODE ENABLED *************\n"); - } - // Initialize compute device - int device = ScaLBL_SetDevice( rank ); - NULL_USE( device ); - ScaLBL_DeviceBarrier(); - comm.barrier(); + if ( rank == 0 ) { + printf( "********************************************************\n" ); + printf( "Running Color LBM \n" ); + printf( "********************************************************\n" ); + if (SimulationMode == "development") + printf("**** DEVELOPMENT MODE ENABLED *************\n"); + } + // Initialize compute device + int device = ScaLBL_SetDevice( rank ); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); - PROFILE_ENABLE( 1 ); - // PROFILE_ENABLE_TRACE(); - // PROFILE_ENABLE_MEMORY(); - PROFILE_SYNCHRONIZE(); - PROFILE_START( "Main" ); - Utilities::setErrorHandlers(); + PROFILE_ENABLE( 1 ); + // PROFILE_ENABLE_TRACE(); + // PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START( "Main" ); + Utilities::setErrorHandlers(); - auto filename = argv[1]; - ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); - ColorModel.ReadParams( filename ); - ColorModel.SetDomain(); - ColorModel.ReadInput(); - ColorModel.Create(); // creating the model will create data structure to match the pore - // structure and allocate variables - ColorModel.Initialize(); // initializing the model will set initial conditions for variables - - if (SimulationMode == "development"){ - double MLUPS=0.0; - int timestep = 0; - /* flow adaptor keys to control */ - auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" ); - int MAX_STEADY_TIME = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); - int SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 100000 ); + auto filename = argv[1]; + ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); + ColorModel.ReadParams( filename ); + ColorModel.SetDomain(); + ColorModel.ReadInput(); + ColorModel.Create(); // creating the model will create data structure to match the pore + // structure and allocate variables + ColorModel.Initialize(); // initializing the model will set initial conditions for variables - int ANALYSIS_INTERVAL = ColorModel.timestepMax; - if (ColorModel.analysis_db->keyExists( "analysis_interval" )){ - ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar( "analysis_interval" ); - } - /* Launch the simulation */ - FlowAdaptor Adapt(ColorModel); - runAnalysis analysis(ColorModel); - - while (ColorModel.timestep < ColorModel.timestepMax){ - /* this will run steady points */ - timestep += MAX_STEADY_TIME; - MLUPS = ColorModel.Run(timestep); - if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); - Adapt.UpdateFractionalFlow(ColorModel); - - /* apply timestep skipping algorithm to accelerate steady-state */ - int skip_time = 0; - timestep = ColorModel.timestep; - while (skip_time < SKIP_TIMESTEPS){ - timestep += ANALYSIS_INTERVAL; - MLUPS = ColorModel.Run(timestep); - Adapt.MoveInterface(ColorModel); - skip_time += ANALYSIS_INTERVAL; - } - //Adapt.Flatten(ColorModel); + if (SimulationMode == "development"){ + double MLUPS=0.0; + int timestep = 0; + bool ContinueSimulation = true; + + /* Variables for simulation protocols */ + auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "none" ); + /* image sequence protocol */ + int IMAGE_INDEX = 0; + int IMAGE_COUNT = 0; + std::vector ImageList; + /* flow adaptor keys to control behavior */ + int SKIP_TIMESTEPS = 0; + int MAX_STEADY_TIME = 1000000; + double ENDPOINT_THRESHOLD = 0.1; + double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled + double SEED_WATER = 0.0; + if (ColorModel.db->keyExists( "FlowAdaptor" )){ + auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" ); + MAX_STEADY_TIME = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); + SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 50000 ); + ENDPOINT_THRESHOLD = flow_db->getWithDefault( "endpoint_threshold", 0.1); + /* protocol specific key values */ + if (PROTOCOL == "fractional flow") + FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault( "fractional_flow_increment", 0.05); + if (PROTOCOL == "seed water") + SEED_WATER = flow_db->getWithDefault( "seed_water", 0.01); + } + /* analysis keys*/ + int ANALYSIS_INTERVAL = ColorModel.timestepMax; + if (ColorModel.analysis_db->keyExists( "analysis_interval" )){ + ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar( "analysis_interval" ); + } + /* Launch the simulation */ + FlowAdaptor Adapt(ColorModel); + runAnalysis analysis(ColorModel); + while (ContinueSimulation){ + /* this will run steady points */ + timestep += MAX_STEADY_TIME; + MLUPS = ColorModel.Run(timestep); + if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); + if (ColorModel.timestep > ColorModel.timestepMax){ + ContinueSimulation = false; + } + + /* Load a new image if image sequence is specified */ + if (PROTOCOL == "image sequence"){ + IMAGE_INDEX++; + 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); + ColorModel.color_db->putScalar("image_index",IMAGE_INDEX); + Adapt.ImageInit(ColorModel, next_image); + } + else{ + if (rank==0) printf("Finished simulating image sequence \n"); + ColorModel.timestep = ColorModel.timestepMax; + } + } + /*********************************************************/ + /* update the fluid configuration with the flow adapter */ + int skip_time = 0; + timestep = ColorModel.timestep; + /* get the averaged flow measures computed internally for the last simulation point*/ + double SaturationChange = 0.0; + double volB = ColorModel.Averages->gwb.V; + double volA = ColorModel.Averages->gnb.V; + double initialSaturation = volB/(volA + volB); + double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; + double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; + double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; + double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; + double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; + double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; + double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); + double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); + /* stop simulation if previous point was sufficiently close to the endpoint*/ + if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false; + if (ContinueSimulation){ + while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){ + timestep += ANALYSIS_INTERVAL; + if (PROTOCOL == "fractional flow") { + Adapt.UpdateFractionalFlow(ColorModel); + } + else if (PROTOCOL == "shell aggregation"){ + double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange; + Adapt.ShellAggregation(ColorModel,target_volume_change); + } + else if (PROTOCOL == "seed water"){ + Adapt.SeedPhaseField(ColorModel,SEED_WATER); + } + /* Run some LBM timesteps to let the system relax a bit */ + MLUPS = ColorModel.Run(timestep); + /* Recompute the volume fraction now that the system has adjusted */ + double volB = ColorModel.Averages->gwb.V; + double volA = ColorModel.Averages->gnb.V; + SaturationChange = volB/(volA + volB) - initialSaturation; + skip_time += ANALYSIS_INTERVAL; + } + if (rank==0) printf(" ********************************************************************* \n"); + if (rank==0) printf(" Updated fractional flow with saturation change = %f \n", SaturationChange); + if (rank==0) printf(" Used protocol = %s \n", PROTOCOL.c_str()); + if (rank==0) printf(" ********************************************************************* \n"); + } + /*********************************************************/ + } + } + else + ColorModel.Run(); - } - ColorModel.WriteDebug(); - } - - else - ColorModel.Run(); - - PROFILE_STOP( "Main" ); - auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); - auto level = db->getWithDefault( "TimerLevel", 1 ); - PROFILE_SAVE( file, level ); - // **************************************************** + PROFILE_STOP( "Main" ); + auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); + auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); + PROFILE_SAVE( file, level ); + // **************************************************** - } // Limit scope so variables that contain communicators will free before MPI_Finialize + } // Limit scope so variables that contain communicators will free before MPI_Finialize - Utilities::shutdown(); - return 0; + Utilities::shutdown(); + return 0; } diff --git a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp index 19d99b9c..30ec4fec 100644 --- a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp +++ b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp @@ -59,6 +59,7 @@ int main( int argc, char **argv ) PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_SingleFluidBGK_simulator" ); auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); PROFILE_SAVE( file,level ); // **************************************************** diff --git a/tests/lbpm_freelee_simulator.cpp b/tests/lbpm_freelee_simulator.cpp index cdba5082..e2ca2d64 100644 --- a/tests/lbpm_freelee_simulator.cpp +++ b/tests/lbpm_freelee_simulator.cpp @@ -83,6 +83,7 @@ int main( int argc, char **argv ) PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_simulator" ); auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); PROFILE_SAVE( file,level ); // ****************************************************