From afc7d6c90e79035061a5bc9326435832bdc29b45 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Wed, 13 Oct 2021 15:15:12 +1100 Subject: [PATCH 1/4] add routines to save advective and electromigrational flux;CPU only;to be verified --- analysis/ElectroChemistry.cpp | 127 ++++++++++++++++++++++++++++++++++ analysis/ElectroChemistry.h | 3 + common/ScaLBL.h | 4 +- cpu/Ion.cpp | 16 ++++- models/IonModel.cpp | 123 ++++++++++++++++++++++++++++++-- models/IonModel.h | 2 + tests/TestIonModel.cpp | 2 + 7 files changed, 268 insertions(+), 9 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 5de909f5..2e926356 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -10,6 +10,9 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr 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 @@ -169,13 +172,17 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P visData[0].meshName = "domain"; visData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); auto ElectricPotential = std::make_shared(); + + //ion concentration std::vector> IonConcentration; for (size_t ion=0; ion()); } + //fluid velocity auto VxVar = std::make_shared(); auto VyVar = std::make_shared(); auto VzVar = std::make_shared(); + // diffusive ion flux std::vector> IonFluxDiffusive; for (size_t ion=0; ion()); IonFluxDiffusive.push_back(std::make_shared()); } + // advective ion flux + std::vector> IonFluxAdvective; + for (size_t ion=0; ion()); + IonFluxAdvective.push_back(std::make_shared()); + IonFluxAdvective.push_back(std::make_shared()); + } + // electro-migrational ion flux + std::vector> IonFluxElectrical; + for (size_t ion=0; ion()); + IonFluxElectrical.push_back(std::make_shared()); + IonFluxElectrical.push_back(std::make_shared()); + } //-------------------------------------------------------------------------------------------------------------------- //-------------------------------------Create Names for Variables------------------------------------------------------ @@ -248,6 +271,58 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]); } } + + if (vis_db->getWithDefault( "save_ion_flux_advective", false )){ + for (size_t ion=0; ionname = 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( "save_ion_flux_electrical", false )){ + for (size_t ion=0; ionname = 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]); + } + } //-------------------------------------------------------------------------------------------------------------------- //------------------------------------Save All Variables-------------------------------------------------------------- @@ -307,7 +382,59 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P fillData.copy(IonFluxDiffusive_z,IonFluxData_z); } } + + if (vis_db->getWithDefault( "save_ion_flux_advective", false )){ + for (size_t ion=0; ionname = 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& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->data; + Array& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->data; + Array& 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( "save_ion_flux_electrical", false )){ + for (size_t ion=0; ionname = 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& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->data; + Array& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->data; + Array& 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( "write_silo", true )) IO::writeData( timestep, visData, Dm->Comm ); //-------------------------------------------------------------------------------------------------------------------- diff --git a/analysis/ElectroChemistry.h b/analysis/ElectroChemistry.h index ee529725..50e44b3f 100644 --- a/analysis/ElectroChemistry.h +++ b/analysis/ElectroChemistry.h @@ -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; diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 6e15b9f1..ba253be6 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 *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); diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index 672ad71c..dc8ef02e 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -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; diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 6fcc3cca..100063ea 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -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; icRegularLayout(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; icRegularLayout(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 Date: Wed, 13 Oct 2021 00:33:08 -0400 Subject: [PATCH 2/4] fix dumb bugs;build passed --- models/IonModel.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/models/IonModel.h b/models/IonModel.h index 1398d9f6..5f3e390c 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -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 &ci_avg_previous); From 354979a395455648f2c9858cd6ef603da02abf2e Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Wed, 13 Oct 2021 15:46:08 +1100 Subject: [PATCH 3/4] update flux-saving routines on GPU;to be verified --- cuda/Ion.cu | 24 ++++++++++++++++++------ hip/Ion.cu | 24 ++++++++++++++++++------ 2 files changed, 36 insertions(+), 12 deletions(-) diff --git a/cuda/Ion.cu b/cuda/Ion.cu index 6e95b02e..49d0b80a 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -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<<>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(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<<>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,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 6428231d..b1d9636e 100644 --- a/hip/Ion.cu +++ b/hip/Ion.cu @@ -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<<>>(neighborList,dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Ion<<>>(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<<>>(dist,Den,FluxDiffusive,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Ion<<>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ From ead582cae0505b7309a6b49a289d28e4a3c1755d Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Wed, 13 Oct 2021 16:36:04 +1100 Subject: [PATCH 4/4] add WriteVis option to save electric field --- analysis/ElectroChemistry.cpp | 48 ++++++++++++++++++++++++++++++----- 1 file changed, 42 insertions(+), 6 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 2e926356..9eb56315 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -171,7 +171,12 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P visData[0].meshName = "domain"; visData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); - auto ElectricPotential = std::make_shared(); + //electric potential + auto ElectricPotentialVar = std::make_shared(); + //electric field + auto ElectricFieldVar_x = std::make_shared(); + auto ElectricFieldVar_y = std::make_shared(); + auto ElectricFieldVar_z = std::make_shared(); //ion concentration std::vector> IonConcentration; @@ -210,11 +215,11 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P //-------------------------------------Create Names for Variables------------------------------------------------------ if (vis_db->getWithDefault( "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( "save_concentration", true )){ @@ -323,6 +328,24 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P visData[0].vars.push_back(IonFluxElectrical[3*ion+2]); } } + + if (vis_db->getWithDefault( "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-------------------------------------------------------------- @@ -435,6 +458,19 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P } } + if (vis_db->getWithDefault( "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& ElectricalFieldxData = visData[0].vars[4+Ion.number_ion_species*(1+9)+0]->data; + Array& ElectricalFieldyData = visData[0].vars[4+Ion.number_ion_species*(1+9)+1]->data; + Array& 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( "write_silo", true )) IO::writeData( timestep, visData, Dm->Comm ); //--------------------------------------------------------------------------------------------------------------------