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

added visualization capability for Lee model

Browse Source
This commit is contained in:
James McClure 2021-03-21 00:15:23 -04:00
parent 9f9b0dbffe
commit d913b9bdc8
6 changed files with 278 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

181
analysis/FreeEnergy.cpp Normal file
View File

@ -0,0 +1,181 @@
#include "analysis/FreeEnergy.h"
FreeEnergyAnalyzer::FreeEnergyAnalyzer(std::shared_ptr <Domain> dm):
Dm(dm)
{
Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0;
ChemicalPotential.resize(Nx,Ny,Nz); ChemicalPotential.fill(0);
Phi.resize(Nx,Ny,Nz); Phi.fill(0);
Pressure.resize(Nx,Ny,Nz); Pressure.fill(0);
Rho.resize(Nx,Ny,Nz); Rho.fill(0);
Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field
Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0);
Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0);
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
if (Dm->rank()==0){
bool WriteHeader=false;
TIMELOG = fopen("free.csv","r");
if (TIMELOG != NULL)
fclose(TIMELOG);
else
WriteHeader=true;
TIMELOG = fopen("free.csv","a+");
if (WriteHeader)
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"timestep\n");
}
}
}
FreeEnergyAnalyzer::~FreeEnergyAnalyzer(){
if (Dm->rank()==0){
fclose(TIMELOG);
}
}
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<Ion.number_ion_species; ion++){
fprintf(TIMELOG,"%.8g ",rho_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_fluctuation_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_fluctuation_global[ion]);
}
*/
fprintf(TIMELOG,"\n");
fflush(TIMELOG);
}
/* else{
fprintf(TIMELOG,"%i ",timestep);
for (int ion=0; ion<Ion.number_ion_species; ion++){
fprintf(TIMELOG,"%.8g ",rho_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_fluctuation_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_fluctuation_local[ion]);
}
fflush(TIMELOG);
} */
}
void FreeEnergyAnalyzer::WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep){
auto vis_db = input_db->getDatabase( "Visualization" );
char VisName[40];
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);
IO::initialize("","silo","false");
// Create the MeshDataStruct
visData.resize(1);
visData[0].meshName = "domain";
visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto VisPhase = std::make_shared<IO::Variable>();
auto VisPressure = std::make_shared<IO::Variable>();
auto VisChemicalPotential = std::make_shared<IO::Variable>();
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
if (vis_db->getWithDefault<bool>( "save_phase_field", true )){
VisPhase->name = "Phase";
VisPhase->type = IO::VariableType::VolumeVariable;
VisPhase->dim = 1;
VisPhase->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VisPhase);
}
if (vis_db->getWithDefault<bool>( "save_potential", true )){
VisPressure->name = "Pressure";
VisPressure->type = IO::VariableType::VolumeVariable;
VisPressure->dim = 1;
VisPressure->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VisPressure);
VisChemicalPotential->name = "ChemicalPotential";
VisChemicalPotential->type = IO::VariableType::VolumeVariable;
VisChemicalPotential->dim = 1;
VisChemicalPotential->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VisChemicalPotential);
}
if (vis_db->getWithDefault<bool>( "save_velocity", false )){
VxVar->name = "Velocity_x";
VxVar->type = IO::VariableType::VolumeVariable;
VxVar->dim = 1;
VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VxVar);
VyVar->name = "Velocity_y";
VyVar->type = IO::VariableType::VolumeVariable;
VyVar->dim = 1;
VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VyVar);
VzVar->name = "Velocity_z";
VzVar->type = IO::VariableType::VolumeVariable;
VzVar->dim = 1;
VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VzVar);
}
if (vis_db->getWithDefault<bool>( "save_phase", true )){
ASSERT(visData[0].vars[0]->name=="Phase");
LeeModel.getPhase(Phi);
Array<double>& PhaseData = visData[0].vars[0]->data;
fillData.copy(Phi,PhaseData);
}
if (vis_db->getWithDefault<bool>( "save_potential", true )){
ASSERT(visData[0].vars[1]->name=="Pressure");
LeeModel.getPotential(Pressure, ChemicalPotential);
Array<double>& PressureData = visData[0].vars[1]->data;
fillData.copy(Pressure,PressureData);
ASSERT(visData[0].vars[2]->name=="ChemicalPotential");
Array<double>& ChemicalPotentialData = visData[0].vars[2]->data;
fillData.copy(ChemicalPotential,ChemicalPotentialData);
}
if (vis_db->getWithDefault<bool>( "save_velocity", false )){
ASSERT(visData[0].vars[3]->name=="Velocity_x");
ASSERT(visData[0].vars[4]->name=="Velocity_y");
ASSERT(visData[0].vars[5]->name=="Velocity_z");
LeeModel.getVelocity(Vel_x,Vel_y,Vel_z);
Array<double>& VelxData = visData[0].vars[3]->data;
Array<double>& VelyData = visData[0].vars[4]->data;
Array<double>& VelzData = visData[0].vars[5]->data;
fillData.copy(Vel_x,VelxData);
fillData.copy(Vel_y,VelyData);
fillData.copy(Vel_z,VelzData);
}
if (vis_db->getWithDefault<bool>( "write_silo", true ))
IO::writeData( timestep, visData, Dm->Comm );
/* if (vis_db->getWithDefault<bool>( "save_8bit_raw", true )){
char CurrentIDFilename[40];
sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
Averages.AggregateLabels(CurrentIDFilename);
}
*/
}

54
analysis/FreeEnergy.h Normal file
View File

@ -0,0 +1,54 @@
/*
* averaging tools for electrochemistry
*/
#ifndef FreeEnergyAnalyzer_INC
#define FreeEnergyAnalyzer_INC
#include <vector>
#include "common/Domain.h"
#include "common/Utilities.h"
#include "common/MPI.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "analysis/distance.h"
#include "analysis/Minkowski.h"
#include "analysis/SubPhase.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#include "models/FreeLeeModel.h"
class FreeEnergyAnalyzer{
public:
std::shared_ptr <Domain> Dm;
double Volume;
// input variables
double rho_n, rho_w;
double nu_n, nu_w;
double gamma_wn, beta;
double Fx, Fy, Fz;
//...........................................................................
int Nx,Ny,Nz;
DoubleArray Rho;
DoubleArray Phi;
DoubleArray ChemicalPotential;
DoubleArray Pressure;
DoubleArray Vel_x;
DoubleArray Vel_y;
DoubleArray Vel_z;
DoubleArray SDs;
FreeEnergyAnalyzer(std::shared_ptr <Domain> Dm);
~FreeEnergyAnalyzer();
void SetParams();
void Basic( ScaLBL_FreeLeeModel &LeeModel, int timestep);
void WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep);
private:
FILE *TIMELOG;
};
#endif

28
models/FreeLeeModel.cpp
View File

@ -31,7 +31,7 @@ void ScaLBL_FreeLeeModel::getPhase(DoubleArray &PhaseValues){
for (int k=1; k<Nzh-1; k++){
for (int j=1; j<Nyh-1; j++){
for (int i=1; i<Nxh-1; i++){
PhaseValues(i-1,j-1,k-1) = PhaseData(i,j,k);
PhaseValues(i-1,j-1,k-1) = PhaseWideHalo(i,j,k);
}
}
}
@ -776,21 +776,17 @@ void ScaLBL_FreeLeeModel::Initialize_SingleFluid(){
}
}
void ScaLBL_FreeLeeModel::Run_TwoFluid(){
double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
if (rank==0){
printf("********************************************************\n");
printf("No. of timesteps: %i \n", timestepMax);
fflush(stdout);
}
int START_TIME = timestep;
int EXIT_TIME = min(returntime, timestepMax);
//************ MAIN ITERATION LOOP ***************************************/
comm.barrier();
auto t1 = std::chrono::system_clock::now();
PROFILE_START("Loop");
while (timestep < timestepMax ) {
while (timestep < EXIT_TIME ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
// *************ODD TIMESTEP*************
@ -890,19 +886,11 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
auto t2 = std::chrono::system_clock::now();
double cputime = std::chrono::duration<double>( t2 - t1 ).count() / timestep;
double cputime = std::chrono::duration<double>( t2 - t1 ).count() / (EXIT_TIME-START_TIME);
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
WriteDebug_TwoFluid();
// ************************************************************************
return MLUPS;
}
void ScaLBL_FreeLeeModel::Run_SingleFluid(){

9
models/FreeLeeModel.h
View File

@ -16,6 +16,9 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model
#include "common/ScaLBL.h"
#include "common/WideHalo.h"
#ifndef ScaLBL_FreeLeeModel_INC
#define ScaLBL_FreeLeeModel_INC
class ScaLBL_FreeLeeModel{
public:
ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM);
@ -28,11 +31,13 @@ public:
void ReadInput();
void Create_TwoFluid();
void Initialize_TwoFluid();
void Run_TwoFluid();
double Run_TwoFluid(int returntime);
void WriteDebug_TwoFluid();
void Create_SingleFluid();
void Initialize_SingleFluid();
void Run_SingleFluid();
void WriteDebug_SingleFluid();
// test utilities
void Create_DummyPhase_MGTest();
@ -97,4 +102,4 @@ private:
void AssignComponentLabels_ChemPotential_ColorGrad();
};
#endif

30
tests/lbpm_freelee_simulator.cpp
View File

@ -8,6 +8,7 @@
#include "common/Utilities.h"
#include "models/FreeLeeModel.h"
#include "analysis/FreeEnergy.h"
//*******************************************************************
// Implementation of Free-Energy Two-Phase LBM (Lee model)
@ -52,10 +53,33 @@ int main( int argc, char **argv )
LeeModel.SetDomain();
LeeModel.ReadInput();
LeeModel.Create_TwoFluid();
FreeEnergyAnalyzer Analysis(LeeModel.Dm);
LeeModel.Initialize_TwoFluid();
LeeModel.Run_TwoFluid();
LeeModel.WriteDebug_TwoFluid();
/*** RUN MAIN TIMESTEPS HERE ************/
double MLUPS=0.0;
int timestep = 0;
int visualization_time = LeeModel.timestepMax;
if (LeeModel.vis_db->keyExists( "visualizataion_interval" )){
visualization_time = LeeModel.vis_db->getScalar<int>( "visualizataion_interval" );
timestep += visualization_time;
}
while (LeeModel.timestep < LeeModel.timestepMax){
MLUPS = LeeModel.Run_TwoFluid(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
Analysis.WriteVis(LeeModel,LeeModel.db, timestep);
timestep += visualization_time;
}
//LeeModel.WriteDebug_TwoFluid();
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
// ************************************************************************
PROFILE_STOP("Main");
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );

2
tests/testGlobalMassFreeLee.cpp
View File

@ -63,7 +63,7 @@ int main( int argc, char **argv )
DoubleArray DensityInit(Nx,Ny,Nz);
LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityInit);
LeeModel.Run_TwoFluid();
double MLUPS = LeeModel.Run_TwoFluid(LeeModel.timestepMax);
DoubleArray DensityFinal(Nx,Ny,Nz);
LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityFinal);