diff --git a/IO/IOHelpers.h b/IO/IOHelpers.h index 231071cd..a7ba6096 100644 --- a/IO/IOHelpers.h +++ b/IO/IOHelpers.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef IO_HELPERS_INC #define IO_HELPERS_INC diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp index e86f749e..5e5eb96a 100644 --- a/IO/Mesh.cpp +++ b/IO/Mesh.cpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "Mesh.h" #include "IO/IOHelpers.h" #include "common/Utilities.h" diff --git a/IO/Mesh.h b/IO/Mesh.h index c2be3d17..9e5f32e6 100644 --- a/IO/Mesh.h +++ b/IO/Mesh.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef MESH_INC #define MESH_INC diff --git a/IO/MeshDatabase.cpp b/IO/MeshDatabase.cpp index f5a5f10c..975c9e27 100644 --- a/IO/MeshDatabase.cpp +++ b/IO/MeshDatabase.cpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "IO/MeshDatabase.h" #include "IO/IOHelpers.h" #include "IO/Mesh.h" diff --git a/IO/MeshDatabase.h b/IO/MeshDatabase.h index 0a14aeca..508f85d8 100644 --- a/IO/MeshDatabase.h +++ b/IO/MeshDatabase.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef MeshDatabase_INC #define MeshDatabase_INC diff --git a/IO/PIO.cpp b/IO/PIO.cpp index a9b6dd43..f959cb49 100644 --- a/IO/PIO.cpp +++ b/IO/PIO.cpp @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "IO/PIO.h" #include "common/MPI.h" #include "common/Utilities.h" diff --git a/IO/PIO.h b/IO/PIO.h index 7906cb85..9b8aeb89 100644 --- a/IO/PIO.h +++ b/IO/PIO.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_PIO #define included_PIO diff --git a/IO/PIO.hpp b/IO/PIO.hpp index e645b2f3..748bf32b 100644 --- a/IO/PIO.hpp +++ b/IO/PIO.hpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_PIO_hpp #define included_PIO_hpp diff --git a/IO/Reader.cpp b/IO/Reader.cpp index 743f0e3f..566c07cb 100644 --- a/IO/Reader.cpp +++ b/IO/Reader.cpp @@ -1,24 +1,8 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "IO/Reader.h" #include "IO/HDF5_IO.h" #include "IO/IOHelpers.h" #include "IO/Mesh.h" #include "IO/MeshDatabase.h" -#include "IO/silo.h" #include "common/Utilities.h" #include diff --git a/IO/Reader.h b/IO/Reader.h index 80afb116..9dfc834b 100644 --- a/IO/Reader.h +++ b/IO/Reader.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef READER_INC #define READER_INC diff --git a/IO/Writer.cpp b/IO/Writer.cpp index a0b6a88a..0fb8d135 100644 --- a/IO/Writer.cpp +++ b/IO/Writer.cpp @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "IO/Writer.h" #include "IO/HDF5_IO.h" #include "IO/IOHelpers.h" diff --git a/IO/Writer.h b/IO/Writer.h index c6db8978..3844f3b2 100644 --- a/IO/Writer.h +++ b/IO/Writer.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef WRITER_INC #define WRITER_INC diff --git a/IO/silo.cpp b/IO/silo.cpp index 333f4831..1a264eb6 100644 --- a/IO/silo.cpp +++ b/IO/silo.cpp @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "IO/silo.h" #include "common/MPI.h" #include "common/Utilities.h" diff --git a/IO/silo.h b/IO/silo.h index e13d5a48..aef994df 100644 --- a/IO/silo.h +++ b/IO/silo.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef SILO_INTERFACE #define SILO_INTERFACE diff --git a/IO/silo.hpp b/IO/silo.hpp index 56e6c985..d86ad17f 100644 --- a/IO/silo.hpp +++ b/IO/silo.hpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef SILO_INTERFACE_HPP #define SILO_INTERFACE_HPP diff --git a/StackTrace/StackTrace.h b/StackTrace/StackTrace.h index 7aab2d72..3773509c 100644 --- a/StackTrace/StackTrace.h +++ b/StackTrace/StackTrace.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_StackTrace #define included_StackTrace diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 4f23f00f..b4c4a504 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -141,6 +141,12 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel) for (size_t i=0; iinlet_layers_z; k < Nz; k++) { for (j = 1; j < Ny; j++) { for (i = 1; i < Nx; i++) { /* electric potential */ @@ -258,29 +279,58 @@ void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion, value_out_local = 0.0; for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { Ion.getIonConcentration(Rho, ion); + Ion.getIonFluxDiffusive(IonFluxDiffusive_x, IonFluxDiffusive_y, IonFluxDiffusive_z, ion); + Ion.getIonFluxAdvective(IonFluxAdvective_x, IonFluxAdvective_y, IonFluxAdvective_z, ion); + Ion.getIonFluxElectrical(IonFluxElectrical_x, IonFluxElectrical_y, IonFluxElectrical_z, ion); + value_membrane_in_local = 0.0; value_membrane_out_local = 0.0; value_in_local = 0.0; value_out_local = 0.0; - for (k = 1; k < Nz; k++) { + + jx_membrane_in_local = jy_membrane_in_local = jz_membrane_in_local = 0.0; + jx_membrane_out_local = jy_membrane_out_local = jz_membrane_out_local = 0.0; + jx_in_local = jy_in_local = jz_in_local = 0.0; + jx_out_local = jy_out_local = jz_out_local = 0.0; + + for (k = Dm->inlet_layers_z; k < Nz; k++) { for (j = 1; j < Ny; j++) { for (i = 1; i < Nx; i++) { /* electric potential */ memdist = Ion.MembraneDistance(i,j,k); value = Rho(i,j,k); + jx = IonFluxDiffusive_x(i,j,k) + IonFluxAdvective_x(i,j,k) + IonFluxElectrical_x(i,j,k); + jy = IonFluxDiffusive_y(i,j,k) + IonFluxAdvective_y(i,j,k) + IonFluxElectrical_y(i,j,k); + jz = IonFluxDiffusive_z(i,j,k) + IonFluxAdvective_z(i,j,k) + IonFluxElectrical_z(i,j,k); + if (memdist < 0.0){ // inside the membrane if (fabs(memdist) < 1.0){ value_membrane_in_local += value; + jx_membrane_in_local += jx; + jy_membrane_in_local += jy; + jz_membrane_in_local += jz; + } value_in_local += value; + jx_in_local += jx; + jy_in_local += jy; + jz_in_local += jz; + } else { // outside the membrane if (fabs(memdist) < 1.0){ value_membrane_out_local += value; + jx_membrane_out_local += jx; + jy_membrane_out_local += jy; + jz_membrane_out_local += jz; + } value_out_local += value; + jx_out_local += jx; + jy_out_local += jy; + jz_out_local += jz; } } } @@ -294,12 +344,57 @@ void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion, value_in_global /= in_global_count; value_membrane_out_global /= membrane_out_global_count; value_membrane_in_global /= membrane_in_global_count; + + jx_out_global = Dm->Comm.sumReduce(jx_out_local); + jx_in_global = Dm->Comm.sumReduce(jx_in_local); + jx_membrane_out_global = Dm->Comm.sumReduce(jx_membrane_out_local); + jx_membrane_in_global = Dm->Comm.sumReduce(jx_membrane_in_local); + + jx_out_global /= out_global_count; + jx_in_global /= in_global_count; + jx_membrane_out_global /= membrane_out_global_count; + jx_membrane_in_global /= membrane_in_global_count; + + jy_out_global = Dm->Comm.sumReduce(jy_out_local); + jy_in_global = Dm->Comm.sumReduce(jy_in_local); + jy_membrane_out_global = Dm->Comm.sumReduce(jy_membrane_out_local); + jy_membrane_in_global = Dm->Comm.sumReduce(jy_membrane_in_local); + + jy_out_global /= out_global_count; + jy_in_global /= in_global_count; + jy_membrane_out_global /= membrane_out_global_count; + jy_membrane_in_global /= membrane_in_global_count; + + jz_out_global = Dm->Comm.sumReduce(jz_out_local); + jz_in_global = Dm->Comm.sumReduce(jz_in_local); + jz_membrane_out_global = Dm->Comm.sumReduce(jz_membrane_out_local); + jz_membrane_in_global = Dm->Comm.sumReduce(jz_membrane_in_local); + + jz_out_global /= out_global_count; + jz_in_global /= in_global_count; + jz_membrane_out_global /= membrane_out_global_count; + jz_membrane_in_global /= membrane_in_global_count; if (Dm->rank() == 0) { fprintf(TIMELOG, "%.8g ", value_out_global); fprintf(TIMELOG, "%.8g ", value_in_global); fprintf(TIMELOG, "%.8g ", value_membrane_out_global); fprintf(TIMELOG, "%.8g ", value_membrane_in_global); + + fprintf(TIMELOG, "%.8g ", jx_out_global); + fprintf(TIMELOG, "%.8g ", jx_in_global); + fprintf(TIMELOG, "%.8g ", jx_membrane_out_global); + fprintf(TIMELOG, "%.8g ", jx_membrane_in_global); + + fprintf(TIMELOG, "%.8g ", jy_out_global); + fprintf(TIMELOG, "%.8g ", jy_in_global); + fprintf(TIMELOG, "%.8g ", jy_membrane_out_global); + fprintf(TIMELOG, "%.8g ", jy_membrane_in_global); + + fprintf(TIMELOG, "%.8g ", jz_out_global); + fprintf(TIMELOG, "%.8g ", jz_in_global); + fprintf(TIMELOG, "%.8g ", jz_membrane_out_global); + fprintf(TIMELOG, "%.8g ", jz_membrane_in_global); } } @@ -427,10 +522,35 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion, ScaLBL_StokesModel &Stokes, std::shared_ptr input_db, int timestep) { - + auto vis_db = input_db->getDatabase("Visualization"); char VisName[40]; auto format = vis_db->getWithDefault( "format", "hdf5" ); + + if (Dm->rank() == 0) { + printf("ElectroChemistryAnalyzer::WriteVis (format = %s)\n", format.c_str()); + if (vis_db->getWithDefault("save_electric_potential", true)){ + printf(" save electric potential \n"); + } + if (vis_db->getWithDefault("save_concentration", true)) { + printf(" save concentration \n"); + } + if (vis_db->getWithDefault("save_velocity", false)) { + printf(" save velocity \n"); + } + if (vis_db->getWithDefault("save_ion_flux_diffusive", false)) { + printf(" save ion flux (diffusive) \n"); + } + if (vis_db->getWithDefault("save_ion_flux_advective", false)) { + printf(" save ion flux (advective) \n"); + } + if (vis_db->getWithDefault("save_ion_flux_electrical", false)) { + printf(" save ion flux (electrical) \n"); + } + if (vis_db->getWithDefault("save_electric_field", false)) { + printf(" save electric field \n"); + } + } std::vector visData; fillHalo fillData(Dm->Comm, Dm->rank_info, @@ -961,6 +1081,31 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion, fillHalo fillData(Dm->Comm, Dm->rank_info, {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1}, 0, 1); + + if (Dm->rank() == 0) { + printf("ElectroChemistryAnalyzer::WriteVis (format = %s)\n", format.c_str()); + if (vis_db->getWithDefault("save_electric_potential", true)){ + printf(" save electric potential \n"); + } + if (vis_db->getWithDefault("save_concentration", true)) { + printf(" save concentration \n"); + } + if (vis_db->getWithDefault("save_velocity", false)) { + printf(" save velocity \n"); + } + if (vis_db->getWithDefault("save_ion_flux_diffusive", false)) { + printf(" save ion flux (diffusive) \n"); + } + if (vis_db->getWithDefault("save_ion_flux_advective", false)) { + printf(" save ion flux (advective) \n"); + } + if (vis_db->getWithDefault("save_ion_flux_electrical", false)) { + printf(" save ion flux (electrical) \n"); + } + if (vis_db->getWithDefault("save_electric_field", false)) { + printf(" save electric field \n"); + } + } IO::initialize("",format,"false"); // Create the MeshDataStruct diff --git a/analysis/FlowAdaptor.h b/analysis/FlowAdaptor.h index 37addd68..b2084726 100644 --- a/analysis/FlowAdaptor.h +++ b/analysis/FlowAdaptor.h @@ -62,7 +62,7 @@ public: * \details Update fractional flow condition. Mass will be preferentially added or removed from * phase regions based on where flow is occurring * @param M ScaLBL_ColorModel - */ + */ double UpdateFractionalFlow(ScaLBL_ColorModel &M); /** diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index e6bfb55d..59700d77 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/Minkowski.h" #include "analysis/pmmc.h" #include "analysis/analysis.h" diff --git a/analysis/Minkowski.h b/analysis/Minkowski.h index 9ecd324a..d14b0665 100644 --- a/analysis/Minkowski.h +++ b/analysis/Minkowski.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Header file for two-phase averaging class #ifndef Minkowski_INC #define Minkowski_INC diff --git a/analysis/PointList.h b/analysis/PointList.h index 21a813ac..397309c9 100644 --- a/analysis/PointList.h +++ b/analysis/PointList.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef PointList_INC #define PointList_INC diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 0b7cd3a5..dafd31fd 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -430,11 +430,11 @@ void SubPhase::Basic() { //double total_flow_rate = water_flow_rate + not_water_flow_rate; //double fractional_flow = water_flow_rate / total_flow_rate; double h = Dm->voxel_length; - double krn = h * h * nu_n * Porosity* Porosity * not_water_flow_rate / force_mag; - double krw = h * h * nu_w * Porosity* Porosity* water_flow_rate / force_mag; + double krn = h * h * nu_n * Porosity * not_water_flow_rate / force_mag; + double krw = h * h * nu_w * Porosity* water_flow_rate / force_mag; /* not counting films */ - double krnf = krn - h * h * nu_n * Porosity* Porosity * not_water_film_flow_rate / force_mag; - double krwf = krw - h * h * nu_w * Porosity* Porosity * water_film_flow_rate / force_mag; + double krnf = krn - h * h * nu_n * Porosity * not_water_film_flow_rate / force_mag; + double krwf = krw - h * h * nu_w * Porosity * water_film_flow_rate / force_mag; double eff_pressure = 1.0 / (krn + krw); // effective pressure drop fprintf(TIMELOG, diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index ebebf222..5670be6b 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/TwoPhase.h" #include "analysis/pmmc.h" diff --git a/analysis/TwoPhase.h b/analysis/TwoPhase.h index 15247431..700b435a 100644 --- a/analysis/TwoPhase.h +++ b/analysis/TwoPhase.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Header file for two-phase averaging class #ifndef TwoPhase_INC #define TwoPhase_INC diff --git a/analysis/analysis.cpp b/analysis/analysis.cpp index 0d2b74fa..0847811d 100644 --- a/analysis/analysis.cpp +++ b/analysis/analysis.cpp @@ -1,20 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - #include "analysis/analysis.h" #include "ProfilerApp.h" diff --git a/analysis/analysis.h b/analysis/analysis.h index 2bb47558..6a729983 100644 --- a/analysis/analysis.h +++ b/analysis/analysis.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef Analysis_H_INC #define Analysis_H_INC diff --git a/analysis/distance.cpp b/analysis/distance.cpp index cf7f1e26..d67193e0 100644 --- a/analysis/distance.cpp +++ b/analysis/distance.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/distance.h" /****************************************************************** diff --git a/analysis/distance.h b/analysis/distance.h index 6043ea56..291bece9 100644 --- a/analysis/distance.h +++ b/analysis/distance.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef Distance_H_INC #define Distance_H_INC diff --git a/analysis/filters.cpp b/analysis/filters.cpp index 22e6a3b4..a262f088 100644 --- a/analysis/filters.cpp +++ b/analysis/filters.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/filters.h" #include "math.h" #include "ProfilerApp.h" diff --git a/analysis/filters.h b/analysis/filters.h index 74d42352..f19a1a90 100644 --- a/analysis/filters.h +++ b/analysis/filters.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef Filters_H_INC #define Filters_H_INC diff --git a/analysis/histogram.h b/analysis/histogram.h index b2cd53d9..fe27f864 100644 --- a/analysis/histogram.h +++ b/analysis/histogram.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ /* * Generate a histogram for volumetric, interfacial and common curve properties * copyright 2014, James E. McClure diff --git a/analysis/imfilter.h b/analysis/imfilter.h index 3390aba8..e6f607b3 100644 --- a/analysis/imfilter.h +++ b/analysis/imfilter.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // These functions mimic the behavior of imfilter in MATLAB #ifndef included_imfilter #define included_imfilter diff --git a/analysis/imfilter.hpp b/analysis/imfilter.hpp index 4272a47e..d2c40d0d 100644 --- a/analysis/imfilter.hpp +++ b/analysis/imfilter.hpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/imfilter.h" #include "ProfilerApp.h" #include diff --git a/analysis/pmmc.h b/analysis/pmmc.h index 5e07cbe8..f2152f0d 100644 --- a/analysis/pmmc.h +++ b/analysis/pmmc.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef pmmc_INC #define pmmc_INC diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 79ff0834..5d1d9c2e 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Run the analysis, blob identification, and write restart files #include "analysis/runAnalysis.h" #include "analysis/analysis.h" diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index bff537f1..15baa00b 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef RunAnalysis_H_INC #define RunAnalysis_H_INC diff --git a/analysis/uCT.cpp b/analysis/uCT.cpp index d7ca0886..feb9bf02 100644 --- a/analysis/uCT.cpp +++ b/analysis/uCT.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "analysis/uCT.h" #include "analysis/analysis.h" #include "analysis/distance.h" diff --git a/analysis/uCT.h b/analysis/uCT.h index 63c3f83b..3d39d806 100644 --- a/analysis/uCT.h +++ b/analysis/uCT.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef uCT_H_INC #define uCT_H_INC diff --git a/common/Array.h b/common/Array.h index 423a33f6..2dd7a785 100644 --- a/common/Array.h +++ b/common/Array.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_ArrayClass #define included_ArrayClass diff --git a/common/Array.hpp b/common/Array.hpp index 971f8663..df56c2b6 100644 --- a/common/Array.hpp +++ b/common/Array.hpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_ArrayClass_hpp #define included_ArrayClass_hpp diff --git a/common/Communication.cpp b/common/Communication.cpp index 22ef57f6..08d98e02 100644 --- a/common/Communication.cpp +++ b/common/Communication.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/Communication.h" /******************************************************** diff --git a/common/Communication.h b/common/Communication.h index 045c20c7..ed549395 100644 --- a/common/Communication.h +++ b/common/Communication.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef COMMUNICATION_H_INC #define COMMUNICATION_H_INC diff --git a/common/Communication.hpp b/common/Communication.hpp index 20af2143..89b5c02e 100644 --- a/common/Communication.hpp +++ b/common/Communication.hpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef COMMUNICATION_HPP_INC #define COMMUNICATION_HPP_INC diff --git a/common/Database.cpp b/common/Database.cpp index 8ef80914..b4c2c905 100644 --- a/common/Database.cpp +++ b/common/Database.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/Database.h" #include "common/Utilities.h" diff --git a/common/Database.h b/common/Database.h index 0f14f360..812b32b6 100644 --- a/common/Database.h +++ b/common/Database.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_Database #define included_Database diff --git a/common/Database.hpp b/common/Database.hpp index a2162b25..a08a80ff 100644 --- a/common/Database.hpp +++ b/common/Database.hpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_Database_hpp #define included_Database_hpp diff --git a/common/Domain.cpp b/common/Domain.cpp index d5411464..6ee0920e 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Created by James McClure // Copyright 2008-2020 #include @@ -355,20 +339,7 @@ void Domain::initialize(std::shared_ptr db) { // Initialize ranks int myrank = Comm.getRank(); rank_info = RankInfoStruct(myrank, nproc[0], nproc[1], nproc[2]); - // inlet layers only apply to lower part of domain - if (rank_info.ix > 0) - inlet_layers_x = 0; - if (rank_info.jy > 0) - inlet_layers_y = 0; - if (rank_info.kz > 0) - inlet_layers_z = 0; - // outlet layers only apply to top part of domain - if (rank_info.ix < nproc[0] - 1) - outlet_layers_x = 0; - if (rank_info.jy < nproc[1] - 1) - outlet_layers_y = 0; - if (rank_info.kz < nproc[2] - 1) - outlet_layers_z = 0; + // Fill remaining variables N = Nx * Ny * Nz; Volume = nx * ny * nz * nproc[0] * nproc[1] * nproc[2] * 1.0; @@ -907,6 +878,22 @@ void Domain::Decomp(const std::string &Filename) { } delete[] SegData; } + /************************/ + // inlet layers only apply to lower part of domain + if (rank_info.ix > 0) + inlet_layers_x = 0; + if (rank_info.jy > 0) + inlet_layers_y = 0; + if (rank_info.kz > 0) + inlet_layers_z = 0; + // outlet layers only apply to top part of domain + if (rank_info.ix < nproc[0] - 1) + outlet_layers_x = 0; + if (rank_info.jy < nproc[1] - 1) + outlet_layers_y = 0; + if (rank_info.kz < nproc[2] - 1) + outlet_layers_z = 0; + /************************/ Comm.barrier(); ComputePorosity(); } diff --git a/common/Domain.h b/common/Domain.h index dfdc4fd8..226c8008 100644 --- a/common/Domain.h +++ b/common/Domain.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef Domain_INC #define Domain_INC diff --git a/common/FunctionTable.h b/common/FunctionTable.h index 7a45880b..c5497848 100644 --- a/common/FunctionTable.h +++ b/common/FunctionTable.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_FunctionTable #define included_FunctionTable diff --git a/common/FunctionTable.hpp b/common/FunctionTable.hpp index e8d6dd33..bcdae59f 100644 --- a/common/FunctionTable.hpp +++ b/common/FunctionTable.hpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_FunctionTable_hpp #define included_FunctionTable_hpp diff --git a/common/MPI.h b/common/MPI.h index d249d661..de3c0534 100644 --- a/common/MPI.h +++ b/common/MPI.h @@ -2,23 +2,6 @@ // Note this is a modified version of the MPI class for the Advanced Multi-Physics Package // Used with permission -/* - -Copyright (c) 2012 UT-Battelle, LLC - -All rights reserved. - -Redistribution and use in source and binary forms, with or without modification, are permitted -provided that the following conditions are met: Redistributions of source code must retain the above -copyright notice, this list of conditions and the following disclaimer. Redistributions in binary -form must reproduce the above copyright notice, this list of conditions and the following disclaimer -in the documentation and/or other materials provided with the distribution. Collection of -administrative costs for redistribution of the source code or binary form is allowed. However, -collection of a royalty or other fee in excess of good faith amount for cost recovery for such -redistribution is prohibited. - -*/ - #ifndef included_LBPM_MPI #define included_LBPM_MPI diff --git a/common/Membrane.cpp b/common/Membrane.cpp index c9943b42..50c7deba 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -12,6 +12,8 @@ Membrane::Membrane(std::shared_ptr sComm, int *dvcNeighbor //...................................................................................... // Create a separate copy of the communicator for the device MPI_COMM_SCALBL = sComm->MPI_COMM_SCALBL.dup(); + int myrank = sComm->MPI_COMM_SCALBL.getRank(); + rank_info = RankInfoStruct(myrank, rank_info.nx, rank_info.ny, rank_info.nz); ScaLBL_CopyToHost(initialNeighborList, dvcNeighborList, 18*Np*sizeof(int)); sComm->MPI_COMM_SCALBL.barrier(); @@ -37,6 +39,8 @@ Membrane::Membrane(std::shared_ptr sComm, int *dvcNeighbor rank_X=sComm->rank_X; rank_Y=sComm->rank_Y; rank_Z=sComm->rank_Z; + + BoundaryCondition = sComm->BoundaryCondition; if (rank == 0){ printf("**** Creating membrane data structure ****** \n"); @@ -828,7 +832,7 @@ void Membrane::SendD3Q7AA(double *dist){ } -void Membrane::RecvD3Q7AA(double *dist){ +void Membrane::RecvD3Q7AA(double *dist, bool BounceBack){ //................................................................................... // Wait for completion of D3Q19 communication @@ -840,7 +844,7 @@ void Membrane::RecvD3Q7AA(double *dist){ // Unpack the distributions on the device //................................................................................... //...Unpacking for x face(q=2)................................ - ScaLBL_D3Q7_Membrane_Unpack(2,dvcRecvDist_x, recvbuf_x,recvCount_x,dist,Np,coefficient_x); + ScaLBL_D3Q7_Membrane_Unpack(2,dvcRecvDist_x, recvbuf_x,recvCount_x,dist,Np,coefficient_x); //................................................................................... //...Packing for X face(q=1)................................ ScaLBL_D3Q7_Membrane_Unpack(1,dvcRecvDist_X, recvbuf_X,recvCount_X,dist,Np,coefficient_X); @@ -851,12 +855,21 @@ void Membrane::RecvD3Q7AA(double *dist){ //...Packing for Y face(q=3)................................. ScaLBL_D3Q7_Membrane_Unpack(3,dvcRecvDist_Y, recvbuf_Y,recvCount_Y,dist,Np,coefficient_Y); //................................................................................... - //...Packing for z face(q=6)................................ - ScaLBL_D3Q7_Membrane_Unpack(6,dvcRecvDist_z, recvbuf_z, recvCount_z,dist,Np,coefficient_z); - //...Packing for Z face(q=5)................................ - ScaLBL_D3Q7_Membrane_Unpack(5,dvcRecvDist_Z, recvbuf_Z,recvCount_Z,dist,Np,coefficient_Z); - //.................................................................................. - + //if (BoundaryCondition > 0 && rank_info.kz == 0) + if (BounceBack && rank_info.kz == 0) + {/* leave the bounce-back distributions in place */} + else { + //...Packing for z face(q=6)................................ + ScaLBL_D3Q7_Membrane_Unpack(6,dvcRecvDist_z, recvbuf_z, recvCount_z,dist,Np,coefficient_z); + } + //if (BoundaryCondition > 0 && rank_info.kz == rank_info.nz-1) + if (BounceBack && rank_info.kz == rank_info.nz-1) + {/* leave the bounce-back distributions in place */} + else { + //...Packing for Z face(q=5)................................ + ScaLBL_D3Q7_Membrane_Unpack(5,dvcRecvDist_Z, recvbuf_Z,recvCount_Z,dist,Np,coefficient_Z); + //.................................................................................. + } MPI_COMM_SCALBL.barrier(); //................................................................................... Lock=false; // unlock the communicator after communications complete diff --git a/common/Membrane.h b/common/Membrane.h index d84c408a..fbfd86cd 100644 --- a/common/Membrane.h +++ b/common/Membrane.h @@ -53,7 +53,8 @@ public: int Np; int Nx,Ny,Nz,N; int membraneLinkCount; - + int BoundaryCondition; + int *initialNeighborList; // original neighborlist int *NeighborList; // modified neighborlist @@ -106,7 +107,7 @@ public: void Read(string filename); void SendD3Q7AA(double *dist); - void RecvD3Q7AA(double *dist); + void RecvD3Q7AA(double *dist, bool BounceBack); void AssignCoefficients(int *Map, double *Psi, double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut); diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index a865b4aa..c513df4f 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2022 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/ScaLBL.h" #include @@ -23,6 +7,9 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ //...................................................................................... // Create a separate copy of the communicator for the device MPI_COMM_SCALBL = Dm->Comm.dup(); + int myrank = MPI_COMM_SCALBL.getRank(); + rank_info = RankInfoStruct(myrank, rank_info.nx, rank_info.ny, rank_info.nz); + //...................................................................................... // Copy the domain size and communication information directly from Dm Nx = Dm->Nx; @@ -1735,18 +1722,7 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){ ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,dist,N); ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Y,4*recvCount_Y,recvCount_Y,recvbuf_Y,dist,N); //................................................................................... - //...Packing for z face(6,12,13,16,17)................................ - ScaLBL_D3Q19_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,dist,N); - ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,dist,N); - ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,dist,N); - ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,recvbuf_z,dist,N); - ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,recvbuf_z,dist,N); - //...Packing for Z face(5,11,14,15,18)................................ - ScaLBL_D3Q19_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,dist,N); - ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); - ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); - ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); - ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); + //.................................................................................. //...Pack the xy edge (8)................................ ScaLBL_D3Q19_Unpack(8,dvcRecvDist_xy,0,recvCount_xy,recvbuf_xy,dist,N); @@ -1756,22 +1732,42 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){ ScaLBL_D3Q19_Unpack(10,dvcRecvDist_xY,0,recvCount_xY,recvbuf_xY,dist,N); //...Pack the XY edge (7)................................ ScaLBL_D3Q19_Unpack(7,dvcRecvDist_XY,0,recvCount_XY,recvbuf_XY,dist,N); - //...Pack the xz edge (12)................................ - ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,dist,N); - //...Pack the xZ edge (14)................................ - ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,recvbuf_xZ,dist,N); - //...Pack the Xz edge (13)................................ - ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,recvbuf_Xz,dist,N); - //...Pack the XZ edge (11)................................ - ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,recvbuf_XZ,dist,N); - //...Pack the yz edge (16)................................ - ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,recvbuf_yz,dist,N); - //...Pack the yZ edge (18)................................ - ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,recvbuf_yZ,dist,N); - //...Pack the Yz edge (17)................................ - ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,dist,N); - //...Pack the YZ edge (15)................................ - ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,dist,N); + + //if (BoundaryCondition == 0 || kproc != 0 ){ + ScaLBL_D3Q19_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,dist,N); + ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,dist,N); + ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,dist,N); + ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,recvbuf_z,dist,N); + ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,recvbuf_z,dist,N); + + //...Pack the xz edge (12)................................ + ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,dist,N); + //...Pack the Xz edge (13)................................ + ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,recvbuf_Xz,dist,N); + //...Pack the yz edge (16)................................ + ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,recvbuf_yz,dist,N); + //...Pack the Yz edge (17)................................ + ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,dist,N); + + //} + //if (BoundaryCondition == 0 || kproc != nprocz-1){ + //...Packing for Z face(5,11,14,15,18)................................ + ScaLBL_D3Q19_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,dist,N); + ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); + ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); + ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); + ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N); + + //...Pack the xZ edge (14)................................ + ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,recvbuf_xZ,dist,N); + //...Pack the XZ edge (11)................................ + ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,recvbuf_XZ,dist,N); + //...Pack the yZ edge (18)................................ + ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,recvbuf_yZ,dist,N); + //...Pack the YZ edge (15)................................ + ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,dist,N); + //} + //................................................................................... Lock=false; // unlock the communicator after communications complete //................................................................................... diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 64d178d5..57b1ef20 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ /** @file ScaLBL.h */ /* \details Header file for Scalable Lattice Boltzmann Library * Separate implementations for GPU and CPU must both follow the conventions defined in this header @@ -217,6 +201,39 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int */ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz); +/** +* \brief BGK collision based on AA even access pattern for D3Q19 +* @param dist - D3Q19 distributions +* @param start - lattice node to start loop +* @param finish - lattice node to finish loop +* @param Np - size of local sub-domain (derived from Domain structure) +* @param rlx - relaxation parameter +* @param Fx - force in x direction +* @param Fy - force in y direction +* @param Fz - force in z direction +*/ +extern "C" void ScaLBL_D3Q19_AAeven_Kubo(double *dist, double *Integral, int start, int finish, int Np); +/** +* \brief Kubo integral function +* @param neighborList - neighbors based on D3Q19 lattice structure +* @param dist - D3Q19 distributions +* @param integral - time integral +* @param start - lattice node to start loop +* @param finish - lattice node to finish loop +* @param Np - size of local sub-domain (derived from Domain structure) +*/ + +extern "C" void ScaLBL_D3Q19_AAodd_Kubo(int *neighborList, double *dist, double *Integral, int start, int finish, int Np); +/** +* \brief Kubo integral function +* @param neighborList - neighbors based on D3Q19 lattice structure +* @param dist - D3Q19 distributions +* @param integral - time integral +* @param start - lattice node to start loop +* @param finish - lattice node to finish loop +* @param Np - size of local sub-domain (derived from Domain structure) +*/ + // MEMBRANE MODEL extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, double *dist, double *Den, int memLinks, int Np); @@ -308,6 +325,15 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np) extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q7_AAeven_pH_ionization( double *dist, + double *Den, double *ElectricField, double * Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np); + +extern "C" void ScaLBL_D3Q7_AAodd_pH_ionization(int *neighborList, double *dist, + double *Den, double *ElectricField, double *Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np); // LBM Poisson solver @@ -383,16 +409,17 @@ extern "C" void ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *Psi, i extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, - double *dist, double *Den_charge, - double *Psi, double *ElectricField, - double tau, double epsilon_LB, bool UseSlippingVelBC, - int start, int finish, int Np); + double *dist, double *Den_charge, + double *Psi, double *ElectricField, + double tau, double Vt, double Cp, + double epsilon_LB, bool UseSlippingVelBC, + int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, - double *Den_charge, double *Psi, - double *ElectricField, double *Error, double tau, - double epsilon_LB, bool UseSlippingVelBC, - int start, int finish, int Np); + double *Den_charge, double *Psi, double *ElectricField, double *Error, + double tau, double Vt, double Cp, + double epsilon_LB, bool UseSlippingVelBC, + int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np); @@ -735,6 +762,8 @@ public: //...................................................................................... unsigned long int CommunicationCount,SendCount,RecvCount; int Nx,Ny,Nz,N; + int iproc,jproc,kproc; + int nprocx,nprocy,nprocz; int n_bb_d3q7, n_bb_d3q19; int BoundaryCondition; //...................................................................................... @@ -838,8 +867,6 @@ private: // only one set of Send requests can be active at any time (per instance) int i,j,k,n; - int iproc,jproc,kproc; - int nprocx,nprocy,nprocz; int sendtag,recvtag; // Give the object it's own MPI communicator RankInfoStruct rank_info; diff --git a/common/SpherePack.cpp b/common/SpherePack.cpp index 12551238..53785fa9 100644 --- a/common/SpherePack.cpp +++ b/common/SpherePack.cpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include #include diff --git a/common/SpherePack.h b/common/SpherePack.h index 2fbdfca3..1ab7fbaa 100644 --- a/common/SpherePack.h +++ b/common/SpherePack.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef SpherePack_INC #define SpherePack_INC diff --git a/common/UnitTest.cpp b/common/UnitTest.cpp index 9ce551e1..71e80464 100644 --- a/common/UnitTest.cpp +++ b/common/UnitTest.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/UnitTest.h" #include "common/Utilities.h" #include diff --git a/common/UnitTest.h b/common/UnitTest.h index 693c9c72..cb169cb2 100644 --- a/common/UnitTest.h +++ b/common/UnitTest.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_UnitTest #define included_UnitTest diff --git a/common/Units.cpp b/common/Units.cpp index 13b2d1e3..d8df428c 100644 --- a/common/Units.cpp +++ b/common/Units.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/Units.h" #include "common/Utilities.h" diff --git a/common/Units.h b/common/Units.h index cdbd25ca..56f587b2 100644 --- a/common/Units.h +++ b/common/Units.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_Units #define included_Units diff --git a/common/Utilities.cpp b/common/Utilities.cpp index 70ca19a4..03ce113f 100644 --- a/common/Utilities.cpp +++ b/common/Utilities.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "common/Utilities.h" #include "StackTrace/StackTrace.h" #include "StackTrace/ErrorHandlers.h" diff --git a/common/Utilities.h b/common/Utilities.h index e9567535..b8d1b760 100644 --- a/common/Utilities.h +++ b/common/Utilities.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_Utilities #define included_Utilities diff --git a/common/UtilityMacros.h b/common/UtilityMacros.h index 631b0836..bde15f23 100644 --- a/common/UtilityMacros.h +++ b/common/UtilityMacros.h @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // This file contains useful macros including ERROR, WARNING, INSIST, ASSERT, etc. #ifndef included_UtilityMacros #define included_UtilityMacros diff --git a/cpu/BGK.cpp b/cpu/BGK.cpp index 90e6bf1a..5f50519a 100644 --- a/cpu/BGK.cpp +++ b/cpu/BGK.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz) { diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 526bcafa..c55ed5ff 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #define STOKES diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index 153f497b..4483453e 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include extern "C" void ScaLBL_D3Q19_Pack(int q, int *list, int start, int count, diff --git a/cpu/D3Q7.cpp b/cpu/D3Q7.cpp index 2c7b909a..e75abe32 100644 --- a/cpu/D3Q7.cpp +++ b/cpu/D3Q7.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // CPU Functions for D3Q7 Lattice Boltzmann Methods extern "C" void ScaLBL_Scalar_Pack(int *list, int count, double *sendbuf, diff --git a/cpu/Extras.cpp b/cpu/Extras.cpp index eaf9de87..3f73d594 100644 --- a/cpu/Extras.cpp +++ b/cpu/Extras.cpp @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Basic cuda functions callable from C/C++ code #include #include diff --git a/cpu/Greyscale.cpp b/cpu/Greyscale.cpp index 3fb9ff33..779b1e45 100644 --- a/cpu/Greyscale.cpp +++ b/cpu/Greyscale.cpp @@ -1,18 +1,3 @@ -/* - Copyright 2020 Equinor ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include extern "C" void diff --git a/cpu/Ion.cpp b/cpu/Ion.cpp index d77af87b..f1e8e689 100644 --- a/cpu/Ion.cpp +++ b/cpu/Ion.cpp @@ -1,6 +1,242 @@ #include #include +/***** pH equilibrium ******/ +extern "C" void ScaLBL_D3Q7_AAodd_pH_ionization(int *neighborList, double *dist, + double *Den, double *ElectricField, double *Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + int n; + double Ex, Ey, Ez; //electrical field + double ux, uy, uz; + double uEPx, uEPy, uEPz; //electrochemical induced velocity + double Ca, Cb; + double A0, A1, A2, A3, A4, A5, A6; + double B0, B1, B2, B3, B4, B5, B6; + double f0, f1, f2, f3, f4, f5, f6; + int nr1, nr2, nr3, nr4, nr5, nr6; + double rhoe, tmp; + // double factor = Di / (Vt *Vt* ionizationEnergy); + + for (n = start; n < finish; n++) { + + //Load data + //Ci = Den[n]; + Ex = ElectricField[n + 0 * Np]; + Ey = ElectricField[n + 1 * Np]; + Ez = ElectricField[n + 2 * Np]; + + ux = Velocity[n + 0 * Np]; + uy = Velocity[n + 1 * Np]; + uz = Velocity[n + 2 * Np]; + + uEPx = Di / Vt * Ex; + uEPy = Di / Vt * Ey; + uEPz = Di / Vt * Ez; + + // q=0 + // q=1 + nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist) + // q=2 + nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist) + // q=3 + nr3 = neighborList[n + 2 * Np]; // neighbor 4 + // q=4 + nr4 = neighborList[n + 3 * Np]; // neighbor 3 + // q=5 + nr5 = neighborList[n + 4 * Np]; + // q=6 + nr6 = neighborList[n + 5 * Np]; + + A0 = dist[pH_ion*7*Np + n]; + A1 = dist[pH_ion*7*Np + nr1]; // reading the A1 data into register Aq + A2 = dist[pH_ion*7*Np + nr2]; // reading the A2 data into register Aq + A3 = dist[pH_ion*7*Np + nr3]; + A4 = dist[pH_ion*7*Np + nr4]; + A5 = dist[pH_ion*7*Np + nr5]; + A6 = dist[pH_ion*7*Np + nr6]; + + /* + B0 = dist[hydroxide*7*Np + n]; + B1 = dist[hydroxide*7*Np + nr1]; // reading the B1 data into register Bq + B2 = dist[hydroxide*7*Np + nr2]; // reading the B2 data into register Bq + B3 = dist[hydroxide*7*Np + nr3]; + B4 = dist[hydroxide*7*Np + nr4]; + B5 = dist[hydroxide*7*Np + nr5]; + B6 = dist[hydroxide*7*Np + nr6]; + */ + + // charge density + rhoe = A0 + A1 + A2 + A3 + A4 + A5 + A6; + //rhoe = Ca - Cb; + // new equilibrium + tmp = sqrt(rhoe*rhoe + 4.04e-14); + Ca = rhoe + tmp; + Cb = Ca - rhoe; + + Den[pH_ion*Np + n] = Ca - Cb; + + // proton production + A1 = 0.125 * Ca * (1.0 + 4.0 * (ux + uEPx)); + A2 = 0.125 * Ca * (1.0 - 4.0 * (ux + uEPx)); + A3 = 0.125 * Ca * (1.0 + 4.0 * (uy) + uEPy); + A4 = 0.125 * Ca * (1.0 - 4.0 * (uy) + uEPy); + A5 = 0.125 * Ca * (1.0 + 4.0 * (uz) + uEPz); + A6 = 0.125 * Ca * (1.0 - 4.0 * (uz) + uEPz); + + A0 = Ca - (A1+A2+A3+A4+A5+A6); + + // hydroxide ions created by water ionization (no net charge increase) + //Cb += (f1 + f2 + f3 + f4 + f5 + f6); + // use relative mass of hydroxide + momentum conservation + B1 = 0.125 * Cb * (1.0 + 4.0 * (ux - uEPx)); + B2 = 0.125 * Cb * (1.0 - 4.0 * (ux - uEPx)); + B3 = 0.125 * Cb * (1.0 + 4.0 * (uy - uEPy)); + B4 = 0.125 * Cb * (1.0 - 4.0 * (uy - uEPy)); + B5 = 0.125 * Cb * (1.0 + 4.0 * (uz - uEPz)); + B6 = 0.125 * Cb * (1.0 - 4.0 * (uz - uEPz)); + + B0 = Cb - (B1 + B2 + B3 + B4 + B5 + B6); + + B0 = Cb - (B1 + B2 + B3 + B4 + B5 + B6); + + f0 = A0 - B0; + f1 = A1 - B1; + f2 = A2 - B2; + f3 = A3 - B3; + f4 = A4 - B4; + f5 = A5 - B5; + f6 = A6 - B6; + + dist[pH_ion*7*Np + n] = f0; + dist[pH_ion*7*Np + nr2] = f1; + dist[pH_ion*7*Np + nr1] = f2; + dist[pH_ion*7*Np + nr4] = f3; + dist[pH_ion*7*Np + nr3] = f4; + dist[pH_ion*7*Np + nr6] = f5; + dist[pH_ion*7*Np + nr5] = f6; + /* + dist[pH_ion*7*Np + n] = f0; + dist[pH_ion*7*Np + nr1] = f1; + dist[pH_ion*7*Np + nr2] = f2; + dist[pH_ion*7*Np + nr3] = f3; + dist[pH_ion*7*Np + nr4] = f4; + dist[pH_ion*7*Np + nr5] = f5; + dist[pH_ion*7*Np + nr6] = f6; + */ + + + } +} + +extern "C" void ScaLBL_D3Q7_AAeven_pH_ionization( double *dist, + double *Den, double *ElectricField, double * Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + + int n; + double Ex, Ey, Ez; //electrical field + double ux, uy, uz; + double uEPx, uEPy, uEPz; //electrochemical induced velocity + double Ca, Cb; + double A0, A1, A2, A3, A4, A5, A6; + double B0, B1, B2, B3, B4, B5, B6; + double f0, f1, f2, f3, f4, f5, f6; + double rhoe, tmp; + // double factor = Di / (Vt *Vt* ionizationEnergy); + + for (n = start; n < finish; n++) { + //Load data + //Ci = Den[n]; + Ex = ElectricField[n + 0 * Np]; + Ey = ElectricField[n + 1 * Np]; + Ez = ElectricField[n + 2 * Np]; + + ux = Velocity[n + 0 * Np]; + uy = Velocity[n + 1 * Np]; + uz = Velocity[n + 2 * Np]; + + uEPx = Di / Vt * Ex; + uEPy = Di / Vt * Ey; + uEPz = Di / Vt * Ez; + + A0 = dist[pH_ion*7*Np + n]; + A1 = dist[pH_ion*7*Np +2 * Np + n]; + A2 = dist[pH_ion*7*Np +1 * Np + n]; + A3 = dist[pH_ion*7*Np +4 * Np + n]; + A4 = dist[pH_ion*7*Np +3 * Np + n]; + A5 = dist[pH_ion*7*Np +6 * Np + n]; + A6 = dist[pH_ion*7*Np +5 * Np + n]; + + // charge density + rhoe = A0 + A1 + A2 + A3 + A4 + A5 + A6; + //rhoe = Ca - Cb; + // new equilibrium + tmp = sqrt(rhoe*rhoe + 4.04e-14); + Ca = rhoe + tmp; + Cb = Ca - rhoe; + if (Ca < 0.0) printf("Error in hydronium concentration, %f (charge density = %f) \n", Ca, rhoe); + if (Cb < 0.0) printf("Error in hydroxide concentration, %f \n", Cb); + + Den[pH_ion*Np + n] = Ca - Cb; + + // proton production + A1 = 0.125 * Ca * (1.0 + 4.0 * (ux + uEPx)); + A2 = 0.125 * Ca * (1.0 - 4.0 * (ux + uEPx)); + A3 = 0.125 * Ca * (1.0 + 4.0 * (uy) + uEPy); + A4 = 0.125 * Ca * (1.0 - 4.0 * (uy) + uEPy); + A5 = 0.125 * Ca * (1.0 + 4.0 * (uz) + uEPz); + A6 = 0.125 * Ca * (1.0 - 4.0 * (uz) + uEPz); + + A0 = Ca - (A1+A2+A3+A4+A5+A6); + + // hydroxide ions created by water ionization (no net charge increase) + //Cb += (f1 + f2 + f3 + f4 + f5 + f6); + // use relative mass of hydroxide + momentum conservation + B1 = 0.125 * Cb * (1.0 + 4.0 * (ux - uEPx)); + B2 = 0.125 * Cb * (1.0 - 4.0 * (ux - uEPx)); + B3 = 0.125 * Cb * (1.0 + 4.0 * (uy - uEPy)); + B4 = 0.125 * Cb * (1.0 - 4.0 * (uy - uEPy)); + B5 = 0.125 * Cb * (1.0 + 4.0 * (uz - uEPz)); + B6 = 0.125 * Cb * (1.0 - 4.0 * (uz - uEPz)); + + B0 = Cb - (B1 + B2 + B3 + B4 + B5 + B6); + + f0 = A0 - B0; + f1 = A1 - B1; + f2 = A2 - B2; + f3 = A3 - B3; + f4 = A4 - B4; + f5 = A5 - B5; + f6 = A6 - B6; + + if (Ez > 0.0 && n == start){ + printf("Ca = %.5g, Cb = %.5g \n", Ca, Cb); + printf(" charge density = %.5g \n", rhoe); + printf(" Ez = %.5g, A5 = %.5g, A6 = %.5g \n", Ez, f5, f6); + } + + dist[pH_ion*7*Np + n] = f0; + dist[pH_ion*7*Np +1 * Np + n] = f1; + dist[pH_ion*7*Np +2 * Np + n] = f2; + dist[pH_ion*7*Np +3 * Np + n] = f3; + dist[pH_ion*7*Np +4 * Np + n] = f4; + dist[pH_ion*7*Np +5 * Np + n] = f5; + dist[pH_ion*7*Np +6 * Np + n] = f6; + /* + dist[pH_ion*7*Np +2 * Np + n] = f1; + dist[pH_ion*7*Np +1 * Np + n] = f2; + dist[pH_ion*7*Np +4 * Np + n] = f3; + dist[pH_ion*7*Np +3 * Np + n] = f4; + dist[pH_ion*7*Np +6 * Np + n] = f5; + dist[pH_ion*7*Np +5 * Np + n] = f6; + */ + + } +} +/**** end of pH equlibrium model ********/ + + extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef, double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, int memLinks, int Nx, int Ny, int Nz, int Np){ @@ -38,59 +274,59 @@ extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo( double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut, int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count, const int N, const int Nx, const int Ny, const int Nz) { - //.................................................................................... - // Unack distribution from the recv buffer - // Distribution q matche Cqx, Cqy, Cqz - // swap rule means that the distributions in recvbuf are OPPOSITE of q - // dist may be even or odd distributions stored by stream layout - //.................................................................................... - int n, idx, label, nqm, npm, i, j, k; - double distanceLocal;//, distanceNonlocal; - double psiLocal, psiNonlocal, membranePotential; - double ap,aq; // coefficient + //.................................................................................... + // Unack distribution from the recv buffer + // Distribution q matche Cqx, Cqy, Cqz + // swap rule means that the distributions in recvbuf are OPPOSITE of q + // dist may be even or odd distributions stored by stream layout + //.................................................................................... + int n, idx, label, nqm, npm, i, j, k; + double distanceLocal;//, distanceNonlocal; + double psiLocal, psiNonlocal, membranePotential; + double ap,aq; // coefficient - - for (idx = 0; idx < count; idx++) { - n = d3q7_recvlist[idx]; - label = d3q7_linkList[idx]; - ap = 1.0; // regular streaming rule - aq = 1.0; - if (label > 0 && !(n < 0)){ - nqm = Map[n]; - distanceLocal = Distance[nqm]; - psiLocal = Psi[nqm]; + + for (idx = 0; idx < count; idx++) { + n = d3q7_recvlist[idx]; + label = d3q7_linkList[idx]; + ap = 1.0; // regular streaming rule + aq = 1.0; + if (label > 0 && !(n < 0)){ + nqm = Map[n]; + distanceLocal = Distance[nqm]; + psiLocal = Psi[nqm]; - // Get the 3-D indices from the send process - k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j; - // Streaming link the non-local distribution - i -= Cqx; j -= Cqy; k -= Cqz; - npm = k*Nx*Ny + j*Nx + i; - //distanceNonlocal = Distance[npm]; - psiNonlocal = Psi[npm]; + // Get the 3-D indices from the send process + k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j; + // Streaming link the non-local distribution + i -= Cqx; j -= Cqy; k -= Cqz; + npm = k*Nx*Ny + j*Nx + i; + //distanceNonlocal = Distance[npm]; + psiNonlocal = Psi[npm]; - membranePotential = psiLocal - psiNonlocal; - aq = MassFractionIn; - ap = MassFractionOut; + membranePotential = psiLocal - psiNonlocal; + aq = MassFractionIn; + ap = MassFractionOut; - /* link is inside membrane */ - if (distanceLocal > 0.0){ - if (membranePotential < Threshold*(-1.0)){ - ap = MassFractionIn; - aq = MassFractionOut; - } - else { - ap = ThresholdMassFractionIn; - aq = ThresholdMassFractionOut; - } - } - else if (membranePotential > Threshold){ - aq = ThresholdMassFractionIn; - ap = ThresholdMassFractionOut; - } - } - coef[2*idx]=aq; - coef[2*idx+1]=ap; - } + /* link is inside membrane */ + if (distanceLocal > 0.0){ + if (membranePotential < Threshold*(-1.0)){ + ap = MassFractionIn; + aq = MassFractionOut; + } + else { + ap = ThresholdMassFractionIn; + aq = ThresholdMassFractionOut; + } + } + else if (membranePotential > Threshold){ + aq = ThresholdMassFractionIn; + ap = ThresholdMassFractionOut; + } + } + coef[2*idx]=aq; + coef[2*idx+1]=ap; + } } @@ -350,7 +586,6 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion_v0( double Ex, Ey, Ez; //electrical field double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double f0, f1, f2, f3, f4, f5, f6; - //double X,Y,Z, factor_x, factor_y, factor_z; for (n = start; n < finish; n++) { @@ -503,47 +738,32 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, Den[n] = Ci; - /* use logistic function to prevent negative distributions*/ - //X = 4.0 * (ux + uEPx); - //Y = 4.0 * (uy + uEPy); - //Z = 4.0 * (uz + uEPz); - //factor_x = X / sqrt(1 + X*X); - //factor_y = Y / sqrt(1 + Y*Y); - //factor_z = Z / sqrt(1 + Z*Z); - // q=0 dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; // q = 1 dist[nr2] = f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); - // f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); - // q=2 dist[nr1] = f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); - // f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); // q = 3 dist[nr4] = f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); - // f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y ); // q = 4 dist[nr3] = f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); - // f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); // q = 5 dist[nr6] = f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); - // f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); // q = 6 dist[nr5] = f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); - // f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); } } @@ -559,7 +779,6 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion( double Ex, Ey, Ez; //electrical field double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double f0, f1, f2, f3, f4, f5, f6; - //double X,Y,Z, factor_x, factor_y, factor_z; for (n = start; n < finish; n++) { @@ -599,14 +818,6 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion( FluxElectrical[n + 2 * Np] = uEPz * Ci; Den[n] = Ci; - - /* use logistic function to prevent negative distributions*/ - //X = 4.0 * (ux + uEPx); - //Y = 4.0 * (uy + uEPy); - //Z = 4.0 * (uz + uEPz); - //factor_x = X / sqrt(1 + X*X); - //factor_y = Y / sqrt(1 + Y*Y); - //factor_z = Z / sqrt(1 + Z*Z); // q=0 dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; @@ -614,32 +825,26 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion( // q = 1 dist[1 * Np + n] = f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); - // f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); // q=2 dist[2 * Np + n] = f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); - // f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); // q = 3 dist[3 * Np + n] = f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); - // f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y); // q = 4 dist[4 * Np + n] = f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); - // f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); // q = 5 dist[5 * Np + n] = f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); - // f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); // q = 6 dist[6 * Np + n] = f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); - // f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); } } diff --git a/cpu/MRT.cpp b/cpu/MRT.cpp index 70a3ee37..e47b7cec 100644 --- a/cpu/MRT.cpp +++ b/cpu/MRT.cpp @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - extern "C" void INITIALIZE(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz) { int n, N; diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp index 97ee20e7..93571445 100644 --- a/cpu/Poisson.cpp +++ b/cpu/Poisson.cpp @@ -1,4 +1,5 @@ #include +#include extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList, int *Map, @@ -101,11 +102,17 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, int n; double psi; //electric potential double Ex, Ey, Ez; //electric field - double rho_e; //local charge density + double rho_e, rho_p; //local charge density double f0, f1, f2, f3, f4, f5, f6; int nr1, nr2, nr3, nr4, nr5, nr6; double rlx = 1.0 / tau; int idx; + + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = electron_charge / (kb * T); // 1 / thermal voltage; unit [Vy] for (n = start; n < finish; n++) { @@ -115,6 +122,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; idx = Map[n]; psi = Psi[idx]; + + /* Compute H30+ OH- charge density from Poisson Boltzmann statistics */ + rho_p = 1.04e-7 * (exp(psi*Vt) - exp((-1.0)*psi*Vt)); + rho_e += rho_p; // q=0 f0 = dist[n]; @@ -182,10 +193,16 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, int n; double psi; //electric potential double Ex, Ey, Ez; //electric field - double rho_e; //local charge density + double rho_e, rho_p; //local charge density double f0, f1, f2, f3, f4, f5, f6; double rlx = 1.0 / tau; int idx; + + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = electron_charge / (kb * T); // 1 / thermal voltage; unit [Vy] for (n = start; n < finish; n++) { @@ -195,6 +212,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; idx = Map[n]; psi = Psi[idx]; + + /* Compute H30+ OH- charge density from Poisson Boltzmann statistics */ + rho_p = 1.04e-7 * (exp(psi*Vt) - exp((-1.0)*psi*Vt)); + rho_e += rho_p; f0 = dist[n]; f1 = dist[2 * Np + n]; @@ -630,13 +651,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_ElectricPotential( Psi[idx] = psi - 0.5*rho_e; } } - - - extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, - double tau, double epsilon_LB, bool UseSlippingVelBC, + double tau, double Vt, double Cp, double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np) { int n; double psi; //electric potential @@ -813,8 +831,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double *Error, double tau, + double Vt, double Cp, double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np) { + int n; double psi; //electric potential double Ex, Ey, Ez; //electric field @@ -833,7 +853,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, //Load data //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral //and thus the net space charge density is zero. - rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; + //rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; + rho_e = Den_charge[n] / epsilon_LB; f0 = dist[n]; f1 = dist[2 * Np + n]; @@ -917,19 +938,476 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, } } -extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { +/** **/ +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Grotthus(int *neighborList, int *Map, + double *dist, double *Den_charge, + double *Psi, double *ElectricField, + double tau, double Vt, double Cp, + double epsilon_LB, bool UseSlippingVelBC, + int start, int finish, int Np) { + int n; + double psi, psit; //electric potential + double Ex, Ey, Ez; //electric field + double rho_i, rho_p, rho_e; //local charge density + double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, + f16, f17, f18; + int nr1, nr2, nr3, nr4, nr5, nr6, nr7, nr8, nr9, nr10, nr11, nr12, nr13, + nr14, nr15, nr16, nr17, nr18; + double sum_q; + double rlx = 1.0 / tau; + int idx; + + double W0 = 0.5; + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + double F,G,Fprime; + double factor = 1.0 / epsilon_LB; + double inVt = 1.0 / Vt; + double expsum, expdiff, term, xv; + + /* exponential series coefficients */ + double a3 = 0.3333333333333333; + double a4 = 0.25; //0.08333333333333333; + double a5 = 0.2; // 0.01666666666666667; + double a6 = 0.1666666666666667;//0.002777777777777778; + double a7 = 0.1428571428571428; //0.0003968253968253968; + double a8 = 0.125; //4.96031746031746e-05; + double a9 = 0.1111111111111111; //5.511463844797179e-06; + double a10 = 0.1; //5.511463844797178e-07; + double a11 = 0.09090909090909091; //5.010421677088344e-08; + double a12 = 0.08333333333333333; //4.17535139757362e-09; + double a13 = 0.07692307692307693; + + for (n = start; n < finish; n++) { + + //Load data + //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral + //and thus the net space charge density is zero. + rho_i = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n]; + + // q=0 + f0 = dist[n]; + // q=1 + nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist) + f1 = dist[nr1]; // reading the f1 data into register fq + + nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist) + f2 = dist[nr2]; // reading the f2 data into register fq + + // q=3 + nr3 = neighborList[n + 2 * Np]; // neighbor 4 + f3 = dist[nr3]; + + // q = 4 + nr4 = neighborList[n + 3 * Np]; // neighbor 3 + f4 = dist[nr4]; + + // q=5 + nr5 = neighborList[n + 4 * Np]; + f5 = dist[nr5]; + + // q = 6 + nr6 = neighborList[n + 5 * Np]; + f6 = dist[nr6]; + + // q=7 + nr7 = neighborList[n + 6 * Np]; + f7 = dist[nr7]; + + // q = 8 + nr8 = neighborList[n + 7 * Np]; + f8 = dist[nr8]; + + // q=9 + nr9 = neighborList[n + 8 * Np]; + f9 = dist[nr9]; + + // q = 10 + nr10 = neighborList[n + 9 * Np]; + f10 = dist[nr10]; + + // q=11 + nr11 = neighborList[n + 10 * Np]; + f11 = dist[nr11]; + + // q=12 + nr12 = neighborList[n + 11 * Np]; + f12 = dist[nr12]; + + // q=13 + nr13 = neighborList[n + 12 * Np]; + f13 = dist[nr13]; + + // q=14 + nr14 = neighborList[n + 13 * Np]; + f14 = dist[nr14]; + + // q=15 + nr15 = neighborList[n + 14 * Np]; + f15 = dist[nr15]; + + // q=16 + nr16 = neighborList[n + 15 * Np]; + f16 = dist[nr16]; + + // q=17 + //fq = dist[18*Np+n]; + nr17 = neighborList[n + 16 * Np]; + f17 = dist[nr17]; + + // q=18 + nr18 = neighborList[n + 17 * Np]; + f18 = dist[nr18]; + + + Ex = (f1 - f2 + 0.5*(f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14))*4.0; //NOTE the unit of electric field here is V/lu + Ey = (f3 - f4 + 0.5*(f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18))*4.0; + Ez = (f5 - f6 + 0.5*(f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18))*4.0; + ElectricField[n + 0 * Np] = Ex; + ElectricField[n + 1 * Np] = Ey; + ElectricField[n + 2 * Np] = Ez; + + sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18; + G = 8.0* sum_q + rho_i*factor; + + /* Use Poisson-Boltzmann for fast proton transport */ + psit = 4.0*f0; + // rho_p = Cp * (exp(psi*inVt) - exp(-psi*inVt)); + // rho_e = rho_i + rho_p; + + /* use semi-implicit scheme */ + //Wt = W0 + Cp*inVt*factor*(1.0 + 0.16666666666666667*(psit*inVt)*(psit*inVt) + 0.00833333333333333*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)); + + + for (int s=0; s<10; s++){ + /* approximate the exponential with Taylor series */ + expsum = 2.0; + xv = (psit*inVt); + expdiff = 2.0*xv; + term = xv*xv; + expsum += term; + term *= a3*xv; + expdiff += term; + term *= a4*xv; + expsum += term; + term *= a5*xv; + expdiff += term; + term *= a6*xv; + expsum += term; + term *= a7*xv; + expdiff += term; + term *= a8*xv; + expsum += term; + term *= a9*xv; + expdiff += term; + term *= a10*xv; + expsum += term; + term *= a11*xv; + expdiff += term; + term *= a12*xv; + expsum += term; + term *= a13*xv; + expdiff += term; + + /* Compare to analytical */ + double truesum = exp(xv) + exp(-1.0*xv); + double truediff = exp(xv) - exp(-1.0*xv); + + expdiff = truediff; + expsum = truesum; + + /* Newton iteration */ + F = Cp*factor*expdiff - 8.0*W0*psit + G; + Fprime = Cp*factor*inVt*expsum - 8.0*W0; + + psit -= (F / Fprime); + /* Newton iteration is successful if F=0 */ + + } + + + /* 1/ 5040 = 0.0001984126984126984 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + /* 1/ 362880 = 2.755731922398589e-06 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + /* 1/ 39916800 = 2.505210838544172e-08 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + + /* compute new psi */ + psi = 2.0*f0*(1.0 - rlx) + rlx*psit; //(1.0 / Wt)*(sum_q + 0.125*rho_i); + + idx = Map[n]; + Psi[idx] = psi; + + // q = 0 + dist[n] = W0*psi; //f0 * (1.0 - rlx) - (1.0-0.5*rlx)*W0*rho_e; + + // q = 1 + dist[nr2] = W1*psi; //f1 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 2 + dist[nr1] = W1*psi; //f2 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 3 + dist[nr4] = W1*psi; //f3 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 4 + dist[nr3] = W1*psi; //f4 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 5 + dist[nr6] = W1*psi; //f5 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 6 + dist[nr5] = W1*psi; //f6 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + //........................................................................ + + // q = 7 + dist[nr8] = W2*psi; //f7 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 8 + dist[nr7] = W2*psi; //f8 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 9 + dist[nr10] = W2*psi; //f9 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 10 + dist[nr9] = W2*psi; //f10 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 11 + dist[nr12] = W2*psi; //f11 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 12 + dist[nr11] = W2*psi; //f12 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 13 + dist[nr14] = W2*psi; //f13 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q= 14 + dist[nr13] = W2*psi; //f14 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 15 + dist[nr16] = W2*psi; //f15 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 16 + dist[nr15] = W2*psi; //f16 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 17 + dist[nr18] = W2*psi; //f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + // q = 18 + dist[nr17] = W2*psi; //f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Grotthus(int *Map, double *dist, + double *Den_charge, double *Psi, double *ElectricField, double *Error, + double tau, double Vt, double Cp, + double epsilon_LB, bool UseSlippingVelBC, + int start, int finish, int Np) { + int n; + double psi, psit; //electric potential + double Ex, Ey, Ez; //electric field + double rho_e, rho_i, rho_p; //local charge density + double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, + f16, f17, f18; + double error,sum_q; + double rlx = 1.0 / tau; + int idx; + double W0 = 0.5; double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nread, nr5; + + double F,G,Fprime; + double factor = 1.0 / epsilon_LB; + double inVt = 1.0 / Vt; + double expsum, expdiff, term, xv; - for (int idx = 0; idx < count; idx++) { - int n = list[idx]; + /* exponential series coefficients */ + double a3 = 0.3333333333333333; + double a4 = 0.25; //0.08333333333333333; + double a5 = 0.2; // 0.01666666666666667; + double a6 = 0.1666666666666667;//0.002777777777777778; + double a7 = 0.1428571428571428; //0.0003968253968253968; + double a8 = 0.125; //4.96031746031746e-05; + double a9 = 0.1111111111111111; //5.511463844797179e-06; + double a10 = 0.1; //5.511463844797178e-07; + double a11 = 0.09090909090909091; //5.010421677088344e-08; + double a12 = 0.08333333333333333; //4.17535139757362e-09; + double a13 = 0.07692307692307693; + + for (n = start; n < finish; n++) { + + //Load data + //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral + //and thus the net space charge density is zero. + rho_i = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n]; + + f0 = dist[n]; + f1 = dist[2 * Np + n]; + f2 = dist[1 * Np + n]; + f3 = dist[4 * Np + n]; + f4 = dist[3 * Np + n]; + f5 = dist[6 * Np + n]; + f6 = dist[5 * Np + n]; + + f7 = dist[8 * Np + n]; + f8 = dist[7 * Np + n]; + f9 = dist[10 * Np + n]; + f10 = dist[9 * Np + n]; + f11 = dist[12 * Np + n]; + f12 = dist[11 * Np + n]; + f13 = dist[14 * Np + n]; + f14 = dist[13 * Np + n]; + f15 = dist[16 * Np + n]; + f16 = dist[15 * Np + n]; + f17 = dist[18 * Np + n]; + f18 = dist[17 * Np + n]; + + /* Ex = (f1 - f2) * rlx * + 4.0; //NOTE the unit of electric field here is V/lu + Ey = (f3 - f4) * rlx * + 4.0; //factor 4.0 is D3Q7 lattice squared speed of sound + Ez = (f5 - f6) * rlx * 4.0; + */ + Ex = (f1 - f2 + 0.5*(f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14))*4.0; //NOTE the unit of electric field here is V/lu + Ey = (f3 - f4 + 0.5*(f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18))*4.0; + Ez = (f5 - f6 + 0.5*(f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18))*4.0; + ElectricField[n + 0 * Np] = Ex; + ElectricField[n + 1 * Np] = Ey; + ElectricField[n + 2 * Np] = Ez; + + sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18; + G = 8.0* sum_q + rho_i*factor; - dist[6 * Np + n] = W1*Vin; - dist[12 * Np + n] = W2*Vin; - dist[13 * Np + n] = W2*Vin; - dist[16 * Np + n] = W2*Vin; - dist[17 * Np + n] = W2*Vin; + /* Use Poisson-Boltzmann for fast proton transport */ + psit = 4.0*f0; + // rho_p = Cp * (exp(psi*inVt) - exp(-psi*inVt)); + // rho_e = rho_i + rho_p; + + /* use semi-implicit scheme */ + //Wt = W0 + Cp*inVt*factor*(1.0 + 0.16666666666666667*(psit*inVt)*(psit*inVt) + 0.00833333333333333*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)); + + for (int s=0; s<10; s++){ + /* approximate the exponential with Taylor series */ + expsum = 2.0; + xv = (psit*inVt); + expdiff = 2.0*xv; + term = xv*xv; + expsum += term; + term *= a3*xv; + expdiff += term; + term *= a4*xv; + expsum += term; + term *= a5*xv; + expdiff += term; + term *= a6*xv; + expsum += term; + term *= a7*xv; + expdiff += term; + term *= a8*xv; + expsum += term; + term *= a9*xv; + expdiff += term; + term *= a10*xv; + expsum += term; + term *= a11*xv; + expdiff += term; + term *= a12*xv; + expsum += term; + term *= a13*xv; + expdiff += term; + + /* Compare to analytical */ + double truesum = exp(xv) + exp(-1.0*xv); + double truediff = exp(xv) - exp(-1.0*xv); + + expdiff = truediff; + expsum = truesum; + + /* iteration */ + F = Cp*factor*expdiff - 8.0*W0*psit + G; + Fprime = Cp*factor*inVt*expsum - 8.0*W0; + + psit -= (F / Fprime); + /* Newton iteration is successful if F=0 */ + + } + + + //if (fabs(expsum - truesum) > 1e-8) printf("Error in sum (psi = %0.5g, Vt =%0.5g): approx = %0.5g, true value = %0.5g \n", psit, Vt, expsum, truesum); + //if (fabs(expdiff - truediff) > 1e-8) printf("Error in diff: approx = %0.5g, true value = %0.5g \n", expdiff, truediff); + + /* 1/ 5040 = 0.0001984126984126984 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + /* 1/ 362880 = 2.755731922398589e-06 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + /* 1/ 39916800 = 2.505210838544172e-08 *(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt)*(psit*inVt) */ + + /* compute new psi */ + psi = 2.0*f0*(1.0 - rlx) + rlx*psit; //(1.0 / Wt)*(sum_q + 0.125*rho_i); + //error = 8.0*(sum_q - f0) + rho_i*factor; + error = Cp*factor*expdiff - 8.0*f0 + G; + Error[n] = error; + + + if (error > 1e-3){ + printf(" Newton's method error (site=%i) = %0.5g \n",n,F); + } + + + idx = Map[n]; + Psi[idx] = psi; + + // q = 0 + dist[n] = W0*psi;// + + // q = 1 + dist[1 * Np + n] = W1*psi;//f1 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 2 + dist[2 * Np + n] = W1*psi;//f2 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 3 + dist[3 * Np + n] = W1*psi;//f3 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 4 + dist[4 * Np + n] = W1*psi;//f4 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 5 + dist[5 * Np + n] = W1*psi;//f5 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + // q = 6 + dist[6 * Np + n] = W1*psi;//f6 * (1.0 - rlx) +W1* (rlx * psi) - (1.0-0.5*rlx)*0.05555555555555555*rho_e; + + dist[7 * Np + n] = W2*psi;//f7 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[8 * Np + n] = W2*psi;//f8* (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[9 * Np + n] = W2*psi;//f9 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[10 * Np + n] = W2*psi;//f10 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[11 * Np + n] = W2*psi;//f11 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[12 * Np + n] = W2*psi;//f12 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[13 * Np + n] = W2*psi;//f13 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[14 * Np + n] = W2*psi;//f14 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[15 * Np + n] = W2*psi;//f15 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[16 * Np + n] = W2*psi;//f16 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[17 * Np + n] = W2*psi;//f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + dist[18 * Np + n] = W2*psi;//f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; + + //........................................................................ + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { + //double W0 = 0.5; + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int n;//nread, nr5; + + double psi = Vin; + + for (int idx = 0; idx < count; idx++) { + n = list[idx]; + + dist[6 * Np + n] = W1*psi; + dist[12 * Np + n] = W2*psi; + dist[13 * Np + n] = W2*psi; + dist[16 * Np + n] = W2*psi; + dist[17 * Np + n] = W2*psi; } } @@ -937,18 +1415,21 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np) { - + //double W0 = 0.5; double W1 = 1.0/24.0; double W2 = 1.0/48.0; + + double psi = Vout; for (int idx = 0; idx < count; idx++) { int n = list[idx]; - dist[5 * Np + n] = W1*Vout; - dist[11 * Np + n] = W2*Vout; - dist[14 * Np + n] = W2*Vout; - dist[15 * Np + n] = W2*Vout; - dist[18 * Np + n] = W2*Vout; + + dist[5 * Np + n] = W1*psi; + dist[11 * Np + n] = W2*psi; + dist[14 * Np + n] = W2*psi; + dist[15 * Np + n] = W2*psi; + dist[18 * Np + n] = W2*psi; } } @@ -959,23 +1440,24 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int Np) { double W1 = 1.0/24.0; double W2 = 1.0/48.0; + int nr5, nr11, nr14, nr15, nr18; + + double psi = Vin; for (int idx = 0; idx < count; idx++) { - int n = list[idx]; - - // Unknown distributions - nr5 = d_neighborList[n + 4 * Np]; - nr11 = d_neighborList[n + 10 * Np]; - nr15 = d_neighborList[n + 14 * Np]; - nr14 = d_neighborList[n + 13 * Np]; - nr18 = d_neighborList[n + 17 * Np]; - - dist[nr5] = W1*Vin; - dist[nr11] = W2*Vin; - dist[nr15] = W2*Vin; - dist[nr14] = W2*Vin; - dist[nr18] = W2*Vin; + int n = list[idx]; + nr5 = d_neighborList[n + 4 * Np]; + nr11 = d_neighborList[n + 10 * Np]; + nr14 = d_neighborList[n + 13 * Np]; + nr15 = d_neighborList[n + 14 * Np]; + nr18 = d_neighborList[n + 17 * Np]; + + dist[nr5] = W1*psi; + dist[nr11] = W2*psi; + dist[nr14] = W2*psi; + dist[nr15] = W2*psi; + dist[nr18] = W2*psi; } } @@ -985,21 +1467,22 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, i double W2 = 1.0/48.0; int nr6, nr12, nr13, nr16, nr17; + double psi = Vout; + for (int idx = 0; idx < count; idx++) { - - int n = list[idx]; - // unknown distributions - nr6 = d_neighborList[n + 5 * Np]; - nr12 = d_neighborList[n + 11 * Np]; - nr16 = d_neighborList[n + 15 * Np]; - nr17 = d_neighborList[n + 16 * Np]; - nr13 = d_neighborList[n + 12 * Np]; - - dist[nr6] = W1*Vout; - dist[nr12] = W2*Vout; - dist[nr16] = W2*Vout; - dist[nr17] = W2*Vout; - dist[nr13] = W2*Vout; + int n = list[idx]; + + nr6 = d_neighborList[n + 5 * Np]; + nr12 = d_neighborList[n + 11 * Np]; + nr13 = d_neighborList[n + 12 * Np]; + nr16 = d_neighborList[n + 15 * Np]; + nr17 = d_neighborList[n + 16 * Np]; + + dist[nr6] = W1*psi; + dist[nr12] = W2*psi; + dist[nr13] = W2*psi; + dist[nr16] = W2*psi; + dist[nr17] = W2*psi; } } diff --git a/cpu/dfh.cpp b/cpu/dfh.cpp index 24c1b7b2..9fef0075 100644 --- a/cpu/dfh.cpp +++ b/cpu/dfh.cpp @@ -1,35 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include diff --git a/cpu/thermal.cpp b/cpu/thermal.cpp index 5dae709d..29e7d2e3 100644 --- a/cpu/thermal.cpp +++ b/cpu/thermal.cpp @@ -1,32 +1,2 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // cpu implementation for thermal lattice boltzmann methods // copyright James McClure, 2014 diff --git a/cuda/Extras.cu b/cuda/Extras.cu index e4ad23c0..8aeedc87 100644 --- a/cuda/Extras.cu +++ b/cuda/Extras.cu @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Basic cuda functions callable from C/C++ code #include #include diff --git a/cuda/Ion.cu b/cuda/Ion.cu index 559ef7ec..c22a1a2a 100644 --- a/cuda/Ion.cu +++ b/cuda/Ion.cu @@ -5,6 +5,218 @@ #define NBLOCKS 1024 #define NTHREADS 512 + +/***** pH equilibrium ******/ +__global__ void dvc_ScaLBL_D3Q7_AAodd_pH_ionization(int *neighborList, double *dist, + double *Den, double *ElectricField, double *Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + int n; + double Ex, Ey, Ez; //electrical field + double ux, uy, uz; + double uEPx, uEPy, uEPz; //electrochemical induced velocity + double Ca, Cb; + double A0, A1, A2, A3, A4, A5, A6; + double B0, B1, B2, B3, B4, B5, B6; + double f0, f1, f2, f3, f4, f5, f6; + int nr1, nr2, nr3, nr4, nr5, nr6; + double rhoe, tmp; + + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s 10Np => odd part of dist) + // q=2 + nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist) + // q=3 + nr3 = neighborList[n + 2 * Np]; // neighbor 4 + // q=4 + nr4 = neighborList[n + 3 * Np]; // neighbor 3 + // q=5 + nr5 = neighborList[n + 4 * Np]; + // q=6 + nr6 = neighborList[n + 5 * Np]; + + A0 = dist[pH_ion*7*Np + n]; + A1 = dist[pH_ion*7*Np + nr1]; // reading the A1 data into register Aq + A2 = dist[pH_ion*7*Np + nr2]; // reading the A2 data into register Aq + A3 = dist[pH_ion*7*Np + nr3]; + A4 = dist[pH_ion*7*Np + nr4]; + A5 = dist[pH_ion*7*Np + nr5]; + A6 = dist[pH_ion*7*Np + nr6]; + + // charge density + rhoe = A0 + A1 + A2 + A3 + A4 + A5 + A6; + //rhoe = Ca - Cb; + // new equilibrium + tmp = sqrt(rhoe*rhoe + 4.04e-14); + Ca = rhoe + tmp; + Cb = Ca - rhoe; + + Den[pH_ion*Np + n] = Ca - Cb; + + // proton production + A1 = 0.125 * Ca * (1.0 + 4.0 * (ux + uEPx)); + A2 = 0.125 * Ca * (1.0 - 4.0 * (ux + uEPx)); + A3 = 0.125 * Ca * (1.0 + 4.0 * (uy) + uEPy); + A4 = 0.125 * Ca * (1.0 - 4.0 * (uy) + uEPy); + A5 = 0.125 * Ca * (1.0 + 4.0 * (uz) + uEPz); + A6 = 0.125 * Ca * (1.0 - 4.0 * (uz) + uEPz); + + A0 = Ca - (A1+A2+A3+A4+A5+A6); + + // hydroxide ions created by water ionization (no net charge increase) + //Cb += (f1 + f2 + f3 + f4 + f5 + f6); + // use relative mass of hydroxide + momentum conservation + B1 = 0.125 * Cb * (1.0 + 4.0 * (ux - uEPx)); + B2 = 0.125 * Cb * (1.0 - 4.0 * (ux - uEPx)); + B3 = 0.125 * Cb * (1.0 + 4.0 * (uy - uEPy)); + B4 = 0.125 * Cb * (1.0 - 4.0 * (uy - uEPy)); + B5 = 0.125 * Cb * (1.0 + 4.0 * (uz - uEPz)); + B6 = 0.125 * Cb * (1.0 - 4.0 * (uz - uEPz)); + + B0 = Cb - (B1 + B2 + B3 + B4 + B5 + B6); + + B0 = Cb - (B1 + B2 + B3 + B4 + B5 + B6); + + f0 = A0 - B0; + f1 = A1 - B1; + f2 = A2 - B2; + f3 = A3 - B3; + f4 = A4 - B4; + f5 = A5 - B5; + f6 = A6 - B6; + + dist[pH_ion*7*Np + n] = f0; + dist[pH_ion*7*Np + nr2] = f1; + dist[pH_ion*7*Np + nr1] = f2; + dist[pH_ion*7*Np + nr4] = f3; + dist[pH_ion*7*Np + nr3] = f4; + dist[pH_ion*7*Np + nr6] = f5; + dist[pH_ion*7*Np + nr5] = f6; + + } + } +} + +__global__ void dvc_ScaLBL_D3Q7_AAeven_pH_ionization( double *dist, + double *Den, double *ElectricField, double * Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + + int n; + double Ex, Ey, Ez; //electrical field + double ux, uy, uz; + double uEPx, uEPy, uEPz; //electrochemical induced velocity + double Ca, Cb; + double A0, A1, A2, A3, A4, A5, A6; + double B0, B1, B2, B3, B4, B5, B6; + double f0, f1, f2, f3, f4, f5, f6; + double rhoe, tmp; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s0) + CD_tmp; } } } @@ -980,3 +1161,33 @@ extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, } } +extern "C" void ScaLBL_D3Q7_AAodd_pH_ionization(int *neighborList, double *dist, + double *Den, double *ElectricField, double *Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + + dvc_ScaLBL_D3Q7_AAodd_pH_ionization<<>>(neighborList,dist,Den,ElectricField, + Velocity,Di,Vt,pH_ion,start,finish,Np); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_AAodd_pH_ionization: %s \n",cudaGetErrorString(err)); + } + +} + +extern "C" void ScaLBL_D3Q7_AAeven_pH_ionization( double *dist, + double *Den, double *ElectricField, double * Velocity, + double Di, double Vt, + int pH_ion, int start, int finish, int Np) { + + dvc_ScaLBL_D3Q7_AAeven_pH_ionization<<>>(dist,Den,ElectricField, + Velocity,Di,Vt,pH_ion,start,finish,Np); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in dvc_ScaLBL_D3Q7_AAeven_pH_ionization: %s \n",cudaGetErrorString(err)); + } + +} + diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index c86043a5..00c9e3da 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -522,7 +522,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, //Load data //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral //and thus the net space charge density is zero. - rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; + //rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; + rho_e = Den_charge[n] / epsilon_LB; f0 = dist[n]; f1 = dist[2 * Np + n]; @@ -681,16 +682,15 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborLi double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nr5, nr11, nr14, nr15, nr18; - + int nr5, nr11, nr14, nr15, nr18; + int idx = blockIdx.x*blockDim.x + threadIdx.x; if (idx < count){ int n = list[idx]; - - // Unknown distributions - nr5 = d_neighborList[n + 4 * Np]; + // Unknown distributions + nr5 = d_neighborList[n + 4 * Np]; nr11 = d_neighborList[n + 10 * Np]; nr15 = d_neighborList[n + 14 * Np]; nr14 = d_neighborList[n + 13 * Np]; @@ -769,8 +769,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, i extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, - double *Psi, double *ElectricField, - double tau, double epsilon_LB, bool UseSlippingVelBC, + double *Psi, double *ElectricField, + double tau, double Vt, double Cp, + double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np) { //cudaProfilerStart(); dvc_ScaLBL_D3Q19_AAodd_Poisson<<>>(neighborList, Map, @@ -783,8 +784,8 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, } extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, - double *Den_charge, double *Psi, - double *ElectricField, double *Error, double tau, + double *Den_charge, double *Psi, double *ElectricField, double *Error, + double tau, double Vt, double Cp, double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np) { diff --git a/hip/BGK.hip b/hip/BGK.hip index aef19bfe..765ca536 100644 --- a/hip/BGK.hip +++ b/hip/BGK.hip @@ -1,19 +1,4 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include "hip/hip_runtime.h" diff --git a/hip/Color.hip b/hip/Color.hip index 0199043f..d009cb7f 100644 --- a/hip/Color.hip +++ b/hip/Color.hip @@ -1,19 +1,4 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include #include "hip/hip_runtime.h" diff --git a/hip/CudaExtras.hip b/hip/CudaExtras.hip index 39ac1d26..7be9c1ac 100644 --- a/hip/CudaExtras.hip +++ b/hip/CudaExtras.hip @@ -1,19 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Basic hip functions callable from C/C++ code #include "hip/hip_runtime.h" diff --git a/hip/D3Q19.hip b/hip/D3Q19.hip index 86ab36bd..3d1150da 100644 --- a/hip/D3Q19.hip +++ b/hip/D3Q19.hip @@ -1,19 +1,4 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include "hip/hip_runtime.h" #include "hip/hip_cooperative_groups.h" diff --git a/hip/D3Q7.hip b/hip/D3Q7.hip index 0dc3dc82..8fb3328b 100644 --- a/hip/D3Q7.hip +++ b/hip/D3Q7.hip @@ -1,19 +1,4 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // GPU Functions for D3Q7 Lattice Boltzmann Methods #include "hip/hip_runtime.h" diff --git a/hip/Greyscale.hip b/hip/Greyscale.hip index 17d3c584..c1b77235 100644 --- a/hip/Greyscale.hip +++ b/hip/Greyscale.hip @@ -1,18 +1,3 @@ -/* - Copyright 2020 Equinor ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include "hip/hip_runtime.h" diff --git a/hip/MRT.hip b/hip/MRT.hip index e4ef3bfb..771b6e30 100644 --- a/hip/MRT.hip +++ b/hip/MRT.hip @@ -1,20 +1,6 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ //************************************************************************* -// CUDA kernels for single-phase ScaLBL_D3Q19_MRT code +// HIP kernels for single-phase ScaLBL_D3Q19_MRT code // James McClure //************************************************************************* #include "hip/hip_runtime.h" diff --git a/hip/dfh.hip b/hip/dfh.hip index d09578b0..ad0a51ef 100644 --- a/hip/dfh.hip +++ b/hip/dfh.hip @@ -1,35 +1,4 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - Copyright Equnior ASA - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include #include #include "hip/hip_runtime.h" diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 099b034e..1b42d0a2 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -7,7 +7,7 @@ #include "common/ReadMicroCT.h" ScaLBL_IonModel::ScaLBL_IonModel(int RANK, int NP, const Utilities::MPI &COMM) - : rank(RANK), nprocs(NP), timestep(0), timestepMax(0), time_conv(0), kb(0), + : rank(RANK), nprocs(NP), timestep(0), timestepGlobal(0), timestepMax(0), time_conv(0), kb(0), electron_charge(0), T(0), Vt(0), k2_inv(0), h(0), tolerance(0), number_ion_species(0), Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), fluidVelx_dummy(0), fluidVely_dummy(0), @@ -30,74 +30,7 @@ ScaLBL_IonModel::~ScaLBL_IonModel() { void ScaLBL_IonModel::ReadParams(string filename, vector &num_iter) { - // read the input database - db = std::make_shared(filename); - domain_db = db->getDatabase("Domain"); - ion_db = db->getDatabase("Ions"); - - // Universal constant - kb = 1.38e-23; //Boltzmann constant;unit [J/K] - electron_charge = 1.6e-19; //electron charge;unit [C] - - //---------------------- Default model parameters --------------------------// - T = 300.0; //temperature; unit [K] - Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] - k2_inv = 4.0; //speed of sound for D3Q7 lattice - h = 1.0; //resolution; unit: um/lu - tolerance = 1.0e-8; - number_ion_species = 1; - tau.push_back(1.0); - IonDiffusivity.push_back( - 1.0e-9); //user-input diffusivity has physical unit [m^2/sec] - IonValence.push_back(1); //algebraic valence charge - IonConcentration.push_back( - 1.0e-3); //user-input ion concentration has physical unit [mol/m^3] - //tau.push_back(0.5+k2_inv*time_conv/(h*1.0e-6)/(h*1.0e-6)*IonDiffusivity[0]); - time_conv.push_back((tau[0] - 0.5) / k2_inv * (h * h * 1.0e-12) / - IonDiffusivity[0]); - fluidVelx_dummy = 0.0; //for debugging, unit [m/sec] - fluidVely_dummy = 0.0; //for debugging, unit [m/sec] - fluidVelz_dummy = 0.0; //for debugging, unit [m/sec] - Ex_dummy = 0.0; //for debugging, unit [V/m] - Ey_dummy = 0.0; //for debugging, unit [V/m] - Ez_dummy = 0.0; //for debugging, unit [V/m] - - sprintf(LocalRankString, "%05d", rank); - sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); - //--------------------------------------------------------------------------// - - BoundaryConditionSolid = 0; - - // Read domain parameters - if (domain_db->keyExists("voxel_length")) { //default unit: um/lu - h = domain_db->getScalar("voxel_length"); - } - - // LB-Ion Model parameters - //if (ion_db->keyExists( "timestepMax" )){ - // timestepMax = ion_db->getScalar( "timestepMax" ); - //} - if (ion_db->keyExists("tolerance")) { - tolerance = ion_db->getScalar("tolerance"); - } - if (ion_db->keyExists("temperature")) { - T = ion_db->getScalar("temperature"); - //re-calculate thermal voltage - Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] - } - if (ion_db->keyExists("FluidVelDummy")) { - fluidVelx_dummy = ion_db->getVector("FluidVelDummy")[0]; - fluidVely_dummy = ion_db->getVector("FluidVelDummy")[1]; - fluidVelz_dummy = ion_db->getVector("FluidVelDummy")[2]; - } - if (ion_db->keyExists("ElectricFieldDummy")) { - Ex_dummy = ion_db->getVector("ElectricFieldDummy")[0]; - Ey_dummy = ion_db->getVector("ElectricFieldDummy")[1]; - Ez_dummy = ion_db->getVector("ElectricFieldDummy")[2]; - } - if (ion_db->keyExists("number_ion_species")) { - number_ion_species = ion_db->getScalar("number_ion_species"); - } + //------ Load number of iteration from multiphysics controller ------// if (num_iter.size() != number_ion_species) { ERROR("Error: number_ion_species and num_iter_Ion_List (from " @@ -105,196 +38,9 @@ void ScaLBL_IonModel::ReadParams(string filename, vector &num_iter) { } else { timestepMax.assign(num_iter.begin(), num_iter.end()); } - //-------------------------------------------------------------------// - if (ion_db->keyExists("tauList")) { - tau.clear(); - tau = ion_db->getVector("tauList"); - vector Di = ion_db->getVector( - "IonDiffusivityList"); //temp storing ion diffusivity in physical unit - if (tau.size() != number_ion_species || - Di.size() != number_ion_species) { - ERROR("Error: number_ion_species, tauList and IonDiffusivityList " - "must be of the same length! \n"); - } else { - time_conv.clear(); - for (size_t i = 0; i < tau.size(); i++) { - time_conv.push_back((tau[i] - 0.5) / k2_inv * - (h * h * 1.0e-12) / Di[i]); - } - } - } - //read ion related list - //NOTE: ion diffusivity has INPUT unit: [m^2/sec] - // it must be converted to LB unit: [lu^2/lt] - if (ion_db->keyExists("IonDiffusivityList")) { - IonDiffusivity.clear(); - IonDiffusivity = ion_db->getVector("IonDiffusivityList"); - if (IonDiffusivity.size() != number_ion_species) { - ERROR("Error: number_ion_species and IonDiffusivityList must be " - "the same length! \n"); - } else { - for (size_t i = 0; i < IonDiffusivity.size(); i++) { - IonDiffusivity[i] = - IonDiffusivity[i] * time_conv[i] / - (h * h * 1.0e-12); //LB diffusivity has unit [lu^2/lt] - } - } - } else { - for (size_t i = 0; i < IonDiffusivity.size(); i++) { - //convert ion diffusivity in physical unit to LB unit - IonDiffusivity[i] = - IonDiffusivity[i] * time_conv[i] / - (h * h * 1.0e-12); //LB diffusivity has unit [lu^2/lt] - } - } - // read time relaxation time list - //read ion algebric valence list - if (ion_db->keyExists("IonValenceList")) { - IonValence.clear(); - IonValence = ion_db->getVector("IonValenceList"); - if (IonValence.size() != number_ion_species) { - ERROR("Error: number_ion_species and IonValenceList must be the " - "same length! \n"); - } - } - //read initial ion concentration list; INPUT unit [mol/m^3] - //it must be converted to LB unit [mol/lu^3] - if (ion_db->keyExists("IonConcentrationList")) { - IonConcentration.clear(); - IonConcentration = ion_db->getVector("IonConcentrationList"); - if (IonConcentration.size() != number_ion_species) { - ERROR("Error: number_ion_species and IonConcentrationList must be " - "the same length! \n"); - } else { - for (size_t i = 0; i < IonConcentration.size(); i++) { - IonConcentration[i] = - IonConcentration[i] * - (h * h * h * - 1.0e-18); //LB ion concentration has unit [mol/lu^3] - } - } - } else { - for (size_t i = 0; i < IonConcentration.size(); i++) { - IonConcentration[i] = - IonConcentration[i] * - (h * h * h * - 1.0e-18); //LB ion concentration has unit [mol/lu^3] - } - } - - //Read solid boundary condition specific to Ion model - BoundaryConditionSolid = 0; - if (ion_db->keyExists("BC_Solid")) { - BoundaryConditionSolid = ion_db->getScalar("BC_Solid"); - } - // Read boundary condition for ion transport - // BC = 0: zero-flux bounce-back BC - // BC = 1: fixed ion concentration; unit=[mol/m^3] - // BC = 2: fixed ion flux (inward flux); unit=[mol/m^2/sec] - BoundaryConditionInlet.push_back(0); - BoundaryConditionOutlet.push_back(0); - //Inlet - if (ion_db->keyExists("BC_InletList")) { - BoundaryConditionInlet = ion_db->getVector("BC_InletList"); - if (BoundaryConditionInlet.size() != number_ion_species) { - ERROR("Error: number_ion_species and BC_InletList must be of the " - "same length! \n"); - } - unsigned short int BC_inlet_min = *min_element( - BoundaryConditionInlet.begin(), BoundaryConditionInlet.end()); - unsigned short int BC_inlet_max = *max_element( - BoundaryConditionInlet.begin(), BoundaryConditionInlet.end()); - if (BC_inlet_min == 0 && BC_inlet_max > 0) { - ERROR("Error: BC_InletList: mix of periodic, ion concentration and " - "flux BC is not supported! \n"); - } - if (BC_inlet_min > 0) { - //read in inlet values Cin - if (ion_db->keyExists("InletValueList")) { - Cin = ion_db->getVector("InletValueList"); - if (Cin.size() != number_ion_species) { - ERROR("Error: number_ion_species and InletValueList must " - "be the same length! \n"); - } - } else { - ERROR("Error: Non-periodic BCs are specified but " - "InletValueList cannot be found! \n"); - } - for (size_t i = 0; i < BoundaryConditionInlet.size(); i++) { - switch (BoundaryConditionInlet[i]) { - case 1: //fixed boundary ion concentration [mol/m^3] - Cin[i] = - Cin[i] * - (h * h * h * - 1.0e-18); //LB ion concentration has unit [mol/lu^3] - break; - case 21: //fixed boundary ion flux [mol/m^2/sec] - Cin[i] = Cin[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - case 22: //fixed boundary ion flux [mol/m^2/sec] - Cin[i] = Cin[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - case 23: //fixed boundary ion flux [mol/m^2/sec] - Cin[i] = Cin[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - } - } - } - } - //Outlet - if (ion_db->keyExists("BC_OutletList")) { - BoundaryConditionOutlet = ion_db->getVector("BC_OutletList"); - if (BoundaryConditionOutlet.size() != number_ion_species) { - ERROR("Error: number_ion_species and BC_OutletList must be of the " - "same length! \n"); - } - unsigned short int BC_outlet_min = *min_element( - BoundaryConditionOutlet.begin(), BoundaryConditionOutlet.end()); - unsigned short int BC_outlet_max = *max_element( - BoundaryConditionOutlet.begin(), BoundaryConditionOutlet.end()); - if (BC_outlet_min == 0 && BC_outlet_max > 0) { - ERROR("Error: BC_OutletList: mix of periodic, ion concentration " - "and flux BC is not supported! \n"); - } - if (BC_outlet_min > 0) { - //read in outlet values Cout - if (ion_db->keyExists("OutletValueList")) { - Cout = ion_db->getVector("OutletValueList"); - if (Cout.size() != number_ion_species) { - ERROR("Error: number_ion_species and OutletValueList must " - "be the same length! \n"); - } - } else { - ERROR("Error: Non-periodic BCs are specified but " - "OutletValueList cannot be found! \n"); - } - for (size_t i = 0; i < BoundaryConditionOutlet.size(); i++) { - switch (BoundaryConditionOutlet[i]) { - case 1: //fixed boundary ion concentration [mol/m^3] - Cout[i] = - Cout[i] * - (h * h * h * - 1.0e-18); //LB ion concentration has unit [mol/lu^3] - break; - case 21: //fixed boundary ion flux [mol/m^2/sec] - Cout[i] = Cout[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - case 22: //fixed boundary ion flux [mol/m^2/sec] - Cout[i] = Cout[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - case 23: //fixed boundary ion flux [mol/m^2/sec] - Cout[i] = Cout[i] * (h * h * 1.0e-12) * - time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] - break; - } - } - } - } + + ReadParams(filename); + } void ScaLBL_IonModel::ReadParams(string filename) { @@ -310,6 +56,9 @@ void ScaLBL_IonModel::ReadParams(string filename) { // Universal constant kb = 1.38e-23; //Boltzmann constant;unit [J/K] electron_charge = 1.6e-19; //electron charge;unit [C] + + use_Grotthus = false; + pH_ion = -1; //---------------------- Default model parameters --------------------------// T = 300.0; //temperature; unit [K] @@ -335,7 +84,7 @@ void ScaLBL_IonModel::ReadParams(string filename) { Ez_dummy = 0.0; //for debugging, unit [V/m] sprintf(LocalRankString, "%05d", rank); - sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); + sprintf(LocalRestartFile, "%s%s", "IonModel.", LocalRankString); //--------------------------------------------------------------------------// // Read domain parameters @@ -387,8 +136,13 @@ void ScaLBL_IonModel::ReadParams(string filename) { for (size_t i = 0; i < tau.size(); i++) { time_conv.push_back((tau[i] - 0.5) / k2_inv * (h * h * 1.0e-12) / Di[i]); + if (rank==0) printf(" tauList[%zu] = %.5g \n",i,tau[i]); + if (rank==0) printf(" Di[%zu] = %.5g \n",i,Di[i]); + if (rank==0) printf(" time_conv[%zu] = %.5g \n",i,time_conv[i]); + } } + } //read ion related list //NOTE: ion diffusivity has INPUT unit: [m^2/sec] @@ -397,17 +151,19 @@ void ScaLBL_IonModel::ReadParams(string filename) { IonDiffusivity.clear(); IonDiffusivity = ion_db->getVector("IonDiffusivityList"); if (IonDiffusivity.size() != number_ion_species) { - ERROR("Error: number_ion_species and IonDiffusivityList must be " - "the same length! \n"); + ERROR("Error: number_ion_species and IonDiffusivityList must be " + "the same length! \n"); } else { - for (size_t i = 0; i < IonDiffusivity.size(); i++) { - IonDiffusivity[i] = - IonDiffusivity[i] * time_conv[i] / - (h * h * 1.0e-12); //LB diffusivity has unit [lu^2/lt] - } + for (size_t i = 0; i < IonDiffusivity.size(); i++) { + IonDiffusivity[i] = + IonDiffusivity[i] * time_conv[i] / + (h * h * 1.0e-12); //LB diffusivity has unit [lu^2/lt] + if (rank==0) printf(" IonDiffusivityList[%zu] = %.5g [lu^2/lt] \n",i,IonDiffusivity[i]); + } } + } else { - for (size_t i = 0; i < IonDiffusivity.size(); i++) { + for (size_t i = 0; i < IonDiffusivity.size(); i++) { //convert ion diffusivity in physical unit to LB unit IonDiffusivity[i] = IonDiffusivity[i] * time_conv[i] / @@ -423,6 +179,28 @@ void ScaLBL_IonModel::ReadParams(string filename) { ERROR("Error: number_ion_species and IonValenceList must be the " "same length! \n"); } + for (size_t i = 0; i < IonValence.size(); i++) { + if (rank==0) printf(" IonValenceList[%zu] = %i \n",i,IonValence[i]); + } + } + + if (ion_db->keyExists("WaterIonList")) { + use_Grotthus = true; + vector GrotthusList = ion_db->getVector("WaterIonList"); + IonizationEnergy = ion_db->getWithDefault("IonizationEnergy", 20.24); // in eV + if (GrotthusList.size() != 1 ) { + ERROR("Error: water ionization model only supports one pH ion \n"); + } + else { + pH_ion = GrotthusList[0]; + if (rank==0) printf( " Grotthus mechanism enabled. " + "pH ion is %zu, \n",pH_ion); + } + /* check that user specifies consistent valence charge */ + if (IonValence[pH_ion] != 1){ + ERROR("Error: pH ion must have valence charge of +1 \n"); + } + } //read initial ion concentration list; INPUT unit [mol/m^3] //it must be converted to LB unit [mol/lu^3] @@ -432,22 +210,28 @@ void ScaLBL_IonModel::ReadParams(string filename) { if (IonConcentration.size() != number_ion_species) { ERROR("Error: number_ion_species and IonConcentrationList must be " "the same length! \n"); - } else { - for (size_t i = 0; i < IonConcentration.size(); i++) { - IonConcentration[i] = - IonConcentration[i] * - (h * h * h * - 1.0e-18); //LB ion concentration has unit [mol/lu^3] - } + } + else { + for (size_t i = 0; i < IonConcentration.size(); i++) { + IonConcentration[i] = + IonConcentration[i] * + (h * h * h * + 1.0e-18); //LB ion concentration has unit [mol/lu^3] + if (rank==0) printf(" IonConcentrationList[%zu] = %.5g [mol/lu^3] \n",i,IonConcentration[i]); + } } - } else { + } + else { + if (rank==0) printf(" IonConcentrationList not specified in input database. Setting default values \n"); for (size_t i = 0; i < IonConcentration.size(); i++) { IonConcentration[i] = IonConcentration[i] * (h * h * h * 1.0e-18); //LB ion concentration has unit [mol/lu^3] + if (rank==0) printf(" IonConcentrationList[%zu] = %.5g [mol/lu^3] \n",i,IonConcentration[i]); } } + if (USE_MEMBRANE){ membrane_db = db->getDatabase("Membrane"); @@ -546,7 +330,6 @@ void ScaLBL_IonModel::ReadParams(string filename) { } } } - } //Read solid boundary condition specific to Ion model BoundaryConditionSolid = 0; @@ -557,8 +340,10 @@ void ScaLBL_IonModel::ReadParams(string filename) { // BC = 0: zero-flux bounce-back BC // BC = 1: fixed ion concentration; unit=[mol/m^3] // BC = 2: fixed ion flux (inward flux); unit=[mol/m^2/sec] - BoundaryConditionInlet.push_back(0); - BoundaryConditionOutlet.push_back(0); + for (size_t i = 0; i < number_ion_species; i++) { + BoundaryConditionInlet.push_back(0); + BoundaryConditionOutlet.push_back(0); + } //Inlet if (ion_db->keyExists("BC_InletList")) { BoundaryConditionInlet = ion_db->getVector("BC_InletList"); @@ -574,7 +359,7 @@ void ScaLBL_IonModel::ReadParams(string filename) { ERROR("Error: BC_InletList: mix of periodic, ion concentration and " "flux BC is not supported! \n"); } - if (BC_inlet_min > 0) { + if (BC_inlet_min > 1) { //read in inlet values Cin if (ion_db->keyExists("InletValueList")) { Cin = ion_db->getVector("InletValueList"); @@ -588,21 +373,21 @@ void ScaLBL_IonModel::ReadParams(string filename) { } for (size_t i = 0; i < BoundaryConditionInlet.size(); i++) { switch (BoundaryConditionInlet[i]) { - case 1: //fixed boundary ion concentration [mol/m^3] + case 2: //fixed boundary ion concentration [mol/m^3] Cin[i] = Cin[i] * (h * h * h * 1.0e-18); //LB ion concentration has unit [mol/lu^3] break; - case 21: //fixed boundary ion flux [mol/m^2/sec] + case 3: //fixed boundary ion flux [mol/m^2/sec] Cin[i] = Cin[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; - case 22: //fixed boundary ion flux [mol/m^2/sec] + case 4: //fixed boundary ion flux [mol/m^2/sec] Cin[i] = Cin[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; - case 23: //fixed boundary ion flux [mol/m^2/sec] + case 5: //fixed boundary ion flux [mol/m^2/sec] Cin[i] = Cin[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; @@ -625,7 +410,7 @@ void ScaLBL_IonModel::ReadParams(string filename) { ERROR("Error: BC_OutletList: mix of periodic, ion concentration " "and flux BC is not supported! \n"); } - if (BC_outlet_min > 0) { + if (BC_outlet_min > 1) { //read in outlet values Cout if (ion_db->keyExists("OutletValueList")) { Cout = ion_db->getVector("OutletValueList"); @@ -639,21 +424,21 @@ void ScaLBL_IonModel::ReadParams(string filename) { } for (size_t i = 0; i < BoundaryConditionOutlet.size(); i++) { switch (BoundaryConditionOutlet[i]) { - case 1: //fixed boundary ion concentration [mol/m^3] + case 2: //fixed boundary ion concentration [mol/m^3] Cout[i] = Cout[i] * (h * h * h * 1.0e-18); //LB ion concentration has unit [mol/lu^3] break; - case 21: //fixed boundary ion flux [mol/m^2/sec] + case 3: //fixed boundary ion flux [mol/m^2/sec] Cout[i] = Cout[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; - case 22: //fixed boundary ion flux [mol/m^2/sec] + case 4: //fixed boundary ion flux [mol/m^2/sec] Cout[i] = Cout[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; - case 23: //fixed boundary ion flux [mol/m^2/sec] + case 5: //fixed boundary ion flux [mol/m^2/sec] Cout[i] = Cout[i] * (h * h * 1.0e-12) * time_conv[i]; //LB ion flux has unit [mol/lu^2/lt] break; @@ -661,6 +446,17 @@ void ScaLBL_IonModel::ReadParams(string filename) { } } } + if (ion_db->keyExists("BC_frequency")) { + BC_frequency = ion_db->getVector("BC_frequency"); + } + if (ion_db->keyExists("BC_amplitude")) { + BC_amplitude = ion_db->getVector("BC_amplitude"); + if (BC_amplitude.size() != number_ion_species || + BC_amplitude.size() != number_ion_species) { + ERROR("Error: number_ion_species and boundary condition specification " + "must be of the same length! \n"); + } + } } void ScaLBL_IonModel::SetDomain() { @@ -779,7 +575,7 @@ void ScaLBL_IonModel::ReadInput() { sprintf(LocalRankString, "%05d", Dm->rank()); sprintf(LocalRankFilename, "%s%s", "ID.", LocalRankString); - sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); + sprintf(LocalRestartFile, "%s%s", "IonModel.", LocalRankString); if (domain_db->keyExists("Filename")) { auto Filename = domain_db->getScalar("Filename"); @@ -993,8 +789,8 @@ void ScaLBL_IonModel::Create() { if (rank == 0) printf("LB Ion Solver: Allocating distributions \n"); //......................device distributions................................. - int dist_mem_size = Np * sizeof(double); - int neighborSize = 18 * (Np * sizeof(int)); + size_t dist_mem_size = Np * sizeof(double); + size_t neighborSize = 18 * (Np * sizeof(int)); //........................................................................... ScaLBL_AllocateDeviceMemory((void **)&NeighborList, neighborSize); ScaLBL_AllocateDeviceMemory((void **)&dvcMap, sizeof(int) * Np); @@ -1157,12 +953,11 @@ void ScaLBL_IonModel::Initialize() { Ci_host = new double[number_ion_species * Np]; ifstream File(LocalRestartFile, ios::binary); - int idx; - double value,sum; + double value; // Read the distributions for (size_t ic = 0; ic < number_ion_species; ic++){ for (int n = 0; n < Np; n++) { - sum = 0.0; + double sum = 0.0; for (int q = 0; q < 7; q++) { File.read((char *)&value, sizeof(value)); cDist[ic * 7 * Np + q * Np + n] = value; @@ -1220,25 +1015,31 @@ void ScaLBL_IonModel::Initialize() { i + 1); break; case 1: + if (rank == 0) + printf( + "LB Ion Solver: outlet boundary for Ion %zu is bounce-back \n", + i + 1); + break; + case 2: if (rank == 0) printf("LB Ion Solver: inlet boundary for Ion %zu is " "concentration = %.5g [mol/m^3] \n", i + 1, Cin[i] / (h * h * h * 1.0e-18)); break; - case 21: + case 3: if (rank == 0) printf("LB Ion Solver: inlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive flux only. \n", i + 1, Cin[i] / (h * h * 1.0e-12) / time_conv[i]); break; - case 22: + case 4: if (rank == 0) printf( "LB Ion Solver: inlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive + advective flux. \n", i + 1, Cin[i] / (h * h * 1.0e-12) / time_conv[i]); break; - case 23: + case 5: if (rank == 0) printf("LB Ion Solver: inlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive + advective + " @@ -1255,24 +1056,30 @@ void ScaLBL_IonModel::Initialize() { break; case 1: if (rank == 0) - printf("LB Ion Solver: outlet boundary for Ion %zu is " - "concentration = %.5g [mol/m^3] \n", - i + 1, Cout[i] / (h * h * h * 1.0e-18)); - break; - case 21: - if (rank == 0) + printf( + "LB Ion Solver: outlet boundary for Ion %zu is bounce-back \n", + i + 1); + break; + case 2: + if (rank == 0) + printf("LB Ion Solver: outlet boundary for Ion %zu is " + "concentration = %.5g [mol/m^3] \n", + i + 1, Cout[i] / (h * h * h * 1.0e-18)); + break; + case 3: + if (rank == 0) printf("LB Ion Solver: outlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive flux only. \n", i + 1, Cout[i] / (h * h * 1.0e-12) / time_conv[i]); break; - case 22: + case 4: if (rank == 0) printf( "LB Ion Solver: outlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive + advective flux. \n", i + 1, Cout[i] / (h * h * 1.0e-12) / time_conv[i]); break; - case 23: + case 5: if (rank == 0) printf("LB Ion Solver: outlet boundary for Ion %zu is (inward) " "flux = %.5g [mol/m^2/sec]; Diffusive + advective + " @@ -1332,23 +1139,23 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); //--------------------------------------- Set boundary conditions -------------------------------------// - if (BoundaryConditionInlet[ic] > 0) { + if (BoundaryConditionInlet[ic] > 1) { switch (BoundaryConditionInlet[ic]) { - case 1: + case 2: ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); break; - case 21: + case 3: ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 22: + case 4: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 23: + case 5: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], &ElectricField[2 * Np], @@ -1356,23 +1163,23 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { break; } } - if (BoundaryConditionOutlet[ic] > 0) { + if (BoundaryConditionOutlet[ic] > 1) { switch (BoundaryConditionOutlet[ic]) { - case 1: + case 2: ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); break; - case 21: + case 3: ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 22: + case 4: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 23: + case 5: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], &ElectricField[2 * Np], @@ -1416,23 +1223,23 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); //--------------------------------------- Set boundary conditions -------------------------------------// - if (BoundaryConditionInlet[ic] > 0) { + if (BoundaryConditionInlet[ic] > 1) { switch (BoundaryConditionInlet[ic]) { - case 1: + case 2: ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); break; - case 21: + case 3: ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 22: + case 4: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 23: + case 5: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], &Velocity[2 * Np], &ElectricField[2 * Np], @@ -1440,23 +1247,23 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { break; } } - if (BoundaryConditionOutlet[ic] > 0) { + if (BoundaryConditionOutlet[ic] > 1) { switch (BoundaryConditionOutlet[ic]) { - case 1: + case 2: ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); break; - case 21: + case 3: ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 22: + case 4: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], timestep); break; - case 23: + case 5: ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], &Velocity[2 * Np], &ElectricField[2 * Np], @@ -1496,8 +1303,8 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) { ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, - ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, + 0, ScaLBL_Comm->LastExterior(), Np); } //************************************************************************/ if (rank == 0) @@ -1537,40 +1344,132 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl //double starttime,stoptime,cputime; //ScaLBL_Comm->Barrier(); comm.barrier(); //auto t1 = std::chrono::system_clock::now(); - + for (size_t ic = 0; ic < number_ion_species; ic++) { + + bool BounceBack = false; + if (BoundaryConditionInlet[ic] > 0) + BounceBack = true; /* set the mass transfer coefficients for the membrane */ - IonMembrane->AssignCoefficients(dvcMap, Psi, ThresholdVoltage[ic],MassFractionIn[ic], - MassFractionOut[ic],ThresholdMassFractionIn[ic],ThresholdMassFractionOut[ic]); + if (USE_MEMBRANE) + IonMembrane->AssignCoefficients(dvcMap, Psi, ThresholdVoltage[ic],MassFractionIn[ic], + MassFractionOut[ic],ThresholdMassFractionIn[ic],ThresholdMassFractionOut[ic]); timestep = 0; while (timestep < timestepMax[ic]) { //************************************************************************/ // *************ODD TIMESTEP*************// timestep++; + timestepGlobal++; //LB-Ion collison IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL - ScaLBL_D3Q7_AAodd_Ion( - IonMembrane->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); + if ( ic == pH_ion ){ + ScaLBL_D3Q7_AAodd_pH_ionization(IonMembrane->NeighborList, fq, + Ci, ElectricField, Velocity, + rlx[pH_ion], Vt, pH_ion, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + } + else { + ScaLBL_D3Q7_AAodd_Ion( + IonMembrane->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); + } - IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE + IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7], BounceBack); //WRITE INTO OPPOSITE + /* SET BOUNDARY CONDITIONS */ + /* + //--------------------------------------- Set boundary conditions -------------------------------------// + if ( ic != pH_ion ){ + if (BoundaryConditionInlet[ic] == 2) { + double BC_value = Cin[ic]*(1.0+BC_amplitude[ic]*sin(timestepGlobal*BC_frequency[ic])); + //printf("Setting inlet BC phase = %4.3e \n", BC_value); + IonMembrane->D3Q7_Ion_Concentration_BC_z( + NeighborList, &fq[ic * Np * 7],BC_value, timestep); - ScaLBL_D3Q7_AAodd_Ion( - IonMembrane->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 (BoundaryConditionInlet[ic] == 2) { + double BC_value = Cout[ic]*(1.0-BC_amplitude[ic]*sin(timestepGlobal*BC_frequency[ic])); + //printf("Setting outlet BC phase = %4.3e \n", BC_value); + IonMembrane->D3Q7_Ion_Concentration_BC_Z( + NeighborList, &fq[ic * Np * 7], BC_value, timestep); + } + } + else { + */ + { + if (BoundaryConditionInlet[ic] > 1) { + switch (BoundaryConditionInlet[ic]) { + case 2: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); + break; + case 3: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 4: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 5: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + if (BoundaryConditionOutlet[ic] > 1) { + switch (BoundaryConditionOutlet[ic]) { + case 2: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); + break; + case 3: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 4: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 5: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + } - IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]); + if ( ic == pH_ion ){ + ScaLBL_D3Q7_AAodd_pH_ionization(IonMembrane->NeighborList, fq, + Ci, ElectricField, Velocity, + rlx[pH_ion], Vt, pH_ion, + 0, ScaLBL_Comm->LastExterior(), Np); + } + else { + ScaLBL_D3Q7_AAodd_Ion( + IonMembrane->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 (USE_MEMBRANE) + IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]); /* if (BoundaryConditionSolid == 1) { //TODO IonSolid may also be species-dependent @@ -1581,27 +1480,114 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl */ // *************EVEN TIMESTEP*************// timestep++; + timestepGlobal++; //LB-Ion collison IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL - 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); + if ( ic == pH_ion ){ + ScaLBL_D3Q7_AAeven_pH_ionization( fq, + Ci, ElectricField, Velocity, + rlx[pH_ion], Vt, pH_ion, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + } + else { + 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); + } + IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7], BounceBack); //WRITE INTO OPPOSITE + /* + //--------------------------------------- Set boundary conditions -------------------------------------// + if ( ic != pH_ion ){ + if (BoundaryConditionInlet[ic] == 2) { + double BC_value = Cin[ic]*(1.0+BC_amplitude[ic]*sin(timestepGlobal*BC_frequency[ic])); + //printf("Setting inlet BC phase = %4.3e \n", BC_value); + IonMembrane->D3Q7_Ion_Concentration_BC_z( + NeighborList, &fq[ic * Np * 7],BC_value, timestep); - IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE + } + if (BoundaryConditionInlet[ic] == 2) { + double BC_value = Cout[ic]*(1.0-BC_amplitude[ic]*sin(timestepGlobal*BC_frequency[ic])); + //printf("Setting outlet BC phase = %4.3e \n", BC_value); + IonMembrane->D3Q7_Ion_Concentration_BC_Z( + NeighborList, &fq[ic * Np * 7], BC_value, timestep); - 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); - - IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]); + } + } + else {*/ + { + if (BoundaryConditionInlet[ic] > 1) { + switch (BoundaryConditionInlet[ic]) { + case 2: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], timestep); + break; + case 3: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 4: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 5: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z( + NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + if (BoundaryConditionOutlet[ic] > 1) { + switch (BoundaryConditionOutlet[ic]) { + case 2: + ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], timestep); + break; + case 3: + ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 4: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], timestep); + break; + case 5: + ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z( + NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic], + &Velocity[2 * Np], &ElectricField[2 * Np], + IonDiffusivity[ic], IonValence[ic], Vt, timestep); + break; + } + } + } + // End BC + + if ( ic == pH_ion ){ + ScaLBL_D3Q7_AAeven_pH_ionization( fq, + Ci, ElectricField, Velocity, + rlx[pH_ion], Vt, pH_ion, + 0, ScaLBL_Comm->LastExterior(), Np); + } + else { + 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 (USE_MEMBRANE) + IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]); ScaLBL_Comm->Barrier(); comm.barrier(); @@ -1617,17 +1603,17 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl } } - //Compute charge density for Poisson equation - for (size_t ic = 0; ic < number_ion_species; ic++) { - int Valence = IonValence[ic]; - ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, - ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, 0, - ScaLBL_Comm->LastExterior(), Np); - } + //Compute charge density for Poisson equation + for (size_t ic = 0; ic < number_ion_species; ic++) { + int Valence = IonValence[ic]; + ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, + 0, ScaLBL_Comm->LastExterior(), Np); + } - /* DoubleArray Charge(Nx,Ny,Nz); + /* DoubleArray Charge(Nx,Ny,Nz); ScaLBL_Comm->RegularLayout(Map, ChargeDensity, Charge); double charge_sum=0.0; double charge_sum_total=0.0; @@ -1662,13 +1648,76 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl //if (rank==0) printf("********************************************************\n"); } +void ScaLBL_IonModel::TestGrotthus() { + /* Set up dummy electric field */ + double * Psi; + double * ElectricField; + double * Velocity; + ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz); + ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); + + if (rank == 0) { + printf(" TestGrotthus: create dummy electric field... \n"); + } + + double ALPHA = 0.001; + double *PSI, *EF, *VEL; + EF = new double [3*Np]; + VEL = new double [3*Np]; + PSI = new double [Nx*Ny*Nz]; + + if (rank == 0) { + printf(" TestGrotthus: initialize variables (%i, %i, %i)... \n",Nx,Ny,Nz); + } + + int z = ScaLBL_Comm->kproc * (Nz-2); + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + int idx = Map(i, j, k); + int n = k * Nx * Ny + j * Nx + i; + double VALUE = (z + k - 1)*ALPHA; + Psi[n] = VALUE; + if (!(idx < 0)) { + EF[idx] = 0.0; + EF[Np + idx] = 0.0; + EF[2*Np + idx] = ALPHA; + + VEL[idx] = 0.0; + VEL[Np + idx] = 0.0; + VEL[2*Np + idx] = 0.0; + } + } + } + } + if (rank == 0) { + printf(" TestGrotthus: copy to GPU... \n"); + } + ScaLBL_CopyToDevice(Psi, PSI, sizeof(double) * Nx * Ny * Nz); + ScaLBL_CopyToDevice(ElectricField, EF, sizeof(double) * 3 * Np); + ScaLBL_CopyToDevice(Velocity, VEL, sizeof(double) * 3 * Np); + + if (rank == 0) { + printf(" TestGrotthus: run model... \n"); + } + RunMembrane(Velocity, ElectricField, Psi); + + delete(PSI); + delete(EF); + delete(VEL); + ScaLBL_FreeDeviceMemory(Psi); + ScaLBL_FreeDeviceMemory(Velocity); + ScaLBL_FreeDeviceMemory(ElectricField); +} + + void ScaLBL_IonModel::Checkpoint(){ if (rank == 0) { printf(" ION MODEL: Writing restart file! \n"); } - int idx; - double value,sum; + double value; double*cDist; cDist = new double[7 * number_ion_species * Np]; ScaLBL_CopyToHost(cDist, fq, 7 * Np * number_ion_species *sizeof(double)); diff --git a/models/IonModel.h b/models/IonModel.h index 9a6756e6..498d8ef8 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -50,10 +50,14 @@ public: void DummyFluidVelocity(); void DummyElectricField(); void Checkpoint(); + + void TestGrotthus(); + double CalIonDenConvergence(vector &ci_avg_previous); bool Restart; int timestep; + int timestepGlobal; vector timestepMax; int BoundaryConditionSolid; double h; //domain resolution, unit [um/lu] @@ -62,6 +66,10 @@ public: double tolerance; double fluidVelx_dummy, fluidVely_dummy, fluidVelz_dummy; double Ex_dummy, Ey_dummy, Ez_dummy; + + bool use_Grotthus; + size_t pH_ion; + double IonizationEnergy; size_t number_ion_species; vector BoundaryConditionInlet; @@ -79,6 +87,8 @@ public: vector Cout; //outlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec] vector tau; vector time_conv; + vector BC_frequency; + vector BC_amplitude; int Nx, Ny, Nz, N, Np; int rank, nprocx, nprocy, nprocz, nprocs; diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 1af357db..2f3e419b 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -1,5 +1,5 @@ /* - * Multi-relaxation time LBM Model + * Gauss's Law solver */ #include "models/PoissonSolver.h" #include "analysis/distance.h" @@ -36,10 +36,11 @@ ScaLBL_Poisson::~ScaLBL_Poisson() ScaLBL_FreeDeviceMemory(dvcMap); ScaLBL_FreeDeviceMemory(Psi); ScaLBL_FreeDeviceMemory(Psi_BCLabel); + ScaLBL_FreeDeviceMemory(Permittivity); ScaLBL_FreeDeviceMemory(ElectricField); ScaLBL_FreeDeviceMemory(ResidualError); ScaLBL_FreeDeviceMemory(fq); - + if ( TIMELOG ) fclose( TIMELOG ); } @@ -66,7 +67,7 @@ void ScaLBL_Poisson::ReadParams(string filename){ TestPeriodicTime = 1.0;//unit: [sec] TestPeriodicTimeConv = 0.01; //unit [sec/lt] TestPeriodicSaveInterval = 0.1; //unit [sec] - Restart = "false"; + Restart = false; // LB-Poisson Model parameters if (electric_db->keyExists( "Restart" )){ @@ -141,7 +142,7 @@ void ScaLBL_Poisson::ReadParams(string filename){ /* restart string */ sprintf(LocalRankString, "%05d", rank); - sprintf(LocalRestartFile, "%s%s", "Psi.", LocalRankString); + sprintf(LocalRestartFile, "%s%s", "PoissonSolver.", LocalRankString); if (rank==0) printf("***********************************************************************************\n"); if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance); @@ -211,7 +212,7 @@ void ScaLBL_Poisson::ReadInput(){ sprintf(LocalRankString,"%05d",Dm->rank()); sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); - sprintf(LocalRestartFile,"%s%s","Psi.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","PoissonSolver.",LocalRankString); if (domain_db->keyExists( "Filename" )){ @@ -277,6 +278,44 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid, int *poisson_sol if (NLABELS != AffinityList.size() || NLABELS != BoundaryConditionSolidList.size()){ ERROR("Error: LB-Poisson Solver: BC_SolidList, SolidLabels and SolidValues all must be of the same length! \n"); } + + if (electric_db->keyExists( "PermittivityValues" )) + { + /* assign the permittivity based on the material*/ + double *Permittivity_host; + Permittivity_host = new double[Nx*Ny*Nz]; + + double PERMITTIVITY = epsilon_LB; + auto LabelList = electric_db->getVector( "SolidLabels" ); + auto PermittivityList = electric_db->getVector( "PermittivityValues" ); + size_t NLABELS = LabelList.size(); + if (NLABELS != PermittivityList.size()){ + ERROR("Error: LB-Poisson Solver: SolidLabels and PermittivityList all must be of the same length! \n"); + } + + for (int k=0;kid[n]; + PERMITTIVITY=epsilon_LB; + // Assign the affinity from the paired list + for (unsigned int idx=0; idx < NLABELS; idx++){ + if (VALUE == LabelList[idx]){ + PERMITTIVITY=PermittivityList[idx]; + //label_count[idx] += 1.0; + idx = NLABELS; + } + } + int idx=Map(i,j,k); + if (!(idx<0)) Permittivity_host[n] = PERMITTIVITY; + } + } + } + + ScaLBL_CopyToDevice(Permittivity, Permittivity_host, sizeof(double)*Nx*Ny*Nz); + delete [] Permittivity_host; + } std::vector label_count( NLABELS, 0.0 ); std::vector label_count_global( NLABELS, 0.0 ); @@ -368,13 +407,14 @@ void ScaLBL_Poisson::Create(){ // LBM variables if (rank==0) printf ("LB-Poisson Solver: Allocating distributions \n"); //......................device distributions................................. - int dist_mem_size = Np*sizeof(double); - int neighborSize=18*(Np*sizeof(int)); + size_t dist_mem_size = Np*sizeof(double); + size_t neighborSize=18*(Np*sizeof(int)); //........................................................................... ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); //ScaLBL_AllocateDeviceMemory((void **) &dvcID, sizeof(signed char)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz); + ScaLBL_AllocateDeviceMemory((void **) &Permittivity, sizeof(double)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &Psi_BCLabel, sizeof(int)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &ResidualError, sizeof(double)*Np); @@ -443,8 +483,12 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ Vin0 = Vout0 = 1.0; //unit: [V] freqIn = freqOut = 50.0; //unit: [Hz] PhaseShift_In = PhaseShift_Out = 0.0; //unit: [radian] - Vin = 1.0; //Boundary-z (inlet) electric potential - Vout = 1.0; //Boundary-Z (outlet) electric potential + Vin = 0.0; //Boundary-z (inlet) electric potential + Vout = 0.0; //Boundary-Z (outlet) electric potential + + /* Assign permittivity value to the solid */ + signed char VALUE=0; + double AFFINITY=0.f; if (BoundaryConditionInlet==0 && BoundaryConditionOutlet==0){ @@ -453,6 +497,7 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ auto LabelList = electric_db->getVector( "InitialValueLabels" ); auto AffinityList = electric_db->getVector( "InitialValues" ); + auto PermittivityList = electric_db->getVector( "PermittivityValues" ); size_t NLABELS = LabelList.size(); if (NLABELS != AffinityList.size()){ @@ -549,26 +594,36 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ if (BoundaryConditionOutlet==2) Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,0); //initialize a linear electrical potential between inlet and outlet - double slope = (Vout-Vin)/(Nz-2); - double psi_linearized; + double slope = (Vout-Vin)/((Nz-2)*Dm->nprocz()); + double psi_linearized = Vin; for (int k=0;kid[n]>0){ - psi_init[n] = psi_linearized; - } - } - } + if (Dm->kproc() == 0){ + if (k==0 || k==1){ + psi_linearized = Vin; + } + else{ + psi_linearized = slope*(Dm->kproc()*(Nz-2) + (k-1)) + Vin; + } + } + if (Dm->kproc() == Dm->nprocz()-1){ + if (k==Nz-1 || k==Nz-2){ + psi_linearized = Vout; + } + else{ + psi_linearized = slope*(Dm->kproc()*(Nz-2) + (k-1)) + Vin; + } + } + else{ + psi_linearized = slope*(Dm->kproc()*(Nz-2) + (k-1)) + Vin; + } + for (int j=0;jid[n]>0){ + psi_init[n] = psi_linearized; + } + } + } } } else{//mixed periodic and non-periodic BCs are not supported! @@ -643,79 +698,6 @@ void ScaLBL_Poisson::Initialize(double time_conv_from_Study){ //} } -//void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){ -// -// //.......create and start timer............ -// //double starttime,stoptime,cputime; -// //comm.barrier(); -// //auto t1 = std::chrono::system_clock::now(); -// double *host_Error; -// host_Error = new double [Np]; -// -// timestep=0; -// double error = 1.0; -// while (timestep < timestepMax && error > tolerance) { -// //************************************************************************/ -// // *************ODD TIMESTEP*************// -// timestep++; -// //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential -// SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision -// ScaLBL_Comm->Barrier(); comm.barrier(); -// -// // *************EVEN TIMESTEP*************// -// timestep++; -// //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential -// SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision -// ScaLBL_Comm->Barrier(); comm.barrier(); -// //************************************************************************/ -// -// -// // Check convergence of steady-state solution -// if (timestep==2){ -// //save electric potential for convergence check -// } -// if (timestep%analysis_interval==0){ -// /* get the elecric potential */ -// ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); -// if (rank==0) printf(" ... getting Poisson solver error \n"); -// double err = 0.0; -// double max_error = 0.0; -// ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np); -// for (int idx=0; idx max_error ){ -// max_error = err; -// } -// } -// error=Dm->Comm.maxReduce(max_error); -// -// /* compute the eletric field */ -// //ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np); -// -// } -// } -// if(WriteLog==true){ -// getConvergenceLog(timestep,error); -// } -// -// //************************************************************************/ -// ////if (rank==0) printf("LB-Poission Solver: a steady-state solution is obtained\n"); -// ////if (rank==0) printf("---------------------------------------------------------------------------\n"); -// //// Compute the walltime per timestep -// //auto t2 = std::chrono::system_clock::now(); -// //double cputime = std::chrono::duration( t2 - t1 ).count() / timestep; -// //// Performance obtained from each node -// //double MLUPS = double(Np)/cputime/1000000; -// -// //if (rank==0) printf("******************* LB-Poisson Solver Statistics ********************\n"); -// //if (rank==0) printf("CPU time = %f \n", cputime); -// //if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); -// //MLUPS *= nprocs; -// //if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); -// //if (rank==0) printf("*********************************************************************\n"); -// -//} - void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){ double error = 1.0; @@ -794,44 +776,120 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times timestep=0; auto t1 = std::chrono::system_clock::now(); while (timestep < timestepMax && error > tolerance) { - //************************************************************************/ - // *************ODD TIMESTEP*************// - timestep++; - //SolveElectricPotentialAAodd(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision - ScaLBL_Comm->Barrier(); comm.barrier(); + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] + double Cp = 1.014e-7; // proton concentration - // *************EVEN TIMESTEP*************// - timestep++; - //SolveElectricPotentialAAeven(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision - ScaLBL_Comm->Barrier(); comm.barrier(); - //************************************************************************/ + timestep=0; + auto t1 = std::chrono::system_clock::now(); + while (timestep < timestepMax && error > tolerance) { + //************************************************************************/ + // *************ODD TIMESTEP*************// + // Set boundary conditions + timestep++; - // Check convergence of steady-state solution - if (timestep==2){ - //save electric potential for convergence check - } - if (timestep%analysis_interval==0){ - /* get the elecric potential */ - ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); - if (rank==0) printf(" ... getting Poisson solver error \n"); - double err = 0.0; - double max_error = 0.0; - ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np); - for (int idx=0; idx max_error ){ - max_error = err; - } - } - error=Dm->Comm.maxReduce(max_error); + //SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC, timestep);//perform collision + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - /* compute the eletric field */ - //ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np); + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } - } + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); + + // *************EVEN TIMESTEP*************// + + timestep++; + + //SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC, timestep);//perform collision + //ScaLBL_Comm->Barrier(); comm.barrier(); + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + + // Set boundary conditions + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } + + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); comm.barrier(); + //************************************************************************/ + } + // Check convergence of steady-state solution + if (timestep==2){ + //save electric potential for convergence check + } + if (timestep%analysis_interval==0){ + /* get the elecric potential */ + ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); + if (rank==0) printf(" ... getting Poisson solver error \n"); + double err = 0.0; + double max_error = 0.0; + ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np); + for (int idx=0; idx max_error ){ + max_error = err; + } + } + error=Dm->Comm.maxReduce(max_error); + if (rank==0) printf(" error = %0.5g \n",error); + } } + // SetMeanZeroVoltage(); + if (rank == 0) printf("---------------------------------------------------------------" "----\n"); @@ -840,7 +898,10 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times double cputime = std::chrono::duration(t2 - t1).count() / timestep; // Performance obtained from each node double MLUPS = double(Np) / cputime / 1000000; - + + if(WriteLog==true){ + getConvergenceLog(timestep,error); + } if (rank == 0) printf("********************************************************\n"); if (rank == 0) @@ -876,6 +937,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo double *host_Error; host_Error = new double [Np]; + + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] + double Cp = 1.014e-7; // proton concentration timestep=0; auto t1 = std::chrono::system_clock::now(); @@ -883,14 +951,80 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo //************************************************************************/ // *************ODD TIMESTEP*************// timestep++; - //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision - ScaLBL_Comm->Barrier(); comm.barrier(); + //SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC, timestep);//perform collision + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + // Set boundary conditions + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); + // *************EVEN TIMESTEP*************// timestep++; - //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision + + + //SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC, timestep);//perform collision + //ScaLBL_Comm->Barrier(); comm.barrier(); + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + + // Set boundary conditions + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); + ScaLBL_Comm->D3Q19_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, + tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); + + ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ @@ -972,6 +1106,8 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo //ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np); } } + // SetMeanZeroVoltage(); + if (rank == 0) printf("---------------------------------------------------------------" "----\n"); @@ -1000,6 +1136,43 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo } } +void ScaLBL_Poisson::SetMeanZeroVoltage(){ + /* get the elecric potential */ + ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); + + double local_mean_voltage = 0.0; + double global_mean_voltage = 0.0; + double local_count = 0.0; + double global_count = 0.0; + for (int k=1; kComm.sumReduce(local_mean_voltage); + global_count = Dm->Comm.sumReduce(local_count); + global_mean_voltage /= global_count; + // rescale the far-field electric potential + for (int k=0; k 0){ switch (BoundaryConditionInlet){ case 1: - ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; case 2: Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); - ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; } } if (BoundaryConditionOutlet > 0){ switch (BoundaryConditionOutlet){ case 1: - ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vin, timestep); break; case 2: Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); - ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vin, timestep); break; } } @@ -1147,6 +1320,14 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC, int timestep){ + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] + double Cp = 1.014e-7; // proton concentration + + if (lattice_scheme.compare("D3Q7")==0){ ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); @@ -1163,35 +1344,35 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe } else if (lattice_scheme.compare("D3Q19")==0){ ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - //ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE // Set boundary conditions - if (BoundaryConditionInlet > 0){ + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ switch (BoundaryConditionInlet){ case 1: - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; case 2: Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep); - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; } } - if (BoundaryConditionOutlet > 0){ + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ switch (BoundaryConditionOutlet){ case 1: - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep); break; case 2: Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep); - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep); break; } } - ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes. //something like: @@ -1200,6 +1381,13 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe } void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC, int timestep){ + + // Universal constant + double kb = 1.38e-23; //Boltzmann constant;unit [J/K] + double electron_charge = 1.6e-19; //electron charge;unit [C] + double T = 300.0; //temperature; unit [K] + double Vt = kb * T / electron_charge; //thermal voltage; unit [Vy] + double Cp = 1.014e-7; // proton concentration if (lattice_scheme.compare("D3Q7")==0){ ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1214,35 +1402,34 @@ void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingV } else if (lattice_scheme.compare("D3Q19")==0){ ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - // ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE // Set boundary conditions - if (BoundaryConditionInlet > 0){ + if (BoundaryConditionInlet > 0 && Dm->kproc()==0){ switch (BoundaryConditionInlet){ case 1: - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; case 2: Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep); - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep); break; } } - if (BoundaryConditionOutlet > 0){ + if (BoundaryConditionOutlet > 0 && Dm->kproc() == nprocz-1){ switch (BoundaryConditionOutlet){ case 1: - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep); break; case 2: Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep); - ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep); break; } } - ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, Vt, Cp, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); @@ -1279,7 +1466,8 @@ void ScaLBL_Poisson::DummyChargeDensity(){ for (int i=0; iRegularLayout(Map,ResidualError,ReturnValues); +} + void ScaLBL_Poisson::getElectricPotential(DoubleArray &ReturnValues){ //This function wirte out the data in a normal layout (by aggregating all decomposed domains) //ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField); @@ -1373,12 +1566,15 @@ void ScaLBL_Poisson::WriteVis( int timestep) { auto vis_db = db->getDatabase("Visualization"); auto format = vis_db->getWithDefault( "format", "hdf5" ); - DoubleArray ElectricalPotential(Nx, Ny, Nz); std::vector visData; fillHalo fillData(Dm->Comm, Dm->rank_info, {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1}, 0, 1); + DoubleArray ElectricalPotential(Nx, Ny, Nz); + + DoubleArray SolverError(Nx, Ny, Nz); + IO::initialize("",format,"false"); // Create the MeshDataStruct visData.resize(1); @@ -1389,9 +1585,23 @@ void ScaLBL_Poisson::WriteVis( int timestep) { Dm->Nz - 2, Dm->Lx, Dm->Ly, Dm->Lz); //electric potential auto ElectricPotentialVar = std::make_shared(); + auto SolverErrorVar = std::make_shared(); //-------------------------------------------------------------------------------------------------------------------- + DoubleArray Analytical(Nx, Ny, Nz); + for (int k=0; kgetWithDefault("save_electric_potential", true)) { ElectricPotentialVar->name = "ElectricPotential"; @@ -1400,6 +1610,14 @@ void ScaLBL_Poisson::WriteVis( int timestep) { ElectricPotentialVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); visData[0].vars.push_back(ElectricPotentialVar); } + + if (vis_db->getWithDefault("save_error", true)) { + SolverErrorVar->name = "SolverError"; + SolverErrorVar->type = IO::VariableType::VolumeVariable; + SolverErrorVar->dim = 1; + SolverErrorVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); + visData[0].vars.push_back(SolverErrorVar); + } //-------------------------------------------------------------------------------------------------------------------- //------------------------------------Save All Variables-------------------------------------------------------------- @@ -1410,140 +1628,16 @@ void ScaLBL_Poisson::WriteVis( int timestep) { fillData.copy(ElectricalPotential, ElectricPotentialData); } + //------------------------------------Save All Variables-------------------------------------------------------------- + if (vis_db->getWithDefault("save_error", true)) { + ASSERT(visData[0].vars[1]->name == "SolverError"); + getSolverError(SolverError); + Array &SolverErrorData = visData[0].vars[1]->data; + fillData.copy(SolverError, SolverErrorData); + } + if (vis_db->getWithDefault("write_silo", true)) IO::writeData(timestep, visData, Dm->Comm); //-------------------------------------------------------------------------------------------------------------------- } -//void ScaLBL_Poisson::SolveElectricField(){ -// ScaLBL_Comm_Regular->SendHalo(Psi); -// ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid, -// Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); -// ScaLBL_Comm_Regular->RecvHalo(Psi); -// ScaLBL_Comm->Barrier(); -// if (BoundaryCondition == 1){ -// ScaLBL_Comm->Poisson_D3Q7_BC_z(dvcMap,Psi,Vin); -// ScaLBL_Comm->Poisson_D3Q7_BC_Z(dvcMap,Psi,Vout); -// } -// ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); -// -//} - -//void ScaLBL_Poisson::getElectricPotential(){ -// -// DoubleArray PhaseField(Nx,Ny,Nz); -// ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField); -// //ScaLBL_Comm->Barrier(); comm.barrier(); -// FILE *OUTFILE; -// sprintf(LocalRankFilename,"Electric_Potential.%05i.raw",rank); -// OUTFILE = fopen(LocalRankFilename,"wb"); -// fwrite(PhaseField.data(),8,N,OUTFILE); -// fclose(OUTFILE); -//} - -//old version where Psi is of size Np -//void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid) -//{ -// size_t NLABELS=0; -// signed char VALUE=0; -// double AFFINITY=0.f; -// -// auto LabelList = electric_db->getVector( "SolidLabels" ); -// auto AffinityList = electric_db->getVector( "SolidValues" ); -// -// NLABELS=LabelList.size(); -// if (NLABELS != AffinityList.size()){ -// ERROR("Error: LB-Poisson Solver: SolidLabels and SolidValues must be the same length! \n"); -// } -// -// double label_count[NLABELS]; -// double label_count_global[NLABELS]; -// // Assign the labels -// -// for (size_t idx=0; idxid[n]; -// AFFINITY=0.f; -// // Assign the affinity from the paired list -// for (unsigned int idx=0; idx < NLABELS; idx++){ -// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]); -// if (VALUE == LabelList[idx]){ -// AFFINITY=AffinityList[idx]; -// //NOTE need to convert the user input phys unit to LB unit -// if (BoundaryConditionSolid==2){ -// //for BCS=1, i.e. Dirichlet-type, no need for unit conversion -// //TODO maybe there is a factor of gamm missing here ? -// AFFINITY = AFFINITY*(h*h*1.0e-12)/epsilon_LB; -// } -// label_count[idx] += 1.0; -// idx = NLABELS; -// //Mask->id[n] = 0; // set mask to zero since this is an immobile component -// } -// } -// poisson_solid[n] = AFFINITY; -// } -// } -// } -// -// for (size_t idx=0; idxComm.sumReduce( label_count[idx]); -// -// if (rank==0){ -// printf("LB-Poisson Solver: number of Poisson solid labels: %lu \n",NLABELS); -// for (unsigned int idx=0; idxkeyExists( "Vin" )){ -// Vin = electric_db->getScalar( "Vin" ); -// } -// if (electric_db->keyExists( "Vout" )){ -// Vout = electric_db->getScalar( "Vout" ); -// } -// } -// //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis -// double slope = (Vout-Vin)/(Nz-2); -// double psi_linearized; -// for (int k=0;k. -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #include "threadpool/atomic_helpers.h" #include diff --git a/threadpool/atomic_helpers.h b/threadpool/atomic_helpers.h index 7cafcc4c..6d36312e 100644 --- a/threadpool/atomic_helpers.h +++ b/threadpool/atomic_helpers.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Copyright © 2004 Mark Berrill. All Rights Reserved. This work is distributed with permission, // but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. #ifndef included_ThreadPoolAtomicHelpers diff --git a/threadpool/atomic_list.h b/threadpool/atomic_list.h index da632dc6..c60c2869 100644 --- a/threadpool/atomic_list.h +++ b/threadpool/atomic_list.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_AtomicModelAtomicList #define included_AtomicModelAtomicList diff --git a/threadpool/atomic_list.hpp b/threadpool/atomic_list.hpp index bc65ae35..3a4df598 100644 --- a/threadpool/atomic_list.hpp +++ b/threadpool/atomic_list.hpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ #ifndef included_AtomicList_hpp #define included_AtomicList_hpp diff --git a/threadpool/thread_pool.h b/threadpool/thread_pool.h index 1273f405..9cb5b21a 100644 --- a/threadpool/thread_pool.h +++ b/threadpool/thread_pool.h @@ -1,18 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // Copyright © 2004 Mark Berrill. All Rights Reserved. This work is distributed with permission, // but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A // PARTICULAR PURPOSE. diff --git a/threadpool/thread_pool.hpp b/threadpool/thread_pool.hpp index 34838446..394e5619 100644 --- a/threadpool/thread_pool.hpp +++ b/threadpool/thread_pool.hpp @@ -1,33 +1,3 @@ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -/* - Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University - - This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ // This file contains the template functions for the thread pool #ifndef included_ThreadPoolTmpl #define included_ThreadPoolTmpl