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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James McClure 2021-10-14 08:29:10 -04:00
commit f99e2dbef9
9 changed files with 350 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

175
analysis/ElectroChemistry.cpp
View File

@ -10,6 +10,9 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr <Domain> dm):
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
@ -168,14 +171,23 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
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 ElectricPotential = std::make_shared<IO::Variable>();
//electric potential
auto ElectricPotentialVar = std::make_shared<IO::Variable>();
//electric field
auto ElectricFieldVar_x = std::make_shared<IO::Variable>();
auto ElectricFieldVar_y = std::make_shared<IO::Variable>();
auto ElectricFieldVar_z = std::make_shared<IO::Variable>();
//ion concentration
std::vector<shared_ptr<IO::Variable>> IonConcentration;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
IonConcentration.push_back(std::make_shared<IO::Variable>());
}
//fluid velocity
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
// diffusive ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxDiffusive;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
@ -183,15 +195,31 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
}
// advective ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxAdvective;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
}
// electro-migrational ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxElectrical;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
}
//--------------------------------------------------------------------------------------------------------------------
//-------------------------------------Create Names for Variables------------------------------------------------------
if (vis_db->getWithDefault<bool>( "save_electric_potential", true )){
ElectricPotential->name = "ElectricPotential";
ElectricPotential->type = IO::VariableType::VolumeVariable;
ElectricPotential->dim = 1;
ElectricPotential->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricPotential);
ElectricPotentialVar->name = "ElectricPotential";
ElectricPotentialVar->type = IO::VariableType::VolumeVariable;
ElectricPotentialVar->dim = 1;
ElectricPotentialVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricPotentialVar);
}
if (vis_db->getWithDefault<bool>( "save_concentration", true )){
@ -248,6 +276,76 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
IonFluxAdvective[3*ion+0]->name = VisName;
IonFluxAdvective[3*ion+0]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+0]->dim = 1;
IonFluxAdvective[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+0]);
// y-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
IonFluxAdvective[3*ion+1]->name = VisName;
IonFluxAdvective[3*ion+1]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+1]->dim = 1;
IonFluxAdvective[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+1]);
// z-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
IonFluxAdvective[3*ion+2]->name = VisName;
IonFluxAdvective[3*ion+2]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+2]->dim = 1;
IonFluxAdvective[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
IonFluxElectrical[3*ion+0]->name = VisName;
IonFluxElectrical[3*ion+0]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+0]->dim = 1;
IonFluxElectrical[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+0]);
// y-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
IonFluxElectrical[3*ion+1]->name = VisName;
IonFluxElectrical[3*ion+1]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+1]->dim = 1;
IonFluxElectrical[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+1]);
// z-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
IonFluxElectrical[3*ion+2]->name = VisName;
IonFluxElectrical[3*ion+2]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+2]->dim = 1;
IonFluxElectrical[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_electric_field", false )){
ElectricFieldVar_x->name = "ElectricField_x";
ElectricFieldVar_x->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_x->dim = 1;
ElectricFieldVar_x->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_x);
ElectricFieldVar_y->name = "ElectricField_y";
ElectricFieldVar_y->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_y->dim = 1;
ElectricFieldVar_y->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_y);
ElectricFieldVar_z->name = "ElectricField_z";
ElectricFieldVar_z->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_z->dim = 1;
ElectricFieldVar_z->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_z);
}
//--------------------------------------------------------------------------------------------------------------------
//------------------------------------Save All Variables--------------------------------------------------------------
@ -307,7 +405,72 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
fillData.copy(IonFluxDiffusive_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
//IonFluxDiffusive[3*ion+0]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->name==VisName);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
//IonFluxDiffusive[3*ion+1]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->name==VisName);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
//IonFluxDiffusive[3*ion+2]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->name==VisName);
Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->data;
Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->data;
Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->data;
Ion.getIonFluxAdvective(IonFluxAdvective_x,IonFluxAdvective_y,IonFluxAdvective_z,ion);
fillData.copy(IonFluxAdvective_x,IonFluxData_x);
fillData.copy(IonFluxAdvective_y,IonFluxData_y);
fillData.copy(IonFluxAdvective_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
//IonFluxDiffusive[3*ion+0]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->name==VisName);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
//IonFluxDiffusive[3*ion+1]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->name==VisName);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
//IonFluxDiffusive[3*ion+2]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->name==VisName);
Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->data;
Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->data;
Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->data;
Ion.getIonFluxElectrical(IonFluxElectrical_x,IonFluxElectrical_y,IonFluxElectrical_z,ion);
fillData.copy(IonFluxElectrical_x,IonFluxData_x);
fillData.copy(IonFluxElectrical_y,IonFluxData_y);
fillData.copy(IonFluxElectrical_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_electric_field", false )){
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+0]->name=="ElectricField_x");
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+1]->name=="ElectricField_y");
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+2]->name=="ElectricField_z");
Poisson.getElectricField(ElectricalField_x, ElectricalField_y, ElectricalField_z);
Array<double>& ElectricalFieldxData = visData[0].vars[4+Ion.number_ion_species*(1+9)+0]->data;
Array<double>& ElectricalFieldyData = visData[0].vars[4+Ion.number_ion_species*(1+9)+1]->data;
Array<double>& ElectricalFieldzData = visData[0].vars[4+Ion.number_ion_species*(1+9)+2]->data;
fillData.copy(ElectricalField_x,ElectricalFieldxData);
fillData.copy(ElectricalField_y,ElectricalFieldyData);
fillData.copy(ElectricalField_z,ElectricalFieldzData);
}
if (vis_db->getWithDefault<bool>( "write_silo", true ))
IO::writeData( timestep, visData, Dm->Comm );
//--------------------------------------------------------------------------------------------------------------------

3
analysis/ElectroChemistry.h
View File

@ -35,6 +35,9 @@ public:
DoubleArray Rho; // density field
DoubleArray ChemicalPotential; // density field
DoubleArray ElectricalPotential; // density field
DoubleArray ElectricalField_x; // density field
DoubleArray ElectricalField_y; // density field
DoubleArray ElectricalField_z; // density field
DoubleArray Pressure; // pressure field
DoubleArray Vel_x; // velocity field
DoubleArray Vel_y;

4
common/ScaLBL.h
View File

@ -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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, 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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, 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);

16
cpu/Ion.cpp
View File

@ -80,7 +80,7 @@ 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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -133,6 +133,12 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -158,7 +164,7 @@ 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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -197,6 +203,12 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDi
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;

24
cuda/Ion.cu
View File

@ -97,7 +97,7 @@ __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 *FluxDiffusive, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -154,6 +154,12 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -186,7 +192,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -229,6 +235,12 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -348,10 +360,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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@ -360,10 +372,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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){

24
hip/Ion.cu
View File

@ -98,7 +98,7 @@ __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 *FluxDiffusive, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -155,6 +155,12 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -187,7 +193,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, ddouble *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
@ -230,6 +236,12 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -349,10 +361,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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
@ -361,10 +373,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 *FluxDiffusive, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){

123
models/IonModel.cpp
View File

@ -692,7 +692,9 @@ void ScaLBL_IonModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &fq, number_ion_species*7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Ci, number_ion_species*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &ChargeDensity, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxDiffusive, number_ion_species*3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxDiffusive, number_ion_species*3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxAdvective, number_ion_species*3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxElectrical, number_ion_species*3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("LB Ion Solver: Setting up device map and neighbor list \n");
@ -878,9 +880,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//LB-Ion collison
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
@ -934,9 +936,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//LB-Ion collison
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[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],&FluxDiffusive[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
@ -992,6 +994,38 @@ void ScaLBL_IonModel::getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &Io
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_IonModel::getIonFluxAdvective(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,&FluxAdvective[ic*3*Np+0*Np],IonFlux_x);
IonFlux_LB_to_Phys(IonFlux_x,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],IonFlux_y);
IonFlux_LB_to_Phys(IonFlux_y,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],IonFlux_z);
IonFlux_LB_to_Phys(IonFlux_z,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_IonModel::getIonFluxElectrical(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,&FluxElectrical[ic*3*Np+0*Np],IonFlux_x);
IonFlux_LB_to_Phys(IonFlux_x,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],IonFlux_y);
IonFlux_LB_to_Phys(IonFlux_y,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[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){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
@ -1048,6 +1082,85 @@ void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){
}
}
void ScaLBL_IonModel::getIonFluxAdvective_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (size_t ic=0; ic<number_ion_species; ic++){
//x-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+0*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_X;
sprintf(LocalRankFilename,"IonFluxAdvective_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_X = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_X);
fclose(OUTFILE_X);
//y-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Y;
sprintf(LocalRankFilename,"IonFluxAdvective_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Y = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Y);
fclose(OUTFILE_Y);
//z-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Z;
sprintf(LocalRankFilename,"IonFluxAdvective_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Z = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Z);
fclose(OUTFILE_Z);
}
}
void ScaLBL_IonModel::getIonFluxElectrical_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (size_t ic=0; ic<number_ion_species; ic++){
//x-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+0*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_X;
sprintf(LocalRankFilename,"IonFluxElectrical_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_X = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_X);
fclose(OUTFILE_X);
//y-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Y;
sprintf(LocalRankFilename,"IonFluxElectrical_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Y = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Y);
fclose(OUTFILE_Y);
//z-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+2*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Z;
sprintf(LocalRankFilename,"IonFluxElectrical_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Z = 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<Nz;k++){

6
models/IonModel.h
View File

@ -37,7 +37,11 @@ public:
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 getIonFluxAdvective(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic);
void getIonFluxElectrical(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic);
void getIonFluxDiffusive_debug(int timestep);
void getIonFluxAdvective_debug(int timestep);
void getIonFluxElectrical_debug(int timestep);
void DummyFluidVelocity();
void DummyElectricField();
double CalIonDenConvergence(vector<double> &ci_avg_previous);
@ -86,6 +90,8 @@ public:
double *FluidVelocityDummy;
double *ElectricFieldDummy;
double *FluxDiffusive;
double *FluxAdvective;
double *FluxElectrical;
private:
Utilities::MPI comm;

2
tests/TestIonModel.cpp
View File

@ -76,6 +76,8 @@ int main(int argc, char **argv)
}
IonModel.getIonConcentration_debug(timestep);
IonModel.getIonFluxDiffusive_debug(timestep);
IonModel.getIonFluxAdvective_debug(timestep);
IonModel.getIonFluxElectrical_debug(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");