From 6c262b53e5104df60b42bbfd2595123552d64fe9 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 23 Sep 2021 21:43:07 +1000 Subject: [PATCH 1/5] initialize the work of saving various ion fluxes --- analysis/ElectroChemistry.cpp | 12 +++++++++++- analysis/ElectroChemistry.h | 9 +++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index a762fd82..1a9a986a 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -16,7 +16,16 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm): 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){ bool WriteHeader=false; TIMELOG = fopen("electrokinetic.csv","r"); @@ -167,6 +176,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P auto VxVar = std::make_shared(); auto VyVar = std::make_shared(); auto VzVar = std::make_shared(); + if (vis_db->getWithDefault( "save_electric_potential", true )){ ElectricPotential->name = "ElectricPotential"; diff --git a/analysis/ElectroChemistry.h b/analysis/ElectroChemistry.h index 90874ca0..ee529725 100644 --- a/analysis/ElectroChemistry.h +++ b/analysis/ElectroChemistry.h @@ -40,6 +40,15 @@ public: DoubleArray Vel_y; DoubleArray Vel_z; DoubleArray SDs; + DoubleArray IonFluxDiffusive_x; //ion diffusive flux components + DoubleArray IonFluxDiffusive_y; + DoubleArray IonFluxDiffusive_z; + DoubleArray IonFluxAdvective_x; //ion advective flux components + DoubleArray IonFluxAdvective_y; + DoubleArray IonFluxAdvective_z; + DoubleArray IonFluxElectrical_x; //ion electromigration flux components + DoubleArray IonFluxElectrical_y; + DoubleArray IonFluxElectrical_z; ElectroChemistryAnalyzer(std::shared_ptr Dm); ~ElectroChemistryAnalyzer(); From 01499e672d298ddbfad3d20c8199511145e13f37 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Wed, 29 Sep 2021 15:41:38 +1000 Subject: [PATCH 2/5] continue the work of storing ion diffusive flux;incomplete. --- analysis/ElectroChemistry.cpp | 68 +++++++++++++++++++++++++++-- common/ScaLBL.h | 4 +- cpu/Ion.cpp | 22 +++++++++- cpu/Poisson.cpp | 4 +- models/IonModel.cpp | 80 ++++++++++++++++++++++++++++++++--- models/IonModel.h | 1 + tests/TestIonModel.cpp | 1 + 7 files changed, 166 insertions(+), 14 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 1a9a986a..b46a7473 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -176,8 +176,16 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P auto VxVar = std::make_shared(); auto VyVar = std::make_shared(); auto VzVar = std::make_shared(); + std::vector> IonFluxDiffusive; + for (size_t ion=0; ion()); + IonFluxDiffusive.push_back(std::make_shared()); + IonFluxDiffusive.push_back(std::make_shared()); + } + //-------------------------------------------------------------------------------------------------------------------- - + //-------------------------------------Create Names for Variables------------------------------------------------------ if (vis_db->getWithDefault( "save_electric_potential", true )){ ElectricPotential->name = "ElectricPotential"; ElectricPotential->type = IO::VariableType::VolumeVariable; @@ -214,7 +222,35 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); visData[0].vars.push_back(VzVar); } + + if (vis_db->getWithDefault( "save_ion_flux_diffusive", false )){ + for (size_t ion=0; ionname = VisName; + IonFluxDiffusive[3*ion+0]->type = IO::VariableType::VolumeVariable; + IonFluxDiffusive[3*ion+0]->dim = 1; + IonFluxDiffusive[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(IonFluxDiffusive[3*ion+0]); + // y-component of diffusive flux + sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1); + IonFluxDiffusive[3*ion+1]->name = VisName; + IonFluxDiffusive[3*ion+1]->type = IO::VariableType::VolumeVariable; + IonFluxDiffusive[3*ion+1]->dim = 1; + IonFluxDiffusive[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(IonFluxDiffusive[3*ion+1]); + // z-component of diffusive flux + sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1); + IonFluxDiffusive[3*ion+2]->name = VisName; + IonFluxDiffusive[3*ion+2]->type = IO::VariableType::VolumeVariable; + IonFluxDiffusive[3*ion+2]->dim = 1; + IonFluxDiffusive[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]); + } + } + //-------------------------------------------------------------------------------------------------------------------- + //------------------------------------Save All Variables-------------------------------------------------------------- if (vis_db->getWithDefault( "save_electric_potential", true )){ ASSERT(visData[0].vars[0]->name=="ElectricPotential"); Poisson.getElectricPotential(ElectricalPotential); @@ -225,7 +261,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P if (vis_db->getWithDefault( "save_concentration", true )){ for (size_t ion=0; ionname = VisName; + //IonConcentration[ion]->name = VisName; ASSERT(visData[0].vars[1+ion]->name==VisName); Array& IonConcentrationData = visData[0].vars[1+ion]->data; Ion.getIonConcentration(Rho,ion); @@ -245,10 +281,36 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P fillData.copy(Vel_y,VelyData); fillData.copy(Vel_z,VelzData); } + + if (vis_db->getWithDefault( "save_ion_flux_diffusive", false )){ + for (size_t ion=0; ionname = VisName; + ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->name==VisName); + // y-component of diffusive flux + sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1); + //IonFluxDiffusive[3*ion+1]->name = VisName; + ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->name==VisName); + // z-component of diffusive flux + sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1); + //IonFluxDiffusive[3*ion+2]->name = VisName; + ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->name==VisName); + + Array& IonFluxData_x = visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->data; + Array& IonFluxData_y = visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->data; + Array& IonFluxData_z = visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->data; + Ion.getIonFluxDiffusive(IonFluxDiffusive_x,IonFluxDiffusive_y,IonFluxDiffusive_z,ion); + fillData.copy(IonFluxDiffusive_x,IonFluxData_x); + fillData.copy(IonFluxDiffusive_y,IonFluxData_y); + fillData.copy(IonFluxDiffusive_z,IonFluxData_z); + } + } if (vis_db->getWithDefault( "write_silo", true )) IO::writeData( timestep, visData, Dm->Comm ); - + //-------------------------------------------------------------------------------------------------------------------- /* if (vis_db->getWithDefault( "save_8bit_raw", true )){ char CurrentIDFilename[40]; sprintf(CurrentIDFilename,"id_t%d.raw",timestep); diff --git a/common/ScaLBL.h b/common/ScaLBL.h index d37723c2..6e15b9f1 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -253,10 +253,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *di extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np); extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np); diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index 5d7ef87a..672ad71c 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -80,13 +80,14 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i } } -extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; int nr1,nr2,nr3,nr4,nr5,nr6; @@ -124,6 +125,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D // q=6 nr6 = neighborList[n+5*Np]; f6 = dist[nr6]; + + // compute diffusive flux + flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci); + flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci); + flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci); + FluxDiffusive[n+0*Np] = flux_diffusive_x; + FluxDiffusive[n+1*Np] = flux_diffusive_y; + FluxDiffusive[n+2*Np] = flux_diffusive_z; // q=0 dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; @@ -149,13 +158,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D } } -extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; for (n=start; nFirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], + ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np); if (BoundaryConditionSolid==1){ @@ -933,9 +934,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){ //LB-Ion collison - ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], + ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], + ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic], rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np); if (BoundaryConditionSolid==1){ @@ -973,7 +974,22 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const s ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration); ScaLBL_Comm->Barrier(); comm.barrier(); IonConcentration_LB_to_Phys(IonConcentration); +} +void ScaLBL_IonModel::getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic){ + //This function wirte out the data in a normal layout (by aggregating all decomposed domains) + + ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],IonFlux_x); + IonFlux_LB_to_Phys(IonFlux_x,ic); + ScaLBL_Comm->Barrier(); comm.barrier(); + + ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],IonFlux_y); + IonFlux_LB_to_Phys(IonFlux_y,ic); + ScaLBL_Comm->Barrier(); comm.barrier(); + + ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],IonFlux_z); + IonFlux_LB_to_Phys(IonFlux_z,ic); + ScaLBL_Comm->Barrier(); comm.barrier(); } void ScaLBL_IonModel::getIonConcentration_debug(int timestep){ @@ -992,13 +1008,67 @@ void ScaLBL_IonModel::getIonConcentration_debug(int timestep){ } } +void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){ + //This function write out decomposed data + + DoubleArray PhaseField(Nx,Ny,Nz); + for (size_t ic=0; icRegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],PhaseField); + ScaLBL_Comm->Barrier(); comm.barrier(); + IonFlux_LB_to_Phys(PhaseField); + + FILE *OUTFILE_X; + sprintf(LocalRankFilename,"IonFluxDiffusive_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,OUTFILE_X); + fclose(OUTFILE_X); + + //y-component + ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],PhaseField); + ScaLBL_Comm->Barrier(); comm.barrier(); + IonFlux_LB_to_Phys(PhaseField); + + FILE *OUTFILE_Y; + sprintf(LocalRankFilename,"IonFluxDiffusive_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,OUTFILE_Y); + fclose(OUTFILE_Y); + + //z-component + ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],PhaseField); + ScaLBL_Comm->Barrier(); comm.barrier(); + IonFlux_LB_to_Phys(PhaseField); + + FILE *OUTFILE_Z; + sprintf(LocalRankFilename,"IonFluxDiffusive_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,OUTFILE_Z); + fclose(OUTFILE_Z); + } +} + + void ScaLBL_IonModel::IonConcentration_LB_to_Phys(DoubleArray &Den_reg){ for (int k=0;k Date: Wed, 29 Sep 2021 02:10:35 -0400 Subject: [PATCH 3/5] build passes after fixing various dumb syntax bugs;model to be verified --- analysis/ElectroChemistry.cpp | 12 ++++++------ models/IonModel.cpp | 12 ++++++------ models/IonModel.h | 3 +++ 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index b46a7473..5de909f5 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -288,19 +288,19 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P // x-component of diffusive flux sprintf(VisName,"Ion%zu_FluxDiffusive_x",ion+1); //IonFluxDiffusive[3*ion+0]->name = VisName; - ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->name==VisName); + ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+0]->name==VisName); // y-component of diffusive flux sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1); //IonFluxDiffusive[3*ion+1]->name = VisName; - ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->name==VisName); + ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+1]->name==VisName); // z-component of diffusive flux sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1); //IonFluxDiffusive[3*ion+2]->name = VisName; - ASSERT(visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->name==VisName); + ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+2]->name==VisName); - Array& IonFluxData_x = visData[0].vars[4+Ion.numer_ion_species+3*ion+0]->data; - Array& IonFluxData_y = visData[0].vars[4+Ion.numer_ion_species+3*ion+1]->data; - Array& IonFluxData_z = visData[0].vars[4+Ion.numer_ion_species+3*ion+2]->data; + Array& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species+3*ion+0]->data; + Array& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species+3*ion+1]->data; + Array& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species+3*ion+2]->data; Ion.getIonFluxDiffusive(IonFluxDiffusive_x,IonFluxDiffusive_y,IonFluxDiffusive_z,ion); fillData.copy(IonFluxDiffusive_x,IonFluxData_x); fillData.copy(IonFluxDiffusive_y,IonFluxData_y); diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 1f08186a..6fcc3cca 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -1016,33 +1016,33 @@ void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){ //x-component ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],PhaseField); ScaLBL_Comm->Barrier(); comm.barrier(); - IonFlux_LB_to_Phys(PhaseField); + IonFlux_LB_to_Phys(PhaseField,ic); FILE *OUTFILE_X; sprintf(LocalRankFilename,"IonFluxDiffusive_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); - OUTFILE = fopen(LocalRankFilename,"wb"); + OUTFILE_X = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,OUTFILE_X); fclose(OUTFILE_X); //y-component ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],PhaseField); ScaLBL_Comm->Barrier(); comm.barrier(); - IonFlux_LB_to_Phys(PhaseField); + IonFlux_LB_to_Phys(PhaseField,ic); FILE *OUTFILE_Y; sprintf(LocalRankFilename,"IonFluxDiffusive_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); - OUTFILE = fopen(LocalRankFilename,"wb"); + OUTFILE_Y = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,OUTFILE_Y); fclose(OUTFILE_Y); //z-component ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],PhaseField); ScaLBL_Comm->Barrier(); comm.barrier(); - IonFlux_LB_to_Phys(PhaseField); + IonFlux_LB_to_Phys(PhaseField,ic); FILE *OUTFILE_Z; sprintf(LocalRankFilename,"IonFluxDiffusive_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); - OUTFILE = fopen(LocalRankFilename,"wb"); + OUTFILE_Z = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,OUTFILE_Z); fclose(OUTFILE_Z); } diff --git a/models/IonModel.h b/models/IonModel.h index 33e3eb4b..ee39a1c8 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -36,6 +36,8 @@ public: void Run(double *Velocity, double *ElectricField); void getIonConcentration(DoubleArray &IonConcentration, const size_t ic); void getIonConcentration_debug(int timestep); + void getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic); + void getIonFluxDiffusive_debug(int timestep); void DummyFluidVelocity(); void DummyElectricField(); double CalIonDenConvergence(vector &ci_avg_previous); @@ -99,5 +101,6 @@ private: void AssignSolidBoundary(double *ion_solid); void AssignIonConcentration_FromFile(double *Ci,const vector &File_ion,int ic); void IonConcentration_LB_to_Phys(DoubleArray &Den_reg); + void IonFlux_LB_to_Phys(DoubleArray &Den_reg, const size_t ic); }; #endif From 487478550e183bb671604e10bc0a0201a7a84024 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Wed, 29 Sep 2021 21:05:39 -0400 Subject: [PATCH 4/5] minor revision --- tests/TestIonModel.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/tests/TestIonModel.cpp b/tests/TestIonModel.cpp index 622a0e8f..17199487 100644 --- a/tests/TestIonModel.cpp +++ b/tests/TestIonModel.cpp @@ -21,13 +21,11 @@ int main(int argc, char **argv) { // Initialize MPI and error handlers Utilities::startup( argc, argv ); + 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 - // Initialize MPI - Utilities::startup( argc, argv ); - Utilities::MPI comm( MPI_COMM_WORLD ); - int rank = comm.getRank(); - int nprocs = comm.getSize(); if (rank == 0){ printf("**************************************\n"); From 551bcb172f168c15255423c9e0acc249887a342b Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 30 Sep 2021 01:16:53 -0400 Subject: [PATCH 5/5] update saving diffusive flux in GPU version --- cuda/Ion.cu | 30 ++++++++++++++++++++++++------ hip/Ion.cu | 30 ++++++++++++++++++++++++------ 2 files changed, 48 insertions(+), 12 deletions(-) diff --git a/cuda/Ion.cu b/cuda/Ion.cu index 601d5310..6e95b02e 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -97,13 +97,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *D } } -__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; int nr1,nr2,nr3,nr4,nr5,nr6; @@ -146,6 +147,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub nr6 = neighborList[n+5*Np]; f6 = dist[nr6]; + // compute diffusive flux + flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci); + flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci); + flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci); + FluxDiffusive[n+0*Np] = flux_diffusive_x; + FluxDiffusive[n+1*Np] = flux_diffusive_y; + FluxDiffusive[n+2*Np] = flux_diffusive_z; + // q=0 dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; //dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); @@ -177,13 +186,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub } } -__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; int S = Np/NBLOCKS/NTHREADS + 1; @@ -212,6 +222,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *V f5 = dist[6*Np+n]; f6 = dist[5*Np+n]; + // compute diffusive flux + flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci); + flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci); + flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci); + FluxDiffusive[n+0*Np] = flux_diffusive_x; + FluxDiffusive[n+1*Np] = flux_diffusive_y; + FluxDiffusive[n+2*Np] = flux_diffusive_z; + // q=0 dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; //dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); @@ -330,10 +348,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(neighborList,dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -342,10 +360,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ diff --git a/hip/Ion.cu b/hip/Ion.cu index 2c48858d..6428231d 100644 --- a/hip/Ion.cu +++ b/hip/Ion.cu @@ -98,13 +98,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *D } } -__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; int nr1,nr2,nr3,nr4,nr5,nr6; @@ -147,6 +148,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub nr6 = neighborList[n+5*Np]; f6 = dist[nr6]; + // compute diffusive flux + flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci); + flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci); + flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci); + FluxDiffusive[n+0*Np] = flux_diffusive_x; + FluxDiffusive[n+1*Np] = flux_diffusive_y; + FluxDiffusive[n+2*Np] = flux_diffusive_z; + // q=0 dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; //dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); @@ -178,13 +187,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub } } -__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, ddouble *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ int n; double Ci; double ux,uy,uz; double uEPx,uEPy,uEPz;//electrochemical induced velocity double Ex,Ey,Ez;//electrical field + double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z; double f0,f1,f2,f3,f4,f5,f6; int S = Np/NBLOCKS/NTHREADS + 1; @@ -213,6 +223,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *V f5 = dist[6*Np+n]; f6 = dist[5*Np+n]; + // compute diffusive flux + flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci); + flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci); + flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci); + FluxDiffusive[n+0*Np] = flux_diffusive_x; + FluxDiffusive[n+1*Np] = flux_diffusive_y; + FluxDiffusive[n+2*Np] = flux_diffusive_z; + // q=0 dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci; //dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz))); @@ -331,10 +349,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(neighborList,dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ @@ -343,10 +361,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField, +extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField, double Di, int zi, double rlx, double Vt, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){