added membrane analysis capability
This commit is contained in:
James McClure
parent
b6227dd823
commit
2acaa335aa
|
|
@ -67,6 +67,87 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr<Domain> dm)
|
|||
}
|
||||
}
|
||||
|
||||
ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel)
|
||||
: Dm(IonModel.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;
|
||||
|
||||
if (Dm->rank()==0) printf("Analyze system with sub-domain size = %i x %i x %i \n",Nx,Ny,Nz);
|
||||
|
||||
USE_MEMBRANE = IonModel.USE_MEMBRANE;
|
||||
|
||||
ChemicalPotential.resize(Nx, Ny, Nz);
|
||||
ChemicalPotential.fill(0);
|
||||
ElectricalPotential.resize(Nx, Ny, Nz);
|
||||
ElectricalPotential.fill(0);
|
||||
ElectricalField_x.resize(Nx, Ny, Nz);
|
||||
ElectricalField_x.fill(0);
|
||||
ElectricalField_y.resize(Nx, Ny, Nz);
|
||||
ElectricalField_y.fill(0);
|
||||
ElectricalField_z.resize(Nx, Ny, Nz);
|
||||
ElectricalField_z.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);
|
||||
IonFluxDiffusive_x.resize(Nx, Ny, Nz);
|
||||
IonFluxDiffusive_x.fill(0);
|
||||
IonFluxDiffusive_y.resize(Nx, Ny, Nz);
|
||||
IonFluxDiffusive_y.fill(0);
|
||||
IonFluxDiffusive_z.resize(Nx, Ny, Nz);
|
||||
IonFluxDiffusive_z.fill(0);
|
||||
IonFluxAdvective_x.resize(Nx, Ny, Nz);
|
||||
IonFluxAdvective_x.fill(0);
|
||||
IonFluxAdvective_y.resize(Nx, Ny, Nz);
|
||||
IonFluxAdvective_y.fill(0);
|
||||
IonFluxAdvective_z.resize(Nx, Ny, Nz);
|
||||
IonFluxAdvective_z.fill(0);
|
||||
IonFluxElectrical_x.resize(Nx, Ny, Nz);
|
||||
IonFluxElectrical_x.fill(0);
|
||||
IonFluxElectrical_y.resize(Nx, Ny, Nz);
|
||||
IonFluxElectrical_y.fill(0);
|
||||
IonFluxElectrical_z.resize(Nx, Ny, Nz);
|
||||
IonFluxElectrical_z.fill(0);
|
||||
|
||||
|
||||
if (Dm->rank() == 0) {
|
||||
printf("Set up analysis routines for %i ions \n",IonModel.number_ion_species);
|
||||
|
||||
bool WriteHeader = false;
|
||||
TIMELOG = fopen("electrokinetic.csv", "r");
|
||||
if (TIMELOG != NULL)
|
||||
fclose(TIMELOG);
|
||||
else
|
||||
WriteHeader = true;
|
||||
|
||||
TIMELOG = fopen("electrokinetic.csv", "a+");
|
||||
if (WriteHeader) {
|
||||
// If timelog is empty, write a short header to list the averages
|
||||
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
|
||||
fprintf(TIMELOG, "timestep voltage_out voltage_in ");
|
||||
fprintf(TIMELOG, "voltage_out_membrane voltage_in_membrane ");
|
||||
for (size_t i=0; i<IonModel.number_ion_species; i++){
|
||||
fprintf(TIMELOG, "rho_%lu_out rho_%lu_in ",i, i);
|
||||
fprintf(TIMELOG, "rho_%lu_out_membrane rho_%lu_in_membrane ", i, i);
|
||||
}
|
||||
fprintf(TIMELOG, "count_out count_in ");
|
||||
fprintf(TIMELOG, "count_out_membrane count_in_membrane\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
ElectroChemistryAnalyzer::~ElectroChemistryAnalyzer() {
|
||||
if (Dm->rank() == 0) {
|
||||
fclose(TIMELOG);
|
||||
|
|
@ -75,6 +156,163 @@ ElectroChemistryAnalyzer::~ElectroChemistryAnalyzer() {
|
|||
|
||||
void ElectroChemistryAnalyzer::SetParams() {}
|
||||
|
||||
void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion,
|
||||
ScaLBL_Poisson &Poisson,
|
||||
int timestep) {
|
||||
|
||||
int i, j, k;
|
||||
|
||||
Poisson.getElectricPotential(ElectricalPotential);
|
||||
|
||||
if (Dm->rank() == 0)
|
||||
fprintf(TIMELOG, "%i ", timestep);
|
||||
|
||||
|
||||
/* int iq, ip, nq, np, nqm, npm;
|
||||
Ion.MembraneDistance(i,j,k); // inside (-) or outside (+) the ion
|
||||
for (int link; link<Ion.IonMembrane->membraneLinkCount; link++){
|
||||
int iq = Ion.IonMembrane->membraneLinks[2*link];
|
||||
int ip = Ion.IonMembrane->membraneLinks[2*link+1];
|
||||
|
||||
iq = membrane[2*link]; ip = membrane[2*link+1];
|
||||
nq = iq%Np; np = ip%Np;
|
||||
nqm = Map[nq]; npm = Map[np];
|
||||
}
|
||||
*/
|
||||
unsigned long int in_local_count, out_local_count;
|
||||
unsigned long int in_global_count, out_global_count;
|
||||
|
||||
double value_in_local, value_out_local;
|
||||
double value_in_global, value_out_global;
|
||||
|
||||
double value_membrane_in_local, value_membrane_out_local;
|
||||
double value_membrane_in_global, value_membrane_out_global;
|
||||
|
||||
unsigned long int membrane_in_local_count, membrane_out_local_count;
|
||||
unsigned long int membrane_in_global_count, membrane_out_global_count;
|
||||
|
||||
double memdist,value;
|
||||
in_local_count = 0;
|
||||
out_local_count = 0;
|
||||
membrane_in_local_count = 0;
|
||||
membrane_out_local_count = 0;
|
||||
|
||||
value_membrane_in_local = 0.0;
|
||||
value_membrane_out_local = 0.0;
|
||||
value_in_local = 0.0;
|
||||
value_out_local = 0.0;
|
||||
for (k = 1; k < Nz; k++) {
|
||||
for (j = 1; j < Ny; j++) {
|
||||
for (i = 1; i < Nx; i++) {
|
||||
/* electric potential */
|
||||
memdist = Ion.MembraneDistance(i,j,k);
|
||||
value = ElectricalPotential(i,j,k);
|
||||
if (memdist < 0.0){
|
||||
// inside the membrane
|
||||
if (fabs(memdist) < 1.0){
|
||||
value_membrane_in_local += value;
|
||||
membrane_in_local_count++;
|
||||
}
|
||||
value_in_local += value;
|
||||
in_local_count++;
|
||||
|
||||
}
|
||||
else {
|
||||
// outside the membrane
|
||||
if (fabs(memdist) < 1.0){
|
||||
value_membrane_out_local += value;
|
||||
membrane_out_local_count++;
|
||||
}
|
||||
value_out_local += value;
|
||||
out_local_count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/* these only need to be computed the first time through */
|
||||
out_global_count = Dm->Comm.sumReduce(out_local_count);
|
||||
in_global_count = Dm->Comm.sumReduce(in_local_count);
|
||||
membrane_out_global_count = Dm->Comm.sumReduce(membrane_out_local_count);
|
||||
membrane_in_global_count = Dm->Comm.sumReduce(membrane_in_local_count);
|
||||
|
||||
value_out_global = Dm->Comm.sumReduce(value_out_local);
|
||||
value_in_global = Dm->Comm.sumReduce(value_in_local);
|
||||
value_membrane_out_global = Dm->Comm.sumReduce(value_membrane_out_local);
|
||||
value_membrane_in_global = Dm->Comm.sumReduce(value_membrane_in_local);
|
||||
|
||||
value_out_global /= out_global_count;
|
||||
value_in_global /= in_global_count;
|
||||
value_membrane_out_global /= membrane_out_global_count;
|
||||
value_membrane_in_global /= membrane_in_global_count;
|
||||
|
||||
if (Dm->rank() == 0) {
|
||||
fprintf(TIMELOG, "%.8g ", value_out_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_in_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_membrane_out_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_membrane_in_global);
|
||||
}
|
||||
|
||||
value_membrane_in_local = 0.0;
|
||||
value_membrane_out_local = 0.0;
|
||||
value_in_local = 0.0;
|
||||
value_out_local = 0.0;
|
||||
for (size_t ion = 0; ion < Ion.number_ion_species; ion++) {
|
||||
Ion.getIonConcentration(Rho, ion);
|
||||
value_membrane_in_local = 0.0;
|
||||
value_membrane_out_local = 0.0;
|
||||
value_in_local = 0.0;
|
||||
value_out_local = 0.0;
|
||||
for (k = 1; k < Nz; k++) {
|
||||
for (j = 1; j < Ny; j++) {
|
||||
for (i = 1; i < Nx; i++) {
|
||||
/* electric potential */
|
||||
memdist = Ion.MembraneDistance(i,j,k);
|
||||
value = Rho(i,j,k);
|
||||
if (memdist < 0.0){
|
||||
// inside the membrane
|
||||
if (fabs(memdist) < 1.0){
|
||||
value_membrane_in_local += value;
|
||||
}
|
||||
value_in_local += value;
|
||||
}
|
||||
else {
|
||||
// outside the membrane
|
||||
if (fabs(memdist) < 1.0){
|
||||
value_membrane_out_local += value;
|
||||
}
|
||||
value_out_local += value;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
value_out_global = Dm->Comm.sumReduce(value_out_local);
|
||||
value_in_global = Dm->Comm.sumReduce(value_in_local);
|
||||
value_membrane_out_global = Dm->Comm.sumReduce(value_membrane_out_local);
|
||||
value_membrane_in_global = Dm->Comm.sumReduce(value_membrane_in_local);
|
||||
|
||||
value_out_global /= out_global_count;
|
||||
value_in_global /= in_global_count;
|
||||
value_membrane_out_global /= membrane_out_global_count;
|
||||
value_membrane_in_global /= membrane_in_global_count;
|
||||
|
||||
if (Dm->rank() == 0) {
|
||||
fprintf(TIMELOG, "%.8g ", value_out_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_in_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_membrane_out_global);
|
||||
fprintf(TIMELOG, "%.8g ", value_membrane_in_global);
|
||||
}
|
||||
}
|
||||
|
||||
if (Dm->rank() == 0) {
|
||||
fprintf(TIMELOG, "%lu ", out_global_count);
|
||||
fprintf(TIMELOG, "%lu ", in_global_count);
|
||||
fprintf(TIMELOG, "%lu ", membrane_out_global_count);
|
||||
fprintf(TIMELOG, "%lu\n", membrane_in_global_count);
|
||||
fflush(TIMELOG);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion,
|
||||
ScaLBL_Poisson &Poisson,
|
||||
ScaLBL_StokesModel &Stokes, int timestep) {
|
||||
|
|
|
|||
|
|
@ -56,14 +56,14 @@ public:
|
|||
DoubleArray IonFluxElectrical_z;
|
||||
|
||||
ElectroChemistryAnalyzer(std::shared_ptr<Domain> Dm);
|
||||
ElectroChemistryAnalyzer( ScaLBL_IonModel &IonModel);
|
||||
~ElectroChemistryAnalyzer();
|
||||
|
||||
void SetParams();
|
||||
void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
|
||||
ScaLBL_StokesModel &Stokes, int timestep);
|
||||
void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
|
||||
ScaLBL_StokesModel &Stokes,
|
||||
std::shared_ptr<Database> input_db, int timestep);
|
||||
void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, int timestep);
|
||||
void Membrane(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, int timestep);
|
||||
void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
|
||||
ScaLBL_StokesModel &Stokes,std::shared_ptr<Database> input_db, int timestep);
|
||||
void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, int timestep);
|
||||
void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
|
||||
std::shared_ptr<Database> input_db, int timestep);
|
||||
|
|
|
|||
|
|
@ -58,7 +58,7 @@ public:
|
|||
int *NeighborList; // modified neighborlist
|
||||
|
||||
/* host data structures */
|
||||
int *membraneLinks; // D3Q19 links that cross membrane
|
||||
int *membraneLinks; // D3Q7 links that cross membrane
|
||||
int *membraneTag; // label each link in the membrane
|
||||
double *membraneDist; // distance to membrane for each linked site
|
||||
double *membraneOrientation; // distance to membrane for each linked site
|
||||
|
|
|
|||
|
|
@ -22,148 +22,148 @@ using namespace std;
|
|||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// Initialize MPI and error handlers
|
||||
//Utilities::startup( argc, argv );
|
||||
Utilities::startup( argc, argv, true );
|
||||
// Initialize MPI and error handlers
|
||||
//Utilities::startup( argc, argv );
|
||||
Utilities::startup( argc, argv, true );
|
||||
|
||||
Utilities::MPI comm( MPI_COMM_WORLD );
|
||||
int rank = comm.getRank();
|
||||
int nprocs = comm.getSize();
|
||||
Utilities::MPI comm( MPI_COMM_WORLD );
|
||||
int rank = comm.getRank();
|
||||
int nprocs = comm.getSize();
|
||||
|
||||
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
|
||||
if (rank == 0){
|
||||
printf("********************************************************\n");
|
||||
printf("Running LBPM Nernst-Planck Membrane solver \n");
|
||||
printf("********************************************************\n");
|
||||
}
|
||||
// Initialize compute device
|
||||
int device=ScaLBL_SetDevice(rank);
|
||||
NULL_USE( device );
|
||||
ScaLBL_DeviceBarrier();
|
||||
comm.barrier();
|
||||
if (rank == 0){
|
||||
printf("********************************************************\n");
|
||||
printf("Running LBPM Nernst-Planck Membrane solver \n");
|
||||
printf("********************************************************\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_StokesModel StokesModel(rank,nprocs,comm);
|
||||
ScaLBL_IonModel IonModel(rank,nprocs,comm);
|
||||
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
|
||||
ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
|
||||
auto filename = argv[1];
|
||||
//ScaLBL_StokesModel StokesModel(rank,nprocs,comm);
|
||||
ScaLBL_IonModel IonModel(rank,nprocs,comm);
|
||||
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
|
||||
ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
|
||||
|
||||
bool SlipBC = false;
|
||||
|
||||
// Load controller information
|
||||
Study.ReadParams(filename);
|
||||
bool SlipBC = false;
|
||||
|
||||
// Load user input database files for Navier-Stokes and Ion solvers
|
||||
//StokesModel.ReadParams(filename);
|
||||
// Load controller information
|
||||
Study.ReadParams(filename);
|
||||
|
||||
// Setup other model specific structures
|
||||
//StokesModel.SetDomain();
|
||||
//StokesModel.ReadInput();
|
||||
//StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf("Stokes model setup complete\n");
|
||||
|
||||
IonModel.ReadParams(filename);
|
||||
IonModel.SetDomain();
|
||||
IonModel.ReadInput();
|
||||
IonModel.Create();
|
||||
IonModel.SetMembrane();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Ion model setup complete\n");
|
||||
fflush(stdout);
|
||||
|
||||
// Create analysis object
|
||||
ElectroChemistryAnalyzer Analysis(IonModel.Dm);
|
||||
// Load user input database files for Navier-Stokes and Ion solvers
|
||||
//StokesModel.ReadParams(filename);
|
||||
|
||||
// Get internal iteration number
|
||||
//StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
//StokesModel.Initialize(); // initializing the model will set initial conditions for variables
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf("Stokes model initialized \n");
|
||||
//fflush(stdout);
|
||||
|
||||
//IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
IonModel.timestepMax = Study.getIonNumIter_NernstPlanck_coupling(IonModel.time_conv);
|
||||
IonModel.Initialize();
|
||||
IonModel.DummyFluidVelocity();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Ion model initialized \n");
|
||||
// Get maximal time converting factor based on Sotkes and Ion solvers
|
||||
//Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
Study.time_conv_MainLoop = IonModel.timestepMax[0]*IonModel.time_conv[0];
|
||||
// Setup other model specific structures
|
||||
//StokesModel.SetDomain();
|
||||
//StokesModel.ReadInput();
|
||||
//StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf("Stokes model setup complete\n");
|
||||
|
||||
//----------------------------------- print out for debugging ------------------------------------------//
|
||||
if (rank==0){
|
||||
for (size_t i=0;i<IonModel.timestepMax.size();i++){
|
||||
printf("Main loop time_conv computed from ion %i: %.5g[s/lt]\n",i+1,IonModel.timestepMax[i]*IonModel.time_conv[i]);
|
||||
}
|
||||
}
|
||||
//----------------------------------- print out for debugging ------------------------------------------//
|
||||
IonModel.ReadParams(filename);
|
||||
IonModel.SetDomain();
|
||||
IonModel.ReadInput();
|
||||
IonModel.Create();
|
||||
IonModel.SetMembrane();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Ion model setup complete\n");
|
||||
fflush(stdout);
|
||||
|
||||
// Initialize LB-Poisson model
|
||||
PoissonSolver.ReadParams(filename);
|
||||
PoissonSolver.SetDomain();
|
||||
PoissonSolver.ReadInput();
|
||||
PoissonSolver.Create();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Poisson solver created \n");
|
||||
fflush(stdout);
|
||||
PoissonSolver.Initialize(Study.time_conv_MainLoop);
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Poisson solver initialized \n");
|
||||
fflush(stdout);
|
||||
// Create analysis object
|
||||
ElectroChemistryAnalyzer Analysis(IonModel);
|
||||
|
||||
int timestep=0;
|
||||
while (timestep < Study.timestepMax){
|
||||
|
||||
timestep++;
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf(" Poisson step %i \n",timestep);
|
||||
//StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
|
||||
//fflush(stdout);
|
||||
// Get internal iteration number
|
||||
//StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
//StokesModel.Initialize(); // initializing the model will set initial conditions for variables
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf("Stokes model initialized \n");
|
||||
//fflush(stdout);
|
||||
|
||||
IonModel.RunMembrane(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane
|
||||
//comm.barrier();
|
||||
if (rank == 0) printf(" Membrane step %i \n",timestep);
|
||||
fflush(stdout);
|
||||
//IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
IonModel.timestepMax = Study.getIonNumIter_NernstPlanck_coupling(IonModel.time_conv);
|
||||
IonModel.Initialize();
|
||||
IonModel.DummyFluidVelocity();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Ion model initialized \n");
|
||||
// Get maximal time converting factor based on Sotkes and Ion solvers
|
||||
//Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
|
||||
Study.time_conv_MainLoop = IonModel.timestepMax[0]*IonModel.time_conv[0];
|
||||
|
||||
//----------------------------------- print out for debugging ------------------------------------------//
|
||||
if (rank==0){
|
||||
for (size_t i=0;i<IonModel.timestepMax.size();i++){
|
||||
printf("Main loop time_conv computed from ion %i: %.5g[s/lt]\n",i+1,IonModel.timestepMax[i]*IonModel.time_conv[i]);
|
||||
}
|
||||
}
|
||||
//----------------------------------- print out for debugging ------------------------------------------//
|
||||
|
||||
// Initialize LB-Poisson model
|
||||
PoissonSolver.ReadParams(filename);
|
||||
PoissonSolver.SetDomain();
|
||||
PoissonSolver.ReadInput();
|
||||
PoissonSolver.Create();
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Poisson solver created \n");
|
||||
fflush(stdout);
|
||||
PoissonSolver.Initialize(Study.time_conv_MainLoop);
|
||||
comm.barrier();
|
||||
if (rank == 0) printf("Poisson solver initialized \n");
|
||||
fflush(stdout);
|
||||
|
||||
int timestep=0;
|
||||
while (timestep < Study.timestepMax){
|
||||
|
||||
timestep++;
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
|
||||
//comm.barrier();
|
||||
//if (rank == 0) printf(" Poisson step %i \n",timestep);
|
||||
//StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
|
||||
//fflush(stdout);
|
||||
|
||||
IonModel.RunMembrane(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane
|
||||
//comm.barrier();
|
||||
if (rank == 0) printf(" Membrane step %i \n",timestep);
|
||||
fflush(stdout);
|
||||
|
||||
|
||||
if (timestep%Study.analysis_interval==0){
|
||||
Analysis.Basic(IonModel,PoissonSolver,timestep);
|
||||
}
|
||||
if (timestep%Study.visualization_interval==0){
|
||||
if (timestep%Study.analysis_interval==0){
|
||||
Analysis.Membrane(IonModel,PoissonSolver,timestep);
|
||||
}
|
||||
if (timestep%Study.visualization_interval==0){
|
||||
|
||||
Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep);
|
||||
//StokesModel.getVelocity(timestep);
|
||||
//PoissonSolver.getElectricPotential_debug(timestep);
|
||||
// PoissonSolver.getElectricField_debug(timestep);
|
||||
//IonModel.getIonConcentration_debug(timestep);
|
||||
|
||||
}
|
||||
}
|
||||
Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep);
|
||||
//StokesModel.getVelocity(timestep);
|
||||
//PoissonSolver.getElectricPotential_debug(timestep);
|
||||
// PoissonSolver.getElectricField_debug(timestep);
|
||||
//IonModel.getIonConcentration_debug(timestep);
|
||||
|
||||
if (rank==0) printf("Save simulation raw data at maximum timestep\n");
|
||||
Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep);
|
||||
}
|
||||
}
|
||||
|
||||
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
|
||||
if (rank==0) printf("*************************************************************\n");
|
||||
if (rank==0) printf("Save simulation raw data at maximum timestep\n");
|
||||
Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep);
|
||||
|
||||
PROFILE_STOP("Main");
|
||||
PROFILE_SAVE("lbpm_nernst_planck_membrane_simulator",1);
|
||||
// ****************************************************
|
||||
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
|
||||
if (rank==0) printf("*************************************************************\n");
|
||||
|
||||
} // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
PROFILE_STOP("Main");
|
||||
PROFILE_SAVE("lbpm_nernst_planck_membrane_simulator",1);
|
||||
// ****************************************************
|
||||
|
||||
Utilities::shutdown();
|
||||
} // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
|
||||
Utilities::shutdown();
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user