test only with HDF5

.github/workflows/c-cpp.yml

@@ -102,7 +102,7 @@ jobs:
           -D TIMER_DIRECTORY=$LBPM_TIMER_DIR     \
         -D USE_NETCDF=0                             \
           -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR   \
-        -D USE_SILO=1                               \
+        -D USE_SILO=0                               \
           -D HDF5_DIRECTORY=$LBPM_HDF5_DIR         \
           -D SILO_DIRECTORY=$LBPM_SILO_DIR         \
         -D USE_CUDA=0                               \

CMakeLists.txt

@@ -119,7 +119,6 @@ ADD_DISTCLEAN( analysis null_timer tests liblbpm-wia.* cpu gpu cuda hip example
 # Check for CUDA
 CHECK_ENABLE_FLAG( USE_CUDA 0 )
 CHECK_ENABLE_FLAG( USE_HIP 0 )
-CHECK_ENABLE_FLAG( USE_SYCL 0 )
 NULL_USE( CMAKE_CUDA_FLAGS )
 IF ( USE_CUDA )
     ADD_DEFINITIONS( -DUSE_CUDA )
@@ -127,8 +126,6 @@ IF ( USE_CUDA )
 ELSEIF ( USE_HIP )
     ENABLE_LANGUAGE( HIP )
     ADD_DEFINITIONS( -DUSE_HIP )
-ELSEIF ( USE_SYCL )
-    ADD_DEFINITIONS( -DUSE_SYCL )
 ENDIF()
 
 
@@ -170,8 +167,6 @@ IF ( NOT ONLY_BUILD_DOCS )
         ADD_PACKAGE_SUBDIRECTORY( cuda )
     ELSEIF ( USE_HIP )
         ADD_PACKAGE_SUBDIRECTORY( hip )
-    ELSEIF ( USE_SYCL )
-        ADD_PACKAGE_SUBDIRECTORY( sycl )
     ELSE()
         ADD_PACKAGE_SUBDIRECTORY( cpu )
     ENDIF()

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef IO_HELPERS_INC
 #define IO_HELPERS_INC
 

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 <http://www.gnu.org/licenses/>.
-*/
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #include "Mesh.h"
 #include "IO/IOHelpers.h"
 #include "common/Utilities.h"

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef MESH_INC
 #define MESH_INC
 

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 <http://www.gnu.org/licenses/>.
-*/
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #include "IO/MeshDatabase.h"
 #include "IO/IOHelpers.h"
 #include "IO/Mesh.h"

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef MeshDatabase_INC
 #define MeshDatabase_INC
 

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 <http://www.gnu.org/licenses/>.
-*/
 #include "IO/PIO.h"
 #include "common/MPI.h"
 #include "common/Utilities.h"

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_PIO
 #define included_PIO
 

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 <http://www.gnu.org/licenses/>.
-*/
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_PIO_hpp
 #define included_PIO_hpp
 

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 <http://www.gnu.org/licenses/>.
-*/
 #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 <ProfilerApp.h>

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef READER_INC
 #define READER_INC
 

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 <http://www.gnu.org/licenses/>.
-*/
 #include "IO/Writer.h"
 #include "IO/HDF5_IO.h"
 #include "IO/IOHelpers.h"

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef WRITER_INC
 #define WRITER_INC
 

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 <http://www.gnu.org/licenses/>.
-*/
 #include "IO/silo.h"
 #include "common/MPI.h"
 #include "common/Utilities.h"

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef SILO_INTERFACE
 #define SILO_INTERFACE
 

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 <http://www.gnu.org/licenses/>.
-*/
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef SILO_INTERFACE_HPP
 #define SILO_INTERFACE_HPP
 

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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_StackTrace
 #define included_StackTrace
 

analysis/ElectroChemistry.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "analysis/ElectroChemistry.h"
 
 ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr<Domain> dm)
@@ -49,7 +65,7 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr<Domain> dm)
     IonFluxElectrical_y.fill(0);
     IonFluxElectrical_z.resize(Nx, Ny, Nz);
     IonFluxElectrical_z.fill(0);
-    
+
     if (Dm->rank() == 0) {
         bool WriteHeader = false;
         TIMELOG = fopen("electrokinetic.csv", "r");
@@ -75,9 +91,11 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel)
     Nz = Dm->Nz;
     Volume = (Nx - 2) * (Ny - 2) * (Nz - 2) * Dm->nprocx() * Dm->nprocy() *
              Dm->nprocz() * 1.0;
-    
-    if (Dm->rank()==0) printf("Analyze system with sub-domain size = %i x %i x %i \n",Nx,Ny,Nz);
-    
+
+    if (Dm->rank() == 0)
+        printf("Analyze system with sub-domain size = %i x %i x %i \n", Nx, Ny,
+               Nz);
+
     USE_MEMBRANE = IonModel.USE_MEMBRANE;
 
     ChemicalPotential.resize(Nx, Ny, Nz);
@@ -120,11 +138,11 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel)
     IonFluxElectrical_y.fill(0);
     IonFluxElectrical_z.resize(Nx, Ny, Nz);
     IonFluxElectrical_z.fill(0);
-    
-       
+
     if (Dm->rank() == 0) {
-    	printf("Set up analysis routines for %lu ions \n",IonModel.number_ion_species);
-    	
+        printf("Set up analysis routines for %lu ions \n",
+               IonModel.number_ion_species);
+
         bool WriteHeader = false;
         TIMELOG = fopen("electrokinetic.csv", "r");
         if (TIMELOG != NULL)
@@ -138,9 +156,19 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel)
             //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
             fprintf(TIMELOG, "timestep voltage_out voltage_in ");
             fprintf(TIMELOG, "voltage_out_membrane voltage_in_membrane ");
-            for (size_t i=0; i<IonModel.number_ion_species; i++){
-                fprintf(TIMELOG, "rho_%lu_out rho_%lu_in ",i, i);
-                fprintf(TIMELOG, "rho_%lu_out_membrane rho_%lu_in_membrane ", i, i);
+            for (size_t i = 0; i < IonModel.number_ion_species; i++) {
+                fprintf(TIMELOG, "rho_%lu_out rho_%lu_in ", i, i);
+                fprintf(TIMELOG, "rho_%lu_out_membrane rho_%lu_in_membrane ", i,
+                        i);
+                fprintf(TIMELOG, "jx_%lu_out jx_%lu_in ", i, i);
+                fprintf(TIMELOG, "jx_%lu_out_membrane jx_%lu_in_membrane ", i,
+                        i);
+                fprintf(TIMELOG, "jy_%lu_out jy_%lu_in ", i, i);
+                fprintf(TIMELOG, "jy_%lu_out_membrane jy_%lu_in_membrane ", i,
+                        i);
+                fprintf(TIMELOG, "jz_%lu_out jz_%lu_in ", i, i);
+                fprintf(TIMELOG, "jz_%lu_out_membrane jz_%lu_in_membrane ", i,
+                        i);
             }
             fprintf(TIMELOG, "count_out count_in ");
             fprintf(TIMELOG, "count_out_membrane count_in_membrane\n");
@@ -157,18 +185,16 @@ ElectroChemistryAnalyzer::~ElectroChemistryAnalyzer() {
 void ElectroChemistryAnalyzer::SetParams() {}
 
 void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion,
-                                     ScaLBL_Poisson &Poisson,
-                                     int timestep) {
+                                        ScaLBL_Poisson &Poisson, int timestep) {
 
     int i, j, k;
 
     Poisson.getElectricPotential(ElectricalPotential);
 
-    if (Dm->rank() == 0) 
-    	fprintf(TIMELOG, "%i ", timestep);
+    if (Dm->rank() == 0)
+        fprintf(TIMELOG, "%i ", timestep);
 
-    
- /*   int iq, ip, nq, np, nqm, npm;
+    /*   int iq, ip, nq, np, nqm, npm;
     Ion.MembraneDistance(i,j,k); // inside (-) or outside (+) the ion
     for (int link; link<Ion.IonMembrane->membraneLinkCount; link++){
         int iq = Ion.IonMembrane->membraneLinks[2*link];	
@@ -181,60 +207,74 @@ void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion,
 */
     unsigned long int in_local_count, out_local_count;
     unsigned long int in_global_count, out_global_count;
-    
+
     double value_in_local, value_out_local;
     double value_in_global, value_out_global;
 
     double value_membrane_in_local, value_membrane_out_local;
     double value_membrane_in_global, value_membrane_out_global;
-    
+
+    /* flux terms */
+    double jx_membrane_in_local, jx_membrane_out_local;
+    double jy_membrane_in_local, jy_membrane_out_local;
+    double jz_membrane_in_local, jz_membrane_out_local;
+    double jx_in_local, jx_out_local;
+    double jy_in_local, jy_out_local;
+    double jz_in_local, jz_out_local;
+
+    double jx_membrane_in_global, jx_membrane_out_global;
+    double jy_membrane_in_global, jy_membrane_out_global;
+    double jz_membrane_in_global, jz_membrane_out_global;
+    double jx_in_global, jx_out_global;
+    double jy_in_global, jy_out_global;
+    double jz_in_global, jz_out_global;
+
     unsigned long int membrane_in_local_count, membrane_out_local_count;
     unsigned long int membrane_in_global_count, membrane_out_global_count;
-    
-    double memdist,value;
+
+    double memdist, value, jx, jy, jz;
     in_local_count = 0;
     out_local_count = 0;
     membrane_in_local_count = 0;
     membrane_out_local_count = 0;
-    
+
     value_membrane_in_local = 0.0;
     value_membrane_out_local = 0.0;
     value_in_local = 0.0;
     value_out_local = 0.0;
-    for (k = 1; k < Nz; k++) {
-    	for (j = 1; j < Ny; j++) {
-    		for (i = 1; i < Nx; i++) {
-    			/* electric potential  */
-    			memdist = Ion.MembraneDistance(i,j,k);
-    			value = ElectricalPotential(i,j,k);
-    			if (memdist < 0.0){
-    				// inside the membrane
-    				if (fabs(memdist) < 1.0){
-    					value_membrane_in_local += value;
-    					membrane_in_local_count++;
-    				}
-    				value_in_local += value;
-    				in_local_count++;
+    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 = ElectricalPotential(i, j, k);
+                if (memdist < 0.0) {
+                    // inside the membrane
+                    if (fabs(memdist) < 1.0) {
+                        value_membrane_in_local += value;
+                        membrane_in_local_count++;
+                    }
+                    value_in_local += value;
+                    in_local_count++;
 
-    			}
-    			else {
-    				// outside the membrane
-    				if (fabs(memdist) < 1.0){
-    					value_membrane_out_local += value;
-    					membrane_out_local_count++;
-    				}
-    				value_out_local += value;
-    				out_local_count++;
-    			}
-    		}
-    	}
+                } else {
+                    // outside the membrane
+                    if (fabs(memdist) < 1.0) {
+                        value_membrane_out_local += value;
+                        membrane_out_local_count++;
+                    }
+                    value_out_local += value;
+                    out_local_count++;
+                }
+            }
+        }
     }
     /* these only need to be computed the first time through */
     out_global_count = Dm->Comm.sumReduce(out_local_count);
     in_global_count = Dm->Comm.sumReduce(in_local_count);
     membrane_out_global_count = Dm->Comm.sumReduce(membrane_out_local_count);
     membrane_in_global_count = Dm->Comm.sumReduce(membrane_in_local_count);
-    
+
     value_out_global = Dm->Comm.sumReduce(value_out_local);
     value_in_global = Dm->Comm.sumReduce(value_in_local);
     value_membrane_out_global = Dm->Comm.sumReduce(value_membrane_out_local);
@@ -244,12 +284,12 @@ 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;
-    
+
     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 ", value_out_global);
+        fprintf(TIMELOG, "%.8g ", value_in_global);
+        fprintf(TIMELOG, "%.8g ", value_membrane_out_global);
+        fprintf(TIMELOG, "%.8g ", value_membrane_in_global);
     }
 
     value_membrane_in_local = 0.0;
@@ -258,36 +298,74 @@ 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++) {
-        	for (j = 1; j < Ny; j++) {
-        		for (i = 1; i < Nx; i++) {
-        			/* electric potential  */
-        			memdist = Ion.MembraneDistance(i,j,k);
-        			value = Rho(i,j,k);
-        			if (memdist < 0.0){
-        				// inside the membrane
-        				if (fabs(memdist) < 1.0){
-        					value_membrane_in_local += value;
-        				}
-        				value_in_local += value;
-        			}
-        			else {
-        				// outside the membrane
-        				if (fabs(memdist) < 1.0){
-        					value_membrane_out_local += value;
-        				}
-        				value_out_local += value;
-        			}
-        		}
-        	}
+
+        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;
+                    }
+                }
+            }
         }
         value_out_global = Dm->Comm.sumReduce(value_out_local);
         value_in_global = Dm->Comm.sumReduce(value_in_local);
-        value_membrane_out_global = Dm->Comm.sumReduce(value_membrane_out_local);
+        value_membrane_out_global =
+            Dm->Comm.sumReduce(value_membrane_out_local);
         value_membrane_in_global = Dm->Comm.sumReduce(value_membrane_in_local);
 
         value_out_global /= out_global_count;
@@ -295,24 +373,68 @@ void ElectroChemistryAnalyzer::Membrane(ScaLBL_IonModel &Ion,
         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 ", 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);
         }
     }
 
     if (Dm->rank() == 0) {
-    	fprintf(TIMELOG, "%lu ", out_global_count);
-    	fprintf(TIMELOG, "%lu ", in_global_count);
-    	fprintf(TIMELOG, "%lu ", membrane_out_global_count);
-    	fprintf(TIMELOG, "%lu\n", membrane_in_global_count);
+        fprintf(TIMELOG, "%lu ", out_global_count);
+        fprintf(TIMELOG, "%lu ", in_global_count);
+        fprintf(TIMELOG, "%lu ", membrane_out_global_count);
+        fprintf(TIMELOG, "%lu\n", membrane_in_global_count);
         fflush(TIMELOG);
     }
 }
 
-
 void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion,
                                      ScaLBL_Poisson &Poisson,
                                      ScaLBL_StokesModel &Stokes, int timestep) {
@@ -430,15 +552,41 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion,
 
     auto vis_db = input_db->getDatabase("Visualization");
     char VisName[40];
-    auto format = vis_db->getWithDefault<string>( "format", "hdf5" );
+    auto format = vis_db->getWithDefault<string>("format", "hdf5");
+
+    if (Dm->rank() == 0) {
+        printf("ElectroChemistryAnalyzer::WriteVis (format = %s)\n",
+               format.c_str());
+        if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
+            printf("       save electric potential \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_concentration", true)) {
+            printf("       save concentration \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_velocity", false)) {
+            printf("       save velocity \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_diffusive", false)) {
+            printf("       save ion flux (diffusive) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_advective", false)) {
+            printf("       save ion flux (advective) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_electrical", false)) {
+            printf("       save ion flux (electrical) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_electric_field", false)) {
+            printf("       save electric field \n");
+        }
+    }
 
     std::vector<IO::MeshDataStruct> visData;
     fillHalo<double> fillData(Dm->Comm, Dm->rank_info,
                               {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1},
                               0, 1);
 
-    IO::initialize("",format,"false");
-    // Create the MeshDataStruct    
+    IO::initialize("", format, "false");
+    // Create the MeshDataStruct
     visData.resize(1);
 
     visData[0].meshName = "domain";
@@ -835,8 +983,7 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion,
 }
 
 void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion,
-                                     ScaLBL_Poisson &Poisson,
-                                     int timestep) {
+                                     ScaLBL_Poisson &Poisson, int timestep) {
 
     int i, j, k;
     double Vin = 0.0;
@@ -855,10 +1002,10 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion,
     double *rho_mu_fluctuation_global;
     double *rho_psi_avg_global;
     double *rho_psi_fluctuation_global;
-    
+
     /* Get the distance to the membrane */
-    if (Ion.USE_MEMBRANE){
-    	//Ion.MembraneDistance;
+    if (Ion.USE_MEMBRANE) {
+        //Ion.MembraneDistance;
     }
 
     /* local sub-domain averages */
@@ -955,15 +1102,41 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion,
 
     auto vis_db = input_db->getDatabase("Visualization");
     char VisName[40];
-    auto format = vis_db->getWithDefault<string>( "format", "hdf5" );
+    auto format = vis_db->getWithDefault<string>("format", "hdf5");
 
     std::vector<IO::MeshDataStruct> visData;
     fillHalo<double> fillData(Dm->Comm, Dm->rank_info,
                               {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1},
                               0, 1);
 
-    IO::initialize("",format,"false");
-    // Create the MeshDataStruct    
+    if (Dm->rank() == 0) {
+        printf("ElectroChemistryAnalyzer::WriteVis (format = %s)\n",
+               format.c_str());
+        if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
+            printf("       save electric potential \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_concentration", true)) {
+            printf("       save concentration \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_velocity", false)) {
+            printf("       save velocity \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_diffusive", false)) {
+            printf("       save ion flux (diffusive) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_advective", false)) {
+            printf("       save ion flux (advective) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_ion_flux_electrical", false)) {
+            printf("       save ion flux (electrical) \n");
+        }
+        if (vis_db->getWithDefault<bool>("save_electric_field", false)) {
+            printf("       save electric field \n");
+        }
+    }
+
+    IO::initialize("", format, "false");
+    // Create the MeshDataStruct
     visData.resize(1);
 
     visData[0].meshName = "domain";
@@ -1055,7 +1228,6 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion,
         }
     }
 
-
     if (vis_db->getWithDefault<bool>("save_ion_flux_electrical", false)) {
         for (size_t ion = 0; ion < Ion.number_ion_species; ion++) {
             // x-component of electro-migrational flux

analysis/ElectroChemistry.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /*
  *  averaging tools for electrochemistry
  */
@@ -29,7 +45,7 @@ public:
     double nu_n, nu_w;
     double gamma_wn, beta;
     double Fx, Fy, Fz;
-    
+
     bool USE_MEMBRANE;
 
     //...........................................................................
@@ -56,14 +72,16 @@ public:
     DoubleArray IonFluxElectrical_z;
 
     ElectroChemistryAnalyzer(std::shared_ptr<Domain> Dm);
-    ElectroChemistryAnalyzer( ScaLBL_IonModel &IonModel);
+    ElectroChemistryAnalyzer(ScaLBL_IonModel &IonModel);
     ~ElectroChemistryAnalyzer();
 
     void SetParams();
-    void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, int timestep);
+    void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
+               ScaLBL_StokesModel &Stokes, int timestep);
     void Membrane(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, int timestep);
-    void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, 
-    		ScaLBL_StokesModel &Stokes,std::shared_ptr<Database> input_db, int timestep);
+    void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
+                  ScaLBL_StokesModel &Stokes,
+                  std::shared_ptr<Database> input_db, int timestep);
     void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, int timestep);
     void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson,
                   std::shared_ptr<Database> input_db, int timestep);

analysis/FlowAdaptor.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* Flow adaptor class for multiphase flow methods */
 
 #include "analysis/FlowAdaptor.h"
@@ -17,9 +33,7 @@ FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M) {
     phi_t.fill(0); // time derivative for the phase indicator field
 }
 
-FlowAdaptor::~FlowAdaptor() {
-	
-}
+FlowAdaptor::~FlowAdaptor() {}
 
 double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename) {
     int rank = M.rank;
@@ -237,12 +251,12 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M) {
     //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
     ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7 * Np * sizeof(double));
     ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7 * Np * sizeof(double));
-    
-    delete Aq_tmp; 
+
+    delete Aq_tmp;
     delete Bq_tmp;
-    delete Vel_x; 
-    delete Vel_y; 
-    delete Vel_z; 
+    delete Vel_x;
+    delete Vel_y;
+    delete Vel_z;
     delete Phase;
 
     return (TOTAL_MASS_CHANGE);
@@ -594,8 +608,8 @@ double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M,
     //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
     ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7 * Np * sizeof(double));
     ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7 * Np * sizeof(double));
-    
-    delete Aq_tmp; 
+
+    delete Aq_tmp;
     delete Bq_tmp;
     return (mass_loss);
 }

analysis/FlowAdaptor.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* Flow adaptor class for multiphase flow methods */
 
 #ifndef ScaLBL_FlowAdaptor_INC

analysis/FreeEnergy.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "analysis/FreeEnergy.h"
 
 FreeEnergyAnalyzer::FreeEnergyAnalyzer(std::shared_ptr<Domain> dm) : Dm(dm) {

analysis/FreeEnergy.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /*
  *  averaging tools for electrochemistry
  */

analysis/GreyPhase.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "analysis/GreyPhase.h"
 
 // Constructor

analysis/GreyPhase.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /*
  * Sub-phase averaging tools
  */

analysis/Minkowski.cpp

@@ -161,7 +161,7 @@ void Minkowski::MeasureObject() {
 	 *    1 - labels the rest of the 
 	 */
     //DoubleArray smooth_distance(Nx,Ny,Nz);
-	
+
     for (int k = 0; k < Nz; k++) {
         for (int j = 0; j < Ny; j++) {
             for (int i = 0; i < Nx; i++) {

analysis/SubPhase.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "analysis/SubPhase.h"
 
 // Constructor
@@ -411,7 +427,7 @@ void SubPhase::Basic() {
             dir_z = 1.0;
             force_mag = 1.0;
         }
-        double Porosity = (gwb.V + gnb.V)/Dm->Volume;
+        double Porosity = (gwb.V + gnb.V) / Dm->Volume;
         double saturation = gwb.V / (gwb.V + gnb.V);
         double water_flow_rate =
             gwb.V * (gwb.Px * dir_x + gwb.Py * dir_y + gwb.Pz * dir_z) / gwb.M /
@@ -430,11 +446,13 @@ 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,

analysis/SubPhase.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /*
  * Sub-phase averaging tools
  */

analysis/TwoPhase.cpp

@@ -27,7 +27,6 @@
 #include "IO/Writer.h"
 #include "analysis/filters.h"
 
-
 #include <memory>
 
 #define BLOB_AVG_COUNT 35
@@ -417,7 +416,8 @@ void TwoPhase::UpdateSolid() {
 
 void TwoPhase::UpdateMeshValues() {
     int i, j, k, n;
-    fillHalo<double> fillData(Dm->Comm, Dm->rank_info, {Nx-2,Ny-2,Nz-2}, {1, 1, 1}, 0, 1);
+    fillHalo<double> fillData(Dm->Comm, Dm->rank_info, {Nx - 2, Ny - 2, Nz - 2},
+                              {1, 1, 1}, 0, 1);
 
     //...........................................................................
     //Dm->CommunicateMeshHalo(SDn);
@@ -576,7 +576,7 @@ void TwoPhase::ComputeLocal() {
                     Kwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
-                    
+
                     Jwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
@@ -607,7 +607,7 @@ void TwoPhase::ComputeLocal() {
                     efawns += pmmc_CubeContactAngle(
                         CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
                         SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
-                    
+
                     wwnsdnwn += pmmc_CommonCurveSpeed(
                         CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z, SDs_x,
                         SDs_y, SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
@@ -721,18 +721,19 @@ void TwoPhase::ComputeStatic() {
     kmin = 1;
     kmax = Nz - 1;
     imin = jmin = 1;
-    
+
     /* set fluid isovalue to "grow" NWP for contact angle measurement */
     fluid_isovalue = -1.0;
-    
+
     string FILENAME = "ContactAngle";
-    
+
     char LocalRankString[8];
     char LocalRankFilename[40];
     sprintf(LocalRankString, "%05d", Dm->rank());
-    sprintf(LocalRankFilename, "%s%s%s", "ContactAngle.", LocalRankString,".csv");
+    sprintf(LocalRankFilename, "%s%s%s", "ContactAngle.", LocalRankString,
+            ".csv");
     FILE *ANGLES = fopen(LocalRankFilename, "a+");
-    fprintf(ANGLES,"x y z angle\n");
+    fprintf(ANGLES, "x y z angle\n");
 
     for (k = kmin; k < kmax; k++) {
         for (j = jmin; j < Ny - 1; j++) {
@@ -777,13 +778,13 @@ void TwoPhase::ComputeStatic() {
                     Kwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
-                    
+
                     Jwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
-                    
-                    Xwn += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
-                                                         n_nw_tris, i, j, k);
+
+                    Xwn += geomavg_EulerCharacteristic(
+                        nw_pts, nw_tris, n_nw_pts, n_nw_tris, i, j, k);
 
                     // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
                     pmmc_CubeTrimSurfaceInterpValues(
@@ -801,12 +802,13 @@ void TwoPhase::ComputeStatic() {
                     efawns += pmmc_CubeContactAngle(
                         CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
                         SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
-                    
+
                     for (int p = 0; p < n_local_nws_pts; p++) {
                         // Extract the line segment
                         Point A = local_nws_pts(p);
                         double value = Values(p);
-                        fprintf(ANGLES, "%.8g %.8g %.8g %.8g\n", A.x, A.y, A.z, value);
+                        fprintf(ANGLES, "%.8g %.8g %.8g %.8g\n", A.x, A.y, A.z,
+                                value);
                     }
 
                     pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x,
@@ -816,14 +818,14 @@ void TwoPhase::ComputeStatic() {
 
                     lwns +=
                         pmmc_CubeCurveLength(local_nws_pts, n_local_nws_pts);
-                    
+
                     /* half contribution for vertices / edges at the common line 
                      *   each cube with contact line has a net of undercounting vertices 
                      *   each cube is undercounting edges due to internal counts
                     */
-                    Xwn += 0.25*n_local_nws_pts - 0.5;
-                    Xws += 0.25*n_local_nws_pts - 0.5;
-                    Xns += 0.25*n_local_nws_pts - 0.5;
+                    Xwn += 0.25 * n_local_nws_pts - 0.5;
+                    Xws += 0.25 * n_local_nws_pts - 0.5;
+                    Xns += 0.25 * n_local_nws_pts - 0.5;
                 }
 
                 // Solid interface averagees
@@ -836,12 +838,12 @@ void TwoPhase::ComputeStatic() {
                                                        n_ns_tris);
                     aws += pmmc_CubeSurfaceOrientation(Gws, ws_pts, ws_tris,
                                                        n_ws_tris);
-                    
-                    Xws += geomavg_EulerCharacteristic(ws_pts, ws_tris, n_ws_pts,
-                                                         n_ws_tris, i, j, k);
-                    
-                    Xns += geomavg_EulerCharacteristic(ns_pts, ns_tris, n_ns_pts,
-                                                         n_ns_tris, i, j, k);
+
+                    Xws += geomavg_EulerCharacteristic(
+                        ws_pts, ws_tris, n_ws_pts, n_ws_tris, i, j, k);
+
+                    Xns += geomavg_EulerCharacteristic(
+                        ns_pts, ns_tris, n_ns_pts, n_ns_tris, i, j, k);
                 }
                 //...........................................................................
                 // Compute the integral curvature of the non-wetting phase
@@ -866,11 +868,9 @@ void TwoPhase::ComputeStatic() {
                 Kn += pmmc_CubeSurfaceInterpValue(CubeValues, GaussCurvature,
                                                   nw_pts, nw_tris, Values, i, j,
                                                   k, n_nw_pts, n_nw_tris);
-                
 
                 euler += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
                                                      n_nw_tris, i, j, k);
-
             }
         }
     }
@@ -1538,7 +1538,6 @@ void TwoPhase::Reduce() {
     dEs = dEs * iVol_global;
     lwns_global = lwns_global * iVol_global;
     */
-
 }
 
 void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
@@ -1552,27 +1551,28 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
 
 void TwoPhase::PrintStatic() {
     if (Dm->rank() == 0) {
-    	FILE *STATIC;
+        FILE *STATIC;
         STATIC = fopen("geometry.csv", "a+");
         if (fseek(STATIC, 0, SEEK_SET) == fseek(STATIC, 0, SEEK_CUR)) {
             // If timelog is empty, write a short header to list the averages
             fprintf(STATIC, "sw awn ans aws Jwn Kwn lwns cwns KGws "
-                             "KGwn Xwn Xws Xns "); // Scalar averages
-            fprintf(STATIC,
+                            "KGwn Xwn Xws Xns "); // Scalar averages
+            fprintf(
+                STATIC,
                 "Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
             fprintf(STATIC, "Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
             fprintf(STATIC, "Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
             fprintf(STATIC, "trawn trJwn trRwn "); //trimmed curvature,
-            fprintf(STATIC, "Vw Aw Jw Xw ");      //miknowski measures,
-            fprintf(STATIC, "Vn An Jn Xn\n");     //miknowski measures,
+            fprintf(STATIC, "Vw Aw Jw Xw ");       //miknowski measures,
+            fprintf(STATIC, "Vn An Jn Xn\n");      //miknowski measures,
             //fprintf(STATIC,"Euler Kn2 Jn2 An2\n"); 			//miknowski measures,
         }
 
-        fprintf(STATIC, "%.5g ", sat_w); // saturation 
+        fprintf(STATIC, "%.5g ", sat_w); // saturation
         fprintf(STATIC, "%.5g %.5g %.5g ", awn_global, ans_global,
                 aws_global); // interfacial areas
         fprintf(STATIC, "%.5g %.5g ", Jwn_global,
-                Kwn_global);                      // curvature of wn interface
+                Kwn_global);                     // curvature of wn interface
         fprintf(STATIC, "%.5g ", lwns_global);   // common curve length
         fprintf(STATIC, "%.5g ", efawns_global); // average contact angle
         fprintf(STATIC, "%.5g %.5g ", KNwns_global,
@@ -1592,7 +1592,7 @@ void TwoPhase::PrintStatic() {
                 trRwn_global); // Trimmed curvature
         fprintf(STATIC, "%.5g %.5g %.5g %.5g ", wet_morph->V(), wet_morph->A(),
                 wet_morph->H(), wet_morph->X());
-       fprintf(STATIC, "%.5g %.5g %.5g %.5g\n", nonwet_morph->V(),
+        fprintf(STATIC, "%.5g %.5g %.5g %.5g\n", nonwet_morph->V(),
                 nonwet_morph->A(), nonwet_morph->H(), nonwet_morph->X());
         //fprintf(STATIC,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global);			// minkowski measures
         fclose(STATIC);

analysis/analysis.cpp

@@ -14,7 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-
 #include "analysis/analysis.h"
 #include "ProfilerApp.h"
 

analysis/dcel.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "analysis/dcel.h"
 
 DCEL::DCEL() {}

analysis/dcel.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #ifndef DCEL_INC
 #define DCEL_INC
 

analysis/imfilter.hpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #include "analysis/imfilter.h"
 #include "ProfilerApp.h"
 #include <math.h>

analysis/morphology.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <analysis/morphology.h>
 // Implementation of morphological opening routine
 
@@ -137,7 +153,8 @@ void Morphology::Initialize(std::shared_ptr<Domain> Dm, DoubleArray &Distance) {
     morphRadius.resize(recvLoc);
     //..............................
     /* send the morphological radius */
-    Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 0);
+    Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(),
+                   recvtag + 0);
     Dm->Comm.send(&tmpDistance[0], sendCount, Dm->rank_X(), sendtag + 0);
     /* send the shift values */
     Dm->Comm.Irecv(&xShift[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 1);
@@ -501,7 +518,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
     if (rank == 0)
         printf("Maximum pore size: %f \n", maxdistGlobal);
     final_void_fraction = volume_fraction; //initialize
-    
+
     int ii, jj, kk;
     int imin, jmin, kmin, imax, jmax, kmax;
     int Nx = nx;
@@ -523,28 +540,30 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
 
     int numTry = 0;
     int maxTry = 100;
-    while ( !(void_fraction_new < VoidFraction) && numTry < maxTry) {
+    while (!(void_fraction_new < VoidFraction) && numTry < maxTry) {
         numTry++;
         void_fraction_diff_old = void_fraction_diff_new;
         void_fraction_old = void_fraction_new;
         Rcrit_old = Rcrit_new;
         Rcrit_new -= deltaR * Rcrit_old;
-        if (rank==0) printf("Try %i with radius %f \n", numTry, Rcrit_new);
+        if (rank == 0)
+            printf("Try %i with radius %f \n", numTry, Rcrit_new);
         if (Rcrit_new < 0.5) {
             numTry = maxTry;
         }
         int Window = round(Rcrit_new);
         if (Window == 0)
-            Window =  1; // If Window = 0 at the begining, after the following process will have sw=1.0
+            Window =
+                1; // If Window = 0 at the begining, after the following process will have sw=1.0
         // and sw<Sw will be immediately broken
         double LocalNumber = 0.f;
-        for (int k = 1; k < Nz-1; k++) {
-            for (int j = 1; j < Ny-1; j++) {
-                for (int i = 1; i < Nx-1; i++) {
+        for (int k = 1; k < Nz - 1; k++) {
+            for (int j = 1; j < Ny - 1; j++) {
+                for (int i = 1; i < Nx - 1; i++) {
                     n = k * nx * ny + j * nx + i;
                     if (SignDist(i, j, k) > Rcrit_new) {
                         // loop over the window and update
-                    	//printf("Distance(%i %i %i) = %f \n",i,j,k, SignDist(i,j,k));
+                        //printf("Distance(%i %i %i) = %f \n",i,j,k, SignDist(i,j,k));
                         imin = max(1, i - Window);
                         jmin = max(1, j - Window);
                         kmin = max(1, k - Window);
@@ -611,7 +630,8 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
 
 //***************************************************************************************
 double MorphDrain(DoubleArray &SignDist, signed char *id,
-                  std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius) {
+                  std::shared_ptr<Domain> Dm, double VoidFraction,
+                  double InitialRadius) {
     // SignDist is the distance to the object that you want to constaing the morphological opening
     // VoidFraction is the the empty space where the object inst
     // id is a labeled map
@@ -688,10 +708,10 @@ double MorphDrain(DoubleArray &SignDist, signed char *id,
     double deltaR = 0.05; // amount to change the radius in voxel units
     double Rcrit_old = maxdistGlobal;
     double Rcrit_new = maxdistGlobal;
-    
-    if (InitialRadius < maxdistGlobal){
-    	Rcrit_old = InitialRadius;
-    	Rcrit_new = InitialRadius;
+
+    if (InitialRadius < maxdistGlobal) {
+        Rcrit_old = InitialRadius;
+        Rcrit_new = InitialRadius;
     }
     //if (argc>2){
     //	Rcrit_new = strtod(argv[2],NULL);

analysis/morphology.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // Morphological opening routine
 #include "common/Array.h"
 #include "common/Domain.h"
@@ -7,7 +23,8 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
                  std::shared_ptr<Domain> Dm, double VoidFraction,
                  signed char ErodeLabel, signed char ReplaceLabel);
 double MorphDrain(DoubleArray &SignDist, signed char *id,
-                  std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius);
+                  std::shared_ptr<Domain> Dm, double VoidFraction,
+                  double InitialRadius);
 double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
                  std::shared_ptr<Domain> Dm, double TargetVol,
                  double WallFactor);

analysis/pmmc.h

@@ -4057,7 +4057,7 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues,
                       (A.z - B.z) * (A.z - B.z));
         integral += 0.5 * length * (vA + vB);
     }
-    
+
     return integral;
 }
 //--------------------------------------------------------------------------------------------------------
@@ -4438,12 +4438,12 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
                 fx = f_x(i, j, k);
                 fy = f_y(i, j, k);
                 fz = f_z(i, j, k);
-                
+
                 // Normal to fluid surface
                 Nx.Corners(i - ic, j - jc, k - kc) = fx;
                 Ny.Corners(i - ic, j - jc, k - kc) = fy;
                 Nz.Corners(i - ic, j - jc, k - kc) = fz;
-                
+
                 // Normal to solid surface
                 Sx.Corners(i - ic, j - jc, k - kc) = sx;
                 Sy.Corners(i - ic, j - jc, k - kc) = sy;
@@ -4550,7 +4550,7 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
         nsx /= norm;
         nsy /= norm;
         nsz /= norm;
-        
+
         // Normal vector to the fluid surface
         nwx = Nx.eval(P);
         nwy = Ny.eval(P);
@@ -4578,7 +4578,7 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
             nwsy = -nwsy;
             nwsz = -nwsz;
         }
-        
+
         // common curve normal in the fluid surface tangent plane (rel. geodesic curvature)
         nwnx = twnsy * nwz - twnsz * nwy;
         nwny = twnsz * nwx - twnsx * nwz;
@@ -4596,7 +4596,6 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
             nwnz = -nwnz;
         }
 
-
         if (length > 0.0) {
             // normal curvature component in the direction of the solid surface
             //KNavg += K * (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz) * length;

cmake/macros.cmake

@@ -196,9 +196,6 @@ MACRO( FIND_FILES )
     # Find the HIP sources
     SET( T_HIPSOURCES "" )
     FILE( GLOB T_HIPSOURCES "*.hip" )
-    # Find the SYCL sources
-    SET( T_SYCLSOURCES "" )
-    FILE( GLOB T_SYCLSOURCES "*.dp.cpp" )
     # Find the C sources
     SET( T_CSOURCES "" )
     FILE( GLOB T_CSOURCES "*.c" )
@@ -219,11 +216,10 @@ MACRO( FIND_FILES )
     SET( CXXSOURCES ${CXXSOURCES} ${T_CXXSOURCES} )
     SET( CUDASOURCES ${CUDASOURCES} ${T_CUDASOURCES} )
     SET( HIPSOURCES ${HIPSOURCES} ${T_HIPSOURCES} )
-    SET( SYCLSOURCES ${SYCLSOURCES} ${T_SYCLSOURCES} )
     SET( CSOURCES ${CSOURCES} ${T_CSOURCES} )
     SET( FSOURCES ${FSOURCES} ${T_FSOURCES} )
     SET( M4FSOURCES ${M4FSOURCES} ${T_M4FSOURCES} )
-    SET( SOURCES ${SOURCES} ${T_CXXSOURCES} ${T_CSOURCES} ${T_FSOURCES} ${T_M4FSOURCES} ${CUDASOURCES} ${HIPSOURCES} ${SYCLSOURCES})
+    SET( SOURCES ${SOURCES} ${T_CXXSOURCES} ${T_CSOURCES} ${T_FSOURCES} ${T_M4FSOURCES} ${CUDASOURCES} ${HIPSOURCES} )
 ENDMACRO()
 
 
@@ -238,9 +234,6 @@ MACRO( FIND_FILES_PATH IN_PATH )
     # Find the HIP sources
     SET( T_HIPSOURCES "" )
     FILE( GLOB T_HIPSOURCES "${IN_PATH}/*.hip" )
-    # Find the SYCL sources
-    SET( T_SYCLSOURCES "" )
-    FILE( GLOB T_SYCLSOURCES "${IN_PATH}/*.sycl" )
     # Find the C sources
     SET( T_CSOURCES "" )
     FILE( GLOB T_CSOURCES "${IN_PATH}/*.c" )
@@ -261,10 +254,9 @@ MACRO( FIND_FILES_PATH IN_PATH )
     SET( CXXSOURCES ${CXXSOURCES} ${T_CXXSOURCES} )
     SET( CUDASOURCES ${CUDASOURCES} ${T_CUDASOURCES} )
     SET( HIPSOURCES ${HIPSOURCES} ${T_HIPSOURCES} )
-    SET( SYCLSOURCES ${SYCLSOURCES} ${T_SYCLSOURCES} )
     SET( CSOURCES ${CSOURCES} ${T_CSOURCES} )
     SET( FSOURCES ${FSOURCES} ${T_FSOURCES} )
-    SET( SOURCES ${SOURCES} ${T_CXXSOURCES} ${T_CSOURCES} ${T_FSOURCES} ${CUDASOURCES} ${HIPSOURCES} ${SYCLSOURCES} )
+    SET( SOURCES ${SOURCES} ${T_CXXSOURCES} ${T_CSOURCES} ${T_FSOURCES} ${CUDASOURCES} ${HIPSOURCES} )
 ENDMACRO()
 
 

common/Array.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // clang-format off
 #include "common/Array.h"
 #include "common/Array.hpp"

common/Array.hpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_ArrayClass_hpp
 #define included_ArrayClass_hpp
 

common/ArraySize.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #ifndef included_ArraySizeClass
 #define included_ArraySizeClass
 

common/Communication.h

@@ -208,68 +208,72 @@ inline void CommunicateSendRecvCounts(
 }
 
 //***************************************************************************************
-inline void CommunicateRecvLists( const Utilities::MPI& comm, int sendtag, int recvtag, 
-		int *sendList_x, int *sendList_y, int *sendList_z, int *sendList_X, int *sendList_Y, int *sendList_Z,
-		int *sendList_xy, int *sendList_XY, int *sendList_xY, int *sendList_Xy,
-		int *sendList_xz, int *sendList_XZ, int *sendList_xZ, int *sendList_Xz,
-		int *sendList_yz, int *sendList_YZ, int *sendList_yZ, int *sendList_Yz,
-		int sendCount_x, int sendCount_y, int sendCount_z, int sendCount_X, int sendCount_Y, int sendCount_Z,
-		int sendCount_xy, int sendCount_XY, int sendCount_xY, int sendCount_Xy,
-		int sendCount_xz, int sendCount_XZ, int sendCount_xZ, int sendCount_Xz,
-		int sendCount_yz, int sendCount_YZ, int sendCount_yZ, int sendCount_Yz,
-		int *recvList_x, int *recvList_y, int *recvList_z, int *recvList_X, int *recvList_Y, int *recvList_Z,
-		int *recvList_xy, int *recvList_XY, int *recvList_xY, int *recvList_Xy,
-		int *recvList_xz, int *recvList_XZ, int *recvList_xZ, int *recvList_Xz,
-		int *recvList_yz, int *recvList_YZ, int *recvList_yZ, int *recvList_Yz,
-		int recvCount_x, int recvCount_y, int recvCount_z, int recvCount_X, int recvCount_Y, int recvCount_Z,
-		int recvCount_xy, int recvCount_XY, int recvCount_xY, int recvCount_Xy,
-		int recvCount_xz, int recvCount_XZ, int recvCount_xZ, int recvCount_Xz,
-		int recvCount_yz, int recvCount_YZ, int recvCount_yZ, int recvCount_Yz,
-		int rank_x, int rank_y, int rank_z, int rank_X, int rank_Y, int rank_Z, int rank_xy, int rank_XY, int rank_xY,
-		int rank_Xy, int rank_xz, int rank_XZ, int rank_xZ, int rank_Xz, int rank_yz, int rank_YZ, int rank_yZ, int rank_Yz)
-{
-	MPI_Request req1[18], req2[18];
-	req1[0] = comm.Isend(sendList_x,sendCount_x,rank_x,sendtag+0);
-	req2[0] = comm.Irecv(recvList_X,recvCount_X,rank_X,recvtag+0);
-	req1[1] = comm.Isend(sendList_X,sendCount_X,rank_X,sendtag+1);
-	req2[1] = comm.Irecv(recvList_x,recvCount_x,rank_x,recvtag+1);
-	req1[2] = comm.Isend(sendList_y,sendCount_y,rank_y,sendtag+2);
-	req2[2] = comm.Irecv(recvList_Y,recvCount_Y,rank_Y,recvtag+2);
-	req1[3] = comm.Isend(sendList_Y,sendCount_Y,rank_Y,sendtag+3);
-	req2[3] = comm.Irecv(recvList_y,recvCount_y,rank_y,recvtag+3);
-	req1[4] = comm.Isend(sendList_z,sendCount_z,rank_z,sendtag+4);
-	req2[4] = comm.Irecv(recvList_Z,recvCount_Z,rank_Z,recvtag+4);
-	req1[5] = comm.Isend(sendList_Z,sendCount_Z,rank_Z,sendtag+5);
-	req2[5] = comm.Irecv(recvList_z,recvCount_z,rank_z,recvtag+5);
+inline void CommunicateRecvLists(
+    const Utilities::MPI &comm, int sendtag, int recvtag, int *sendList_x,
+    int *sendList_y, int *sendList_z, int *sendList_X, int *sendList_Y,
+    int *sendList_Z, int *sendList_xy, int *sendList_XY, int *sendList_xY,
+    int *sendList_Xy, int *sendList_xz, int *sendList_XZ, int *sendList_xZ,
+    int *sendList_Xz, int *sendList_yz, int *sendList_YZ, int *sendList_yZ,
+    int *sendList_Yz, int sendCount_x, int sendCount_y, int sendCount_z,
+    int sendCount_X, int sendCount_Y, int sendCount_Z, int sendCount_xy,
+    int sendCount_XY, int sendCount_xY, int sendCount_Xy, int sendCount_xz,
+    int sendCount_XZ, int sendCount_xZ, int sendCount_Xz, int sendCount_yz,
+    int sendCount_YZ, int sendCount_yZ, int sendCount_Yz, int *recvList_x,
+    int *recvList_y, int *recvList_z, int *recvList_X, int *recvList_Y,
+    int *recvList_Z, int *recvList_xy, int *recvList_XY, int *recvList_xY,
+    int *recvList_Xy, int *recvList_xz, int *recvList_XZ, int *recvList_xZ,
+    int *recvList_Xz, int *recvList_yz, int *recvList_YZ, int *recvList_yZ,
+    int *recvList_Yz, int recvCount_x, int recvCount_y, int recvCount_z,
+    int recvCount_X, int recvCount_Y, int recvCount_Z, int recvCount_xy,
+    int recvCount_XY, int recvCount_xY, int recvCount_Xy, int recvCount_xz,
+    int recvCount_XZ, int recvCount_xZ, int recvCount_Xz, int recvCount_yz,
+    int recvCount_YZ, int recvCount_yZ, int recvCount_Yz, int rank_x,
+    int rank_y, int rank_z, int rank_X, int rank_Y, int rank_Z, int rank_xy,
+    int rank_XY, int rank_xY, int rank_Xy, int rank_xz, int rank_XZ,
+    int rank_xZ, int rank_Xz, int rank_yz, int rank_YZ, int rank_yZ,
+    int rank_Yz) {
+    MPI_Request req1[18], req2[18];
+    req1[0] = comm.Isend(sendList_x, sendCount_x, rank_x, sendtag + 0);
+    req2[0] = comm.Irecv(recvList_X, recvCount_X, rank_X, recvtag + 0);
+    req1[1] = comm.Isend(sendList_X, sendCount_X, rank_X, sendtag + 1);
+    req2[1] = comm.Irecv(recvList_x, recvCount_x, rank_x, recvtag + 1);
+    req1[2] = comm.Isend(sendList_y, sendCount_y, rank_y, sendtag + 2);
+    req2[2] = comm.Irecv(recvList_Y, recvCount_Y, rank_Y, recvtag + 2);
+    req1[3] = comm.Isend(sendList_Y, sendCount_Y, rank_Y, sendtag + 3);
+    req2[3] = comm.Irecv(recvList_y, recvCount_y, rank_y, recvtag + 3);
+    req1[4] = comm.Isend(sendList_z, sendCount_z, rank_z, sendtag + 4);
+    req2[4] = comm.Irecv(recvList_Z, recvCount_Z, rank_Z, recvtag + 4);
+    req1[5] = comm.Isend(sendList_Z, sendCount_Z, rank_Z, sendtag + 5);
+    req2[5] = comm.Irecv(recvList_z, recvCount_z, rank_z, recvtag + 5);
 
-	req1[6] = comm.Isend(sendList_xy,sendCount_xy,rank_xy,sendtag+6);
-	req2[6] = comm.Irecv(recvList_XY,recvCount_XY,rank_XY,recvtag+6);
-	req1[7] = comm.Isend(sendList_XY,sendCount_XY,rank_XY,sendtag+7);
-	req2[7] = comm.Irecv(recvList_xy,recvCount_xy,rank_xy,recvtag+7);
-	req1[8] = comm.Isend(sendList_Xy,sendCount_Xy,rank_Xy,sendtag+8);
-	req2[8] = comm.Irecv(recvList_xY,recvCount_xY,rank_xY,recvtag+8);
-	req1[9] = comm.Isend(sendList_xY,sendCount_xY,rank_xY,sendtag+9);
-	req2[9] = comm.Irecv(recvList_Xy,recvCount_Xy,rank_Xy,recvtag+9);
+    req1[6] = comm.Isend(sendList_xy, sendCount_xy, rank_xy, sendtag + 6);
+    req2[6] = comm.Irecv(recvList_XY, recvCount_XY, rank_XY, recvtag + 6);
+    req1[7] = comm.Isend(sendList_XY, sendCount_XY, rank_XY, sendtag + 7);
+    req2[7] = comm.Irecv(recvList_xy, recvCount_xy, rank_xy, recvtag + 7);
+    req1[8] = comm.Isend(sendList_Xy, sendCount_Xy, rank_Xy, sendtag + 8);
+    req2[8] = comm.Irecv(recvList_xY, recvCount_xY, rank_xY, recvtag + 8);
+    req1[9] = comm.Isend(sendList_xY, sendCount_xY, rank_xY, sendtag + 9);
+    req2[9] = comm.Irecv(recvList_Xy, recvCount_Xy, rank_Xy, recvtag + 9);
 
-	req1[10] = comm.Isend(sendList_xz,sendCount_xz,rank_xz,sendtag+10);
-	req2[10] = comm.Irecv(recvList_XZ,recvCount_XZ,rank_XZ,recvtag+10);
-	req1[11] = comm.Isend(sendList_XZ,sendCount_XZ,rank_XZ,sendtag+11);
-	req2[11] = comm.Irecv(recvList_xz,recvCount_xz,rank_xz,recvtag+11);
-	req1[12] = comm.Isend(sendList_Xz,sendCount_Xz,rank_Xz,sendtag+12);
-	req2[12] = comm.Irecv(recvList_xZ,recvCount_xZ,rank_xZ,recvtag+12);
-	req1[13] = comm.Isend(sendList_xZ,sendCount_xZ,rank_xZ,sendtag+13);
-	req2[13] = comm.Irecv(recvList_Xz,recvCount_Xz,rank_Xz,recvtag+13);
+    req1[10] = comm.Isend(sendList_xz, sendCount_xz, rank_xz, sendtag + 10);
+    req2[10] = comm.Irecv(recvList_XZ, recvCount_XZ, rank_XZ, recvtag + 10);
+    req1[11] = comm.Isend(sendList_XZ, sendCount_XZ, rank_XZ, sendtag + 11);
+    req2[11] = comm.Irecv(recvList_xz, recvCount_xz, rank_xz, recvtag + 11);
+    req1[12] = comm.Isend(sendList_Xz, sendCount_Xz, rank_Xz, sendtag + 12);
+    req2[12] = comm.Irecv(recvList_xZ, recvCount_xZ, rank_xZ, recvtag + 12);
+    req1[13] = comm.Isend(sendList_xZ, sendCount_xZ, rank_xZ, sendtag + 13);
+    req2[13] = comm.Irecv(recvList_Xz, recvCount_Xz, rank_Xz, recvtag + 13);
 
-	req1[14] = comm.Isend(sendList_yz,sendCount_yz,rank_yz,sendtag+14);
-	req2[14] = comm.Irecv(recvList_YZ,recvCount_YZ,rank_YZ,recvtag+14);
-	req1[15] = comm.Isend(sendList_YZ,sendCount_YZ,rank_YZ,sendtag+15);
-	req2[15] = comm.Irecv(recvList_yz,recvCount_yz,rank_yz,recvtag+15);
-	req1[16] = comm.Isend(sendList_Yz,sendCount_Yz,rank_Yz,sendtag+16);
-	req2[16] = comm.Irecv(recvList_yZ,recvCount_yZ,rank_yZ,recvtag+16);
-	req1[17] = comm.Isend(sendList_yZ,sendCount_yZ,rank_yZ,sendtag+17);
-	req2[17] = comm.Irecv(recvList_Yz,recvCount_Yz,rank_Yz,recvtag+17);
-    comm.waitAll( 18, req1 );
-    comm.waitAll( 18, req2 );
+    req1[14] = comm.Isend(sendList_yz, sendCount_yz, rank_yz, sendtag + 14);
+    req2[14] = comm.Irecv(recvList_YZ, recvCount_YZ, rank_YZ, recvtag + 14);
+    req1[15] = comm.Isend(sendList_YZ, sendCount_YZ, rank_YZ, sendtag + 15);
+    req2[15] = comm.Irecv(recvList_yz, recvCount_yz, rank_yz, recvtag + 15);
+    req1[16] = comm.Isend(sendList_Yz, sendCount_Yz, rank_Yz, sendtag + 16);
+    req2[16] = comm.Irecv(recvList_yZ, recvCount_yZ, rank_yZ, recvtag + 16);
+    req1[17] = comm.Isend(sendList_yZ, sendCount_yZ, rank_yZ, sendtag + 17);
+    req2[17] = comm.Irecv(recvList_Yz, recvCount_Yz, rank_Yz, recvtag + 17);
+    comm.waitAll(18, req1);
+    comm.waitAll(18, req2);
 }
 
 //***************************************************************************************

common/Communication.hpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef COMMUNICATION_HPP_INC
 #define COMMUNICATION_HPP_INC
 

common/Database.hpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_Database_hpp
 #define included_Database_hpp
 

common/Domain.h

@@ -203,7 +203,7 @@ public: // Public variables (need to create accessors instead)
      * \brief Read domain IDs from file
     */
     void ReadIDs();
-    
+
     /** 
      * \brief Read domain IDs from SWC file
     */

common/FunctionTable.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "FunctionTable.hpp"
 
 /********************************************************
@@ -92,12 +108,11 @@ template<> long double genRand<long double>()
 /********************************************************
  *  axpy                                                 *
  ********************************************************/
-template <>
-void call_axpy<float>(size_t, const float, const float*, float*) {
+template <> void call_axpy<float>(size_t, const float, const float *, float *) {
     ERROR("Not finished");
 }
 template <>
-void call_axpy<double>(size_t, const double, const double*, double*) {
+void call_axpy<double>(size_t, const double, const double *, double *) {
     ERROR("Not finished");
 }
 
@@ -105,22 +120,22 @@ void call_axpy<double>(size_t, const double, const double*, double*) {
  *  Multiply two arrays                                  *
  ********************************************************/
 template <>
-void call_gemv<double>(size_t, size_t, double, double,
-                       const double*, const double*, double*) {
+void call_gemv<double>(size_t, size_t, double, double, const double *,
+                       const double *, double *) {
     ERROR("Not finished");
 }
 template <>
-void call_gemv<float>(size_t, size_t, float, float,
-                      const float*, const float*, float*) {
+void call_gemv<float>(size_t, size_t, float, float, const float *,
+                      const float *, float *) {
     ERROR("Not finished");
 }
 template <>
-void call_gemm<double>(size_t, size_t, size_t, double, double,
-                       const double*, const double*, double*) {
+void call_gemm<double>(size_t, size_t, size_t, double, double, const double *,
+                       const double *, double *) {
     ERROR("Not finished");
 }
 template <>
-void call_gemm<float>(size_t, size_t, size_t, float, float,
-                      const float*, const float*, float*) {
+void call_gemm<float>(size_t, size_t, size_t, float, float, const float *,
+                      const float *, float *) {
     ERROR("Not finished");
 }

common/FunctionTable.hpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #ifndef included_FunctionTable_hpp
 #define included_FunctionTable_hpp
 
@@ -297,10 +281,9 @@ TYPE FunctionTable::sum(const Array<TYPE, FUN, ALLOC> &A) {
 }
 
 template <class TYPE>
-inline void FunctionTable::gemmWrapper(char, char, int, int,
-                                       int, TYPE, const TYPE*,
-                                       int, const TYPE*, int,
-                                       TYPE, TYPE*, int) {
+inline void FunctionTable::gemmWrapper(char, char, int, int, int, TYPE,
+                                       const TYPE *, int, const TYPE *, int,
+                                       TYPE, TYPE *, int) {
     ERROR("Not finished");
 }
 

common/MPI.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // This file impliments a wrapper class for MPI functions
 
 #include "common/MPI.h"
@@ -3691,28 +3707,29 @@ MPI MPI::loadBalance(double local, std::vector<double> work) {
     return split(0, key[getRank()]);
 }
 
-
-
 /****************************************************************************
  * Function Persistent Communication                                         *
  ****************************************************************************/
 template <>
-std::shared_ptr<MPI_Request> MPI::Isend_init<double>(const double *buf, int N, int proc, int tag) const
-{
-    std::shared_ptr<MPI_Request> obj( new MPI_Request, []( MPI_Request *req ) { MPI_Request_free( req ); delete req; } );
-    MPI_Send_init( buf, N, MPI_DOUBLE, proc, tag, communicator, obj.get() );
+std::shared_ptr<MPI_Request> MPI::Isend_init<double>(const double *buf, int N,
+                                                     int proc, int tag) const {
+    std::shared_ptr<MPI_Request> obj(new MPI_Request, [](MPI_Request *req) {
+        MPI_Request_free(req);
+        delete req;
+    });
+    MPI_Send_init(buf, N, MPI_DOUBLE, proc, tag, communicator, obj.get());
     return obj;
 }
-template<>
-std::shared_ptr<MPI_Request> MPI::Irecv_init<double>(double *buf, int N, int proc, int tag) const
-{
-    std::shared_ptr<MPI_Request> obj( new MPI_Request, []( MPI_Request *req ) { MPI_Request_free( req ); delete req; } );
-    MPI_Recv_init( buf, N, MPI_DOUBLE, proc, tag, communicator, obj.get() );
+template <>
+std::shared_ptr<MPI_Request> MPI::Irecv_init<double>(double *buf, int N,
+                                                     int proc, int tag) const {
+    std::shared_ptr<MPI_Request> obj(new MPI_Request, [](MPI_Request *req) {
+        MPI_Request_free(req);
+        delete req;
+    });
+    MPI_Recv_init(buf, N, MPI_DOUBLE, proc, tag, communicator, obj.get());
     return obj;
 }
-void MPI::Start( MPI_Request &request )
-{
-    MPI_Start( &request );
-}
+void MPI::Start(MPI_Request &request) { MPI_Start(&request); }
 
 } // namespace Utilities

common/MPI.h

@@ -1,24 +1,23 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // This file includes a wrapper class for MPI functions
 // 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
 
@@ -745,8 +744,8 @@ public: // Member functions
      *                  need to manually free the request
      */
     template <class type>
-    std::shared_ptr<MPI_Request> Isend_init(const type *buf, int length, int recv_proc,
-                           int tag) const;
+    std::shared_ptr<MPI_Request> Isend_init(const type *buf, int length,
+                                            int recv_proc, int tag) const;
 
     /*!
      * @brief This function sets up an Irecv call (see MPI_Recv_init)
@@ -759,13 +758,14 @@ public: // Member functions
      *                  need to manually free the request
      */
     template <class type>
-    std::shared_ptr<MPI_Request> Irecv_init(type *buf, int length, int send_proc, int tag) const;
+    std::shared_ptr<MPI_Request> Irecv_init(type *buf, int length,
+                                            int send_proc, int tag) const;
 
     /*!
      * @brief Start the MPI communication
      * @param request   Request to start
      */
-    void Start( MPI_Request &request );
+    void Start(MPI_Request &request);
 
     /*!
      * Each processor sends every other processor a single value.

common/Membrane.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* Flow adaptor class for multiphase flow methods */
 
 #ifndef ScaLBL_Membrane_INC
@@ -21,9 +37,9 @@
 * @param dist - memory buffer to hold the distributions
 * @param N - size of the distributions (derived from Domain structure)
 */
-extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,  
-		int *d3q7_recvlist, double *recvbuf, int count,
-		double *dist, int N,  double *coef);
+extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, int *d3q7_recvlist,
+                                            double *recvbuf, int count,
+                                            double *dist, int N, double *coef);
 
 /**
 * \brief Set custom link rules for D3Q19 distribution based on membrane location
@@ -38,8 +54,10 @@ extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
 * @param dist - memory buffer to hold the distributions
 * @param N - size of the distributions (derived from Domain structure)
 */
-extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *coef, int start, int offset, 
-		int linkCount, double *recvbuf, double *dist, int N);
+extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links,
+                                         double *coef, int start, int offset,
+                                         int linkCount, double *recvbuf,
+                                         double *dist, int N);
 
 /**
  * \class Membrane
@@ -51,33 +69,35 @@ extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *c
 class Membrane {
 public:
     int Np;
-    int Nx,Ny,Nz,N;
+    int Nx, Ny, Nz, N;
     int membraneLinkCount;
-    
-    int *initialNeighborList;   // original neighborlist
-    int *NeighborList;			// modified neighborlist
+    int BoundaryCondition;
+
+    int *initialNeighborList; // original neighborlist
+    int *NeighborList;        // modified neighborlist
 
     /* host data structures */
-    int *membraneLinks;    // D3Q7 links that cross membrane
-    int *membraneTag;      // label each link in the membrane
-    double *membraneDist;  // distance to membrane for each linked site
-    double *membraneOrientation;  // distance to membrane for each linked site
+    int *membraneLinks;          // D3Q7 links that cross membrane
+    int *membraneTag;            // label each link in the membrane
+    double *membraneDist;        // distance to membrane for each linked site
+    double *membraneOrientation; // distance to membrane for each linked site
 
     /* 
      * Device data structures
      */
     int *MembraneLinks;
-    double *MembraneCoef;  // mass transport coefficient for the membrane 
+    double *MembraneCoef; // mass transport coefficient for the membrane
     double *MembraneDistance;
     double *MembraneOrientation;
-    
+
     /**
     * \brief Create a flow adaptor to operate on the LB model
     * @param         ScaLBL - originating data structures 
     * @param 		 neighborList - list of neighbors for each site
     */
     //Membrane(std::shared_ptr <Domain> Dm, int *initialNeighborList, int Nsites);
-    Membrane(std::shared_ptr <ScaLBL_Communicator> sComm, int *dvcNeighborList, int Nsites);
+    Membrane(std::shared_ptr<ScaLBL_Communicator> sComm, int *dvcNeighborList,
+             int Nsites);
 
     /**
 	* \brief Destructor
@@ -92,40 +112,44 @@ public:
     * @param Map            - mapping between regular layout and compact layout
     */
     int Create(DoubleArray &Distance, IntArray &Map);
-    
+
     /**
      * \brief Write membrane data to output file
      * @param filename  - name of file to save
     */
     void Write(string filename);
-    
+
     /**
      * \brief Read membrane data from input file
      * @param filename  - name of file to save
     */
     void Read(string filename);
-        
-	void SendD3Q7AA(double *dist);
-	void RecvD3Q7AA(double *dist);
-	void AssignCoefficients(int *Map, double *Psi, double Threshold,
-			double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn,
-			double ThresholdMassFractionOut);
-	void IonTransport(double *dist, double *den);
-	//......................................................................................
-	// Buffers to store data sent and recieved by this MPI process
-	double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z;
-	double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z;
-	//......................................................................................
-    
+
+    void SendD3Q7AA(double *dist);
+    void RecvD3Q7AA(double *dist, bool BounceBack);
+    void AssignCoefficients(int *Map, double *Psi, double Threshold,
+                            double MassFractionIn, double MassFractionOut,
+                            double ThresholdMassFractionIn,
+                            double ThresholdMassFractionOut);
+    void IonTransport(double *dist, double *den);
+    //......................................................................................
+    // Buffers to store data sent and recieved by this MPI process
+    double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y,
+        *sendbuf_Z;
+    double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y,
+        *recvbuf_Z;
+    //......................................................................................
+
 private:
-	bool Lock; 	// use Lock to make sure only one call at a time to protect data in transit
-	int sendtag, recvtag;
-    int iproc,jproc,kproc;
-    int nprocx,nprocy,nprocz;
-	// Give the object it's own MPI communicator
-	RankInfoStruct rank_info;
-	Utilities::MPI MPI_COMM_SCALBL;		// MPI Communicator for this domain
-	MPI_Request req1[18],req2[18];
+    bool
+        Lock; // use Lock to make sure only one call at a time to protect data in transit
+    int sendtag, recvtag;
+    int iproc, jproc, kproc;
+    int nprocx, nprocy, nprocz;
+    // Give the object it's own MPI communicator
+    RankInfoStruct rank_info;
+    Utilities::MPI MPI_COMM_SCALBL; // MPI Communicator for this domain
+    MPI_Request req1[18], req2[18];
     /**
     * \brief   Set up membrane communication
     * \details associate p2p communication links to membrane where necessary
@@ -142,49 +166,59 @@ private:
     * @param dvcMap      	- data structure used to define mapping between dense and sparse representation
     * @param Np      		- number of sites in dense representation
     * */
-	int D3Q7_MapRecv(int Cqx, int Cqy, int Cqz, int *d3q19_recvlist, 
-			int count, int *membraneRecvLabels, DoubleArray &Distance,  int  *dvcMap);
-	//......................................................................................
-	// MPI ranks for all 18 neighbors
-	//......................................................................................
-	// These variables are all private to prevent external things from modifying them!!
-	//......................................................................................
-	int rank;
-	int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
-	int rank_xy,rank_XY,rank_xY,rank_Xy;
-	int rank_xz,rank_XZ,rank_xZ,rank_Xz;
-	int rank_yz,rank_YZ,rank_yZ,rank_Yz;
-	//......................................................................................
-	int SendCount, RecvCount, CommunicationCount;
-	//......................................................................................
-	int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
-	//......................................................................................
-	int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
-	//......................................................................................
-	int linkCount_x[5], linkCount_y[5], linkCount_z[5], linkCount_X[5], linkCount_Y[5], linkCount_Z[5];
-	int linkCount_xy, linkCount_yz, linkCount_xz, linkCount_Xy, linkCount_Yz, linkCount_xZ;
-	int linkCount_xY, linkCount_yZ, linkCount_Xz, linkCount_XY, linkCount_YZ, linkCount_XZ;
-	//......................................................................................
-	// Send buffers that reside on the compute device
-	int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X, *dvcSendList_Y, *dvcSendList_Z;
-	//int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ;
-	//int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ;
-	// Recieve buffers that reside on the compute device
-	int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X, *dvcRecvList_Y, *dvcRecvList_Z;
-	//int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ;
-	//int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ;
-	// Link lists  that reside on the compute device
-	int *dvcRecvLinks_x, *dvcRecvLinks_y, *dvcRecvLinks_z, *dvcRecvLinks_X, *dvcRecvLinks_Y, *dvcRecvLinks_Z;
-	//int *dvcRecvLinks_xy, *dvcRecvLinks_yz, *dvcRecvLinks_xz, *dvcRecvLinks_Xy, *dvcRecvLinks_Yz, *dvcRecvLinks_xZ;
-	//int *dvcRecvLinks_xY, *dvcRecvLinks_yZ, *dvcRecvLinks_Xz, *dvcRecvLinks_XY, *dvcRecvLinks_YZ, *dvcRecvLinks_XZ;
-	// Recieve buffers for the distributions
-	int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z;
-	//int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ;
-	//int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ;
-	//......................................................................................
-	// mass transfer coefficient arrays
-	double *coefficient_x, *coefficient_X, *coefficient_y, *coefficient_Y, *coefficient_z, *coefficient_Z;
-	//......................................................................................
-
+    int D3Q7_MapRecv(int Cqx, int Cqy, int Cqz, int *d3q19_recvlist, int count,
+                     int *membraneRecvLabels, DoubleArray &Distance,
+                     int *dvcMap);
+    //......................................................................................
+    // MPI ranks for all 18 neighbors
+    //......................................................................................
+    // These variables are all private to prevent external things from modifying them!!
+    //......................................................................................
+    int rank;
+    int rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z;
+    int rank_xy, rank_XY, rank_xY, rank_Xy;
+    int rank_xz, rank_XZ, rank_xZ, rank_Xz;
+    int rank_yz, rank_YZ, rank_yZ, rank_Yz;
+    //......................................................................................
+    int SendCount, RecvCount, CommunicationCount;
+    //......................................................................................
+    int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y,
+        sendCount_Z;
+    //......................................................................................
+    int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y,
+        recvCount_Z;
+    //......................................................................................
+    int linkCount_x[5], linkCount_y[5], linkCount_z[5], linkCount_X[5],
+        linkCount_Y[5], linkCount_Z[5];
+    int linkCount_xy, linkCount_yz, linkCount_xz, linkCount_Xy, linkCount_Yz,
+        linkCount_xZ;
+    int linkCount_xY, linkCount_yZ, linkCount_Xz, linkCount_XY, linkCount_YZ,
+        linkCount_XZ;
+    //......................................................................................
+    // Send buffers that reside on the compute device
+    int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X,
+        *dvcSendList_Y, *dvcSendList_Z;
+    //int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ;
+    //int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ;
+    // Recieve buffers that reside on the compute device
+    int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X,
+        *dvcRecvList_Y, *dvcRecvList_Z;
+    //int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ;
+    //int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ;
+    // Link lists  that reside on the compute device
+    int *dvcRecvLinks_x, *dvcRecvLinks_y, *dvcRecvLinks_z, *dvcRecvLinks_X,
+        *dvcRecvLinks_Y, *dvcRecvLinks_Z;
+    //int *dvcRecvLinks_xy, *dvcRecvLinks_yz, *dvcRecvLinks_xz, *dvcRecvLinks_Xy, *dvcRecvLinks_Yz, *dvcRecvLinks_xZ;
+    //int *dvcRecvLinks_xY, *dvcRecvLinks_yZ, *dvcRecvLinks_Xz, *dvcRecvLinks_XY, *dvcRecvLinks_YZ, *dvcRecvLinks_XZ;
+    // Recieve buffers for the distributions
+    int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X,
+        *dvcRecvDist_Y, *dvcRecvDist_Z;
+    //int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ;
+    //int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ;
+    //......................................................................................
+    // mass transfer coefficient arrays
+    double *coefficient_x, *coefficient_X, *coefficient_y, *coefficient_Y,
+        *coefficient_z, *coefficient_Z;
+    //......................................................................................
 };
 #endif

common/ReadMicroCT.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include "common/ReadMicroCT.h"
 #include "common/Utilities.h"
 

common/ReadMicroCT.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #ifndef READMICROCT_H
 #define READMICROCT_H
 

common/SpherePack.cpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #include <stdio.h>
 #include <stdlib.h>
 #include <iostream>

common/Utilities.hpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #ifndef included_Utilities_hpp
 #define included_Utilities_hpp
 
@@ -47,7 +63,7 @@ template <class T> void quicksort(std::vector<T> &x) {
         } else {
             k = (l + ir) /
                 2; // Choose median of left, center and right elements as partitioning
-                // element a. Also rearrange so that a(l) < a(l+1) < a(ir).
+            // element a. Also rearrange so that a(l) < a(l+1) < a(ir).
             tmp_a = arr[k];
             arr[k] = arr[l + 1];
             arr[l + 1] = tmp_a;
@@ -140,7 +156,7 @@ void quicksort(std::vector<T1> &x, std::vector<T2> &y) {
         } else {
             k = (l + ir) /
                 2; // Choose median of left, center and right elements as partitioning
-                // element a. Also rearrange so that a(l) ? a(l+1) ? a(ir).
+            // element a. Also rearrange so that a(l) ? a(l+1) ? a(ir).
             tmp_a = arr[k];
             arr[k] = arr[l + 1];
             arr[l + 1] = tmp_a;

common/WideHalo.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* 
 This class implements support for halo widths larger than 1
  */

common/WideHalo.h

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* 
 This class implements support for halo widths larger than 1
  */

cpu/D3Q19.cpp

@@ -48,7 +48,6 @@ extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count,
     }
 }
 
-
 extern "C" void ScaLBL_D3Q19_AA_Init(double *f_even, double *f_odd, int Np) {
     int n;
     for (n = 0; n < Np; n++) {
@@ -1942,8 +1941,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_Compact(double *dist, int Np) {
     }
 }
 
-extern "C" void ScaLBL_D3Q19_AAodd_Compact(int *neighborList,
-                                           double *dist, int Np) {
+extern "C" void ScaLBL_D3Q19_AAodd_Compact(int *neighborList, double *dist,
+                                           int Np) {
     int nread;
 
     for (int n = 0; n < Np; n++) {

cpu/D3Q7BC.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // CPU Functions for D3Q7 Lattice Boltzmann Methods
 // Boundary Conditions
 
@@ -35,31 +51,39 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue,
     }
 }
 
-extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist,double *BoundaryValue,int* BoundaryLabel,int *BounceBackDist_list,int *BounceBackSolid_list,int N){
+extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(
+    double *dist, double *BoundaryValue, int *BoundaryLabel,
+    int *BounceBackDist_list, int *BounceBackSolid_list, int N) {
 
     int idx;
-    int iq,ib;
-    double value_b,value_b_label,value_q;
-	for (idx=0; idx<N; idx++){
-		iq = BounceBackDist_list[idx];
+    int iq, ib;
+    double value_b, value_b_label, value_q;
+    for (idx = 0; idx < N; idx++) {
+        iq = BounceBackDist_list[idx];
         ib = BounceBackSolid_list[idx];
-		value_b = BoundaryValue[ib];//get boundary value from a solid site
-		value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site
+        value_b = BoundaryValue[ib]; //get boundary value from a solid site
+        value_b_label = BoundaryLabel
+            [ib]; //get boundary label (i.e. type of BC) from a solid site
         value_q = dist[iq];
-        if (value_b_label==1){//Dirichlet BC
-		    dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice
+        if (value_b_label == 1) { //Dirichlet BC
+            dist[iq] =
+                -1.0 * value_q +
+                value_b *
+                    0.25; //NOTE 0.25 is the speed of sound for D3Q7 lattice
         }
-        if (value_b_label==2){//Neumann BC
-		    dist[iq] = value_q + value_b;
+        if (value_b_label == 2) { //Neumann BC
+            dist[iq] = value_q + value_b;
         }
-	}
+    }
 }
 
-extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
-                                                      double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
-                                                      int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
-                                                      double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
-                                                      int count, int Np){
+extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(
+    double *dist, double *zeta_potential, double *ElectricField,
+    double *SolidGrad, double epsilon_LB, double tau, double rho0,
+    double den_scale, double h, double time_conv, int *BounceBackDist_list,
+    int *BounceBackSolid_list, int *FluidBoundary_list, double *lattice_weight,
+    float *lattice_cx, float *lattice_cy, float *lattice_cz, int count,
+    int Np) {
 
     int idx;
     int iq, ib, ifluidBC;
@@ -181,7 +205,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList,
     }
 }
 
-
 extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList,
                                                          int *list,
                                                          double *dist,
@@ -214,8 +237,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList,
     }
 }
 
-
-
 extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi,
                                          double Vin, int count) {
     int idx, n, nm;

cpu/Extras.cpp

@@ -1,5 +1,6 @@
 /*
   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

cpu/FreeLee.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 #include <stdio.h>
 

cpu/Greyscale.cpp

@@ -1,17 +1,18 @@
 /*
-  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 <http://www.gnu.org/licenses/>.       
+  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 <http://www.gnu.org/licenses/>.
 */
 #include <math.h>
 

cpu/GreyscaleColor.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 
 extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(

cpu/Ion.cpp

@@ -1,43 +1,307 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <math.h>
 
-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){
-	
-	int link,iq,ip,nq,np,nqm,npm;
-	double aq, ap, membranePotential;
-	//double dq, dp, dist, orientation;
-	/* Interior Links */
-	for (link=0; link<memLinks; link++){
-		
-		// inside             	//outside
-		aq = MassFractionIn;	ap = MassFractionOut;  
-		iq = membrane[2*link]; 	ip = membrane[2*link+1];
-		nq = iq%Np;				np = ip%Np;
-		nqm = Map[nq];			npm = Map[np]; // strided layout
-		//dq = Distance[nqm];     dp = Distance[npm];
-		/* orientation for link to distance gradient*/
-		//orientation = 1.0/fabs(dq - dp);
-		
-		/* membrane potential for this link */
-		membranePotential = Psi[nqm] - Psi[npm];
-		if (membranePotential > Threshold){
-			aq = ThresholdMassFractionIn;	ap = ThresholdMassFractionOut;  
-		}
-		
-		/* Save the mass transfer coefficients */
-		//coef[2*link] = aq*orientation;		coef[2*link+1] = ap*orientation;
-		coef[2*link] = aq;		coef[2*link+1] = ap;
-	}
+/***** 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) {
+
+    int link, iq, ip, nq, np, nqm, npm;
+    double aq, ap, membranePotential;
+    //double dq, dp, dist, orientation;
+    /* Interior Links */
+    for (link = 0; link < memLinks; link++) {
+
+        // inside             	//outside
+        aq = MassFractionIn;
+        ap = MassFractionOut;
+        iq = membrane[2 * link];
+        ip = membrane[2 * link + 1];
+        nq = iq % Np;
+        np = ip % Np;
+        nqm = Map[nq];
+        npm = Map[np]; // strided layout
+        //dq = Distance[nqm];     dp = Distance[npm];
+        /* orientation for link to distance gradient*/
+        //orientation = 1.0/fabs(dq - dp);
+
+        /* membrane potential for this link */
+        membranePotential = Psi[nqm] - Psi[npm];
+        if (membranePotential > Threshold) {
+            aq = ThresholdMassFractionIn;
+            ap = ThresholdMassFractionOut;
+        }
+
+        /* Save the mass transfer coefficients */
+        //coef[2*link] = aq*orientation;		coef[2*link+1] = ap*orientation;
+        coef[2 * link] = aq;
+        coef[2 * link + 1] = ap;
+    }
 }
 
 extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
-		const int Cqx, const int Cqy, int const Cqz, 
-		int *Map, double *Distance, double *Psi, double Threshold, 
-		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) {
+    const int Cqx, const int Cqy, int const Cqz, int *Map, double *Distance,
+    double *Psi, double Threshold, 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
@@ -45,58 +309,58 @@ extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
     // dist may be even or odd distributions stored by stream layout
     //....................................................................................
     int n, idx, label, nqm, npm, i, j, k;
-    double distanceLocal;//, distanceNonlocal;
+    double distanceLocal; //, distanceNonlocal;
     double psiLocal, psiNonlocal, membranePotential;
-    double ap,aq; // coefficient
+    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];
+        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;
     }
 }
 
-
-extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,  
-		int *d3q7_recvlist, double *recvbuf, int count,
-		double *dist, int N,  double *coef) {
+extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q, int *d3q7_recvlist,
+                                            double *recvbuf, int count,
+                                            double *dist, int N, double *coef) {
     //....................................................................................
     // Unack distribution from the recv buffer
     // Distribution q matche Cqx, Cqy, Cqz
@@ -104,40 +368,50 @@ extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
     // dist may be even or odd distributions stored by stream layout
     //....................................................................................
     int n, idx;
-    double fq,fp,fqq,ap,aq; // coefficient
+    double fq, fp, fqq, ap, aq; // coefficient
     /* First unpack the regular links */
     for (idx = 0; idx < count; idx++) {
-    	n = d3q7_recvlist[idx];
+        n = d3q7_recvlist[idx];
         // update link based on mass transfer coefficients
-        if (!(n < 0)){
-        	aq = coef[2*idx];
-        	ap = coef[2*idx+1];
-        	fq = dist[q * N + n];
-        	fp = recvbuf[idx];
-        	fqq = (1-aq)*fq+ap*fp;
+        if (!(n < 0)) {
+            aq = coef[2 * idx];
+            ap = coef[2 * idx + 1];
+            fq = dist[q * N + n];
+            fp = recvbuf[idx];
+            fqq = (1 - aq) * fq + ap * fp;
             dist[q * N + n] = fqq;
         }
         //printf(" LINK: site=%i, index=%i \n", n, idx);
-    } 
+    }
 }
 
-extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, 
-		double *dist, double *Den, int memLinks, int Np){
-	
-	int link,iq,ip,nq,np;
-	double aq, ap, fq, fp, fqq, fpp, Cq, Cp;
-	for (link=0; link<memLinks; link++){
-		// inside             	//outside
-		aq = coef[2*link];		ap = coef[2*link+1];
-		iq = membrane[2*link]; 	ip = membrane[2*link+1];
-		nq = iq%Np;				np = ip%Np;
-		fq  = dist[iq];			fp = dist[ip];
-		fqq = (1-aq)*fq+ap*fp;	fpp = (1-ap)*fp+aq*fq;
-		Cq = Den[nq];			Cp = Den[np];
-		Cq += fqq - fq;			Cp += fpp - fp;
-		Den[nq] = Cq;			Den[np] = Cp;
-		dist[iq] = fqq;			dist[ip] = fpp;
-	}
+extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
+                                                  double *dist, double *Den,
+                                                  int memLinks, int Np) {
+
+    int link, iq, ip, nq, np;
+    double aq, ap, fq, fp, fqq, fpp, Cq, Cp;
+    for (link = 0; link < memLinks; link++) {
+        // inside             	//outside
+        aq = coef[2 * link];
+        ap = coef[2 * link + 1];
+        iq = membrane[2 * link];
+        ip = membrane[2 * link + 1];
+        nq = iq % Np;
+        np = ip % Np;
+        fq = dist[iq];
+        fp = dist[ip];
+        fqq = (1 - aq) * fq + ap * fp;
+        fpp = (1 - ap) * fp + aq * fq;
+        Cq = Den[nq];
+        Cp = Den[np];
+        Cq += fqq - fq;
+        Cp += fpp - fp;
+        Den[nq] = Cq;
+        Den[np] = Cp;
+        dist[iq] = fqq;
+        dist[ip] = fpp;
+    }
 }
 
 extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList,
@@ -225,13 +499,12 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den,
     }
 }
 
-extern "C" void ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist,
-                                      double *Den, double *FluxDiffusive,
-                                      double *FluxAdvective,
-                                      double *FluxElectrical, double *Velocity,
-                                      double *ElectricField, double Di, int zi,
-                                      double rlx, double Vt, int start,
-                                      int finish, int Np) {
+extern "C" void
+ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist, double *Den,
+                         double *FluxDiffusive, double *FluxAdvective,
+                         double *FluxElectrical, double *Velocity,
+                         double *ElectricField, double Di, int zi, double rlx,
+                         double Vt, int start, int finish, int Np) {
     int n;
     double Ci;
     double ux, uy, uz;
@@ -291,7 +564,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist,
         FluxElectrical[n + 0 * Np] = uEPx * Ci;
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         FluxElectrical[n + 2 * Np] = uEPz * Ci;
-        
+
         //Den[n] = Ci;
 
         /* use logistic function to prevent negative distributions*/
@@ -307,35 +580,33 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist,
 
         // 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 + 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 - 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 + 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 - 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 + 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 - 4.0 * (uz + uEPz));
         //    f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
-
     }
 }
 
@@ -350,7 +621,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++) {
 
@@ -388,9 +658,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion_v0(
         FluxElectrical[n + 0 * Np] = uEPx * Ci;
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         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);
@@ -404,32 +674,32 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion_v0(
 
         // 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 + 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 - 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 + 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 - 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 + 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 - 4.0 * (uz + uEPz));
         //        f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
     }
 }
@@ -500,51 +770,35 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
         FluxElectrical[n + 0 * Np] = uEPx * Ci;
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         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);
+        Den[n] = Ci;
 
         // 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);
-
+            f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
         // 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);
+            f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
         // 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 );
+            f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
         // 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);
+            f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
         // 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);
+            f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
         // 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);
-
+            f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
     }
 }
 
@@ -559,7 +813,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++) {
 
@@ -597,49 +850,35 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
         FluxElectrical[n + 0 * Np] = uEPx * Ci;
         FluxElectrical[n + 1 * Np] = uEPy * Ci;
         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;
 
         // 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);
+            f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
         // 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);
+            f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
         // 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);
+            f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
         // 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);
+            f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
         // 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);
+            f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
         // 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);
+            f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
     }
 }
 
@@ -676,8 +915,9 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den,
 
 extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den,
                                               double *ChargeDensity,
-                                              double IonValence, int ion_component,
-                                              int start, int finish, int Np) {
+                                              double IonValence,
+                                              int ion_component, int start,
+                                              int finish, int Np) {
 
     int n;
     double Ci; //ion concentration of species i

cpu/MRT.cpp

@@ -1,5 +1,6 @@
 /*
   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
@@ -13,7 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-
 extern "C" void INITIALIZE(char *ID, double *f_even, double *f_odd, int Nx,
                            int Ny, int Nz) {
     int n, N;

cpu/MembraneHelper.cpp

@@ -1,6 +1,23 @@
+/*
+  Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University
+  Copyright Equnior ASA
 
-extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int linkCount,
-                                    double *recvbuf, double *dist, int N) {
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start,
+                                      int linkCount, double *recvbuf,
+                                      double *dist, int N) {
     //....................................................................................
     // Unack distribution from the recv buffer
     // Distribution q matche Cqx, Cqy, Cqz
@@ -8,9 +25,9 @@ extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, i
     // dist may be even or odd distributions stored by stream layout
     //....................................................................................
     int n, idx, link;
-    for (link=0; link<linkCount; link++){
+    for (link = 0; link < linkCount; link++) {
 
-    	idx = links[start+link];
+        idx = links[start + link];
         // Get the value from the list -- note that n is the index is from the send (non-local) process
         n = list[start + idx];
         // unpack the distribution to the proper location
@@ -19,8 +36,10 @@ extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, i
     }
 }
 
-extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *coef, int start, int offset, 
-		int linkCount, double *recvbuf, double *dist, int N){
+extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links,
+                                         double *coef, int start, int offset,
+                                         int linkCount, double *recvbuf,
+                                         double *dist, int N) {
     //....................................................................................
     // Unack distribution from the recv buffer
     // Distribution q matche Cqx, Cqy, Cqz
@@ -29,14 +48,14 @@ extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *c
     //....................................................................................
     int n, idx, link;
     double alpha;
-    for (link=offset; link<linkCount; link++){
+    for (link = offset; link < linkCount; link++) {
 
-    	idx = list[start+link];
+        idx = list[start + link];
         // Get the value from the list -- note that n is the index is from the send (non-local) process
         n = list[start + idx];
         alpha = coef[start + idx];
         // unpack the distribution to the proper location
         if (!(n < 0))
-            dist[q * N + n] = alpha*recvbuf[start + idx];
+            dist[q * N + n] = alpha * recvbuf[start + idx];
     }
 }

cpu/MixedGradient.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* Implement Mixed Gradient (Lee et al. JCP 2016)*/
 
 extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi,

cpu/Stokes.cpp

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 
 extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
@@ -38,12 +54,20 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
         Ey = ElectricField[n + 1 * Np];
         Ez = ElectricField[n + 2 * Np];
         //compute total body force, including input body force (Gx,Gy,Gz)
-        Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale; //the extra factors at the end necessarily convert unit from phys to LB
-        Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale;
-        Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale;
+        Fx =
+            (UseSlippingVelBC == 1)
+                ? Gx
+                : Gx +
+                      rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+        Fy = (UseSlippingVelBC == 1)
+                 ? Gy
+                 : Gy + rhoE * Ey * (time_conv * time_conv) /
+                            (h * h * 1.0e-12) / den_scale;
+        Fz = (UseSlippingVelBC == 1)
+                 ? Gz
+                 : Gz + rhoE * Ez * (time_conv * time_conv) /
+                            (h * h * 1.0e-12) / den_scale;
 
         // q=0
         fq = dist[n];
@@ -520,13 +544,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
         //Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
         //              den_scale;
         //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
-        //and body force induced by external efectric field is reduced to slipping velocity BC. 
-        Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale; //the extra factors at the end necessarily convert unit from phys to LB
-        Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale;
-        Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale;
+        //and body force induced by external efectric field is reduced to slipping velocity BC.
+        Fx =
+            (UseSlippingVelBC == 1)
+                ? Gx
+                : Gx +
+                      rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+        Fy = (UseSlippingVelBC == 1)
+                 ? Gy
+                 : Gy + rhoE * Ey * (time_conv * time_conv) /
+                            (h * h * 1.0e-12) / den_scale;
+        Fz = (UseSlippingVelBC == 1)
+                 ? Gz
+                 : Gz + rhoE * Ez * (time_conv * time_conv) /
+                            (h * h * 1.0e-12) / den_scale;
 
         // q=0
         fq = dist[n];

cpu/dfh.cpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  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 <http://www.gnu.org/licenses/>.
-*/
 #include <math.h>
 #include <stdio.h>
 

cpu/thermal.cpp

@@ -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 <http://www.gnu.org/licenses/>.
-*/
-/*
-  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

cuda/BGK.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 
 #define NBLOCKS 1024

cuda/Color.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 #include <stdio.h>
 #include <cuda_profiler_api.h>

cuda/CudaExtras.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // Basic cuda functions callable from C/C++ code
 #include <cuda.h>
 

cuda/D3Q19.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <cooperative_groups.h>
 

cuda/D3Q7.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 // GPU Functions for D3Q7 Lattice Boltzmann Methods
 #include <stdio.h>
 

cuda/D3Q7BC.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 #include <stdio.h>
 #include <cuda_profiler_api.h>

cuda/Extras.cu

@@ -1,5 +1,6 @@
 /*
   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

cuda/FreeLee.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 #include <stdio.h>
 #include <cuda_profiler_api.h>

cuda/Greyscale.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 
 #define NBLOCKS 1024

cuda/GreyscaleColor.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <math.h>
 

cuda/Ion.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <math.h>
 //#include <cuda_profiler_api.h>
@@ -5,6 +21,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<S; s++){
+    	//........Get 1-D index for this thread....................
+    	n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+    	if (n<finish) {
+
+    		//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];
+
+    		// 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; s<S; s++){
+           //........Get 1-D index for this thread....................
+           n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
+           if (n<finish) {
+		
+        //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;     
+
+        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;
+
+        }
+    }
+}
+/**** end of pH equlibrium model ********/
+
 extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int linkCount,
                                     double *recvbuf, double *dist, int N) {
     //....................................................................................
@@ -299,23 +527,25 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *D
 
 __global__  void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField, 
                                            double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
-	int n;
-	double Ci;
-    double ux,uy,uz;
-    double uEPx,uEPy,uEPz;//electrochemical induced velocity
-    double Ex,Ey,Ez;//electrical field
-    double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
-	double f0,f1,f2,f3,f4,f5,f6;
-	double X,Y,Z,factor_x,factor_y,factor_z;
-	int nr1,nr2,nr3,nr4,nr5,nr6;
+    int n;
+    double Ci;
+    double ux, uy, uz;
+    double uEPx, uEPy, uEPz; //electrochemical induced velocity
+    double Ex, Ey, Ez;       //electrical field
+    double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
+    double f0, f1, f2, f3, f4, f5, f6;
+    //double X,Y,Z,factor_x, factor_y, factor_z;
+    int nr1, nr2, nr3, nr4, nr5, nr6;
 
 	int S = Np/NBLOCKS/NTHREADS + 1;
 	for (int s=0; s<S; s++){
 		//........Get 1-D index for this thread....................
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
+	        //Load data
 
 	        //Load data
+	        //Ci = Den[n];
 	        Ex = ElectricField[n + 0 * Np];
 	        Ey = ElectricField[n + 1 * Np];
 	        Ez = ElectricField[n + 2 * Np];
@@ -364,46 +594,32 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
 	        
 	        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 + factor_x);
-	        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
-
+	        f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
 	        // q=2
 	        dist[nr1] =
-	                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-	        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+	        f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
 	        // q = 3
 	        dist[nr4] =
-	                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
-	        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+	        f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
 	        // q = 4
 	        dist[nr3] =
-	                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-	        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+	        f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
 	        // q = 5
 	        dist[nr6] =
-	                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 +  factor_z);
-	        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+	        f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
 	        // q = 6
 	        dist[nr5] =
-	            f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
+	        f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
 
 		}
 	}
@@ -411,14 +627,13 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
 
 __global__  void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField, 
                                             double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
-	int n;
-	double Ci;
-    double ux,uy,uz;
-    double uEPx,uEPy,uEPz;//electrochemical induced velocity
-    double Ex,Ey,Ez;//electrical field
-    double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
-	double f0,f1,f2,f3,f4,f5,f6;
-	double X,Y,Z,factor_x,factor_y,factor_z;
+    int n;
+    double Ci;
+    double ux, uy, uz;
+    double uEPx, uEPy, uEPz; //electrochemical induced velocity
+    double Ex, Ey, Ez;       //electrical field
+    double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
+    double f0, f1, f2, f3, f4, f5, f6;
 
 	int S = Np/NBLOCKS/NTHREADS + 1;
 	for (int s=0; s<S; s++){
@@ -426,83 +641,69 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
 		n =  S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
 		if (n<finish) {
 
-		       //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 = zi * Di / Vt * Ex;
-		        uEPy = zi * Di / Vt * Ey;
-		        uEPz = zi * Di / Vt * Ez;
+	        //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 = zi * Di / Vt * Ex;
+	        uEPy = zi * Di / Vt * Ey;
+	        uEPz = zi * Di / Vt * Ez;
 
-		        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];
+	        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];
 
-		        // compute diffusive flux
-		        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
-		        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
-		        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
-		        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
-		        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
-		        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
-		        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
-		        FluxAdvective[n + 0 * Np] = ux * Ci;
-		        FluxAdvective[n + 1 * Np] = uy * Ci;
-		        FluxAdvective[n + 2 * Np] = uz * Ci;
-		        FluxElectrical[n + 0 * Np] = uEPx * Ci;
-		        FluxElectrical[n + 1 * Np] = uEPy * Ci;
-		        FluxElectrical[n + 2 * Np] = uEPz * Ci;
-		        
-		        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);
+	        // compute diffusive flux
+	        Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
+	        flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
+	        flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
+	        flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
+	        FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
+	        FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
+	        FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
+	        FluxAdvective[n + 0 * Np] = ux * Ci;
+	        FluxAdvective[n + 1 * Np] = uy * Ci;
+	        FluxAdvective[n + 2 * Np] = uz * Ci;
+	        FluxElectrical[n + 0 * Np] = uEPx * Ci;
+	        FluxElectrical[n + 1 * Np] = uEPy * Ci;
+	        FluxElectrical[n + 2 * Np] = uEPz * Ci;
+	        
+	        Den[n] = Ci;
 
-		        // q=0
-		        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
+	        // q=0
+	        dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
 
-		        // q = 1
-		        dist[1 * Np + n] =
-		                f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
-		        //f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
+	        // q = 1
+	        dist[1 * Np + n] =
+	        f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
 
-		        // q=2
-		        dist[2 * Np + n] =
-		                f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
-		        //f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
+	        // q=2
+	        dist[2 * Np + n] =
+	        f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
 
-		        // q = 3
-		        dist[3 * Np + n] =
-		                f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
-		        //f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
+	        // q = 3
+	        dist[3 * Np + n] =
+	        f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
 
-		        // q = 4
-		        dist[4 * Np + n] =
-		                f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
-		        //f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
+	        // q = 4
+	        dist[4 * Np + n] =
+	        f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
 
-		        // q = 5
-		        dist[5 * Np + n] =
-		                f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
-		        //f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
+	        // q = 5
+	        dist[5 * Np + n] =
+	        f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
 
-		        // q = 6
-		        dist[6 * Np + n] =
-		                f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
-		        //f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
+	        // q = 6
+	        dist[6 * Np + n] =
+	        f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
 		}
 	}
 }
@@ -568,10 +769,6 @@ __global__  void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDe
             CD_tmp = F*IonValence*Ci;
             ChargeDensity[n] = CD + CD_tmp;
 
- //           Ci = Den[n+ion_component*Np];
-   //         CD = ChargeDensity[n];
-     //       CD_tmp = F*IonValence*Ci;
-       //     ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
 		}
 	}
 }
@@ -980,3 +1177,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<<<NBLOCKS,NTHREADS >>>(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<<<NBLOCKS,NTHREADS >>>(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));
+	}
+
+}
+

cuda/MRT.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 //*************************************************************************
 // CUDA kernels for single-phase ScaLBL_D3Q19_MRT code
 // James McClure

cuda/MixedGradient.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 /* Implement Mixed Gradient (Lee et al. JCP 2016)*/
 #include <cuda.h>
 #include <stdio.h>

cuda/Poisson.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <math.h>
 //#include <cuda_profiler_api.h>
@@ -522,7 +538,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 +698,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 +785,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<<<NBLOCKS,NTHREADS >>>(neighborList, Map,
@@ -783,8 +800,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) {
 

cuda/Stokes.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <stdio.h>
 #include <math.h>
 //#include <cuda_profiler_api.h>

cuda/dfh.cu

@@ -1,3 +1,19 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
 #include <math.h>
 #include <stdio.h>
 #include <cuda_profiler_api.h>

cuda/exe/CMakeLists.txt

@@ -1,11 +0,0 @@
-#INSTALL_LBPM_EXE( lb1_MRT_mpi )
-INSTALL_LBPM_EXE( lb2_Color )
-INSTALL_LBPM_EXE( lb2_Color_mpi )
-#INSTALL_LBPM_EXE( lb2_Color_pBC_wia_mpi )
-INSTALL_LBPM_EXE( lb2_Color_wia_mpi )
-
-
-# Run the serial ConstrainedBubble inputs as a weekly test
-CONFIGURE_FILE( ${LBPM_SOURCE_DIR}/example/ConstrainedBubble/Color.in ${CMAKE_CURRENT_BINARY_DIR}/Color.in COPYONLY )
-CONFIGURE_FILE( ${LBPM_SOURCE_DIR}/example/ConstrainedBubble/Domain.in ${CMAKE_CURRENT_BINARY_DIR}/Domain.in COPYONLY )
-ADD_LBPM_WEEKLY_TEST( lb2_Color_wia_mpi 1 )

cuda/exe/lb1_MRT-swap.cu

@@ -1,248 +0,0 @@
-#include <stdio.h>
-#include <iostream>
-#include <fstream>
-#include <cuda.h>
-//#include <cutil.h>
-
-using namespace std;
-
-//*************************************************************************
-extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx,
-							  int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-extern "C" void dvc_SwapD3Q19(char *ID, double *f_even, double *f_odd, int Nx,
-							  int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-extern "C" void dvc_MRT(char *ID, double *f_even, double *f_odd, double rlxA, double rlxB, double Fx, double Fy, double Fz,
-		int Nx, int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-
-void Write_Out(double *array, int Nx, int Ny, int Nz){
-	int value;
-	FILE *output;
-	output = fopen("dist.list","w");
-	for (int k=0; k<Nz; k++){
-		for (int j=0; j<Ny; j++){
-			for (int i=0; i<Nx; i++){
-				int index = k*Nx*Ny+j*Nx+i;
-				value = int(array[index]);
-				fprintf(output, "| %i",value);
-			}
-			fprintf(output, " | \n");
-		}
-		fprintf(output,"************************************** \n");	
-	}
-	fclose(output);
-}
-
-//**************************************************************************
-// MRT implementation of the LBM using CUDA
-//**************************************************************************
-
-int main(void)
-{
-
-	int deviceCount;
-	cudaGetDeviceCount(&deviceCount);
-	int device = 1;
-	printf("Number of devices = %i \n", deviceCount);
-	printf("Current device is = %i \n", device);
-	cudaSetDevice(device);
-	
-	// BGK Model parameters
-	string FILENAME;	
-	unsigned int nBlocks, nthreads;
-	int timestepMax, interval;
-	double tau,Fx,Fy,Fz,tol;
-	// Domain variables
-	int Nx,Ny,Nz;
-
-	ifstream input("MRT.in");
-	input >> FILENAME;		// name of the input file
-	input >> Nz;			// number of nodes (x,y,z)
-	input >> nBlocks;				
-	input >> nthreads;				
-	input >> tau;				// relaxation time 
-	input >> Fx;			// External force components (x,y,z)
-	input >> Fy;
-	input >> Fz;
-	input >> timestepMax;			// max no. of timesteps
-	input >> interval;			// error interval
-	input >> tol;				// error tolerance
-	
-	double rlx_setA = 1.f/tau;
-	double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
-	
-	printf("tau = %f \n", tau);
-	printf("Set A = %f \n", rlx_setA);
-	printf("Set B = %f \n", rlx_setB);
-	printf("Force(x) = %f \n", Fx);
-	printf("Force(y) = %f \n", Fy);
-	printf("Force(z) = %f \n", Fz);
-
-	Nx = Ny = Nz;	// Cubic domain
-	
-	int N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
-	
-//	unsigned int nBlocks = 32;
-//	int nthreads = 128;
-	int S = N/nthreads/nBlocks;
-	
-//	unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1);
-	dim3 grid(nBlocks,1,1);
-		
-	printf("Number of blocks = %i \n", nBlocks);
-	printf("Threads per block = %i \n", nthreads);
-	printf("Sweeps per thread = %i \n", S);
-	printf("Number of nodes per side = %i \n", Nx);
-	printf("Total Number of nodes = %i \n", N);
-	
-	//.......................................................................
-	printf("Read input media... \n");
-	// .......... READ THE INPUT FILE .......................................
-	int n;
-	char value;
-	char *id;
-	id = new char[N];	
-	int sum = 0;
-	double porosity;
-	ifstream PM(FILENAME.c_str(),ios::binary);
-	for (int k=0;k<Nz;k++){
-		for (int j=0;j<Ny;j++){
-			for (int i=0;i<Nx;i++){
-				PM.read((char *) (&value), sizeof(value));
-				n = k*Nx*Ny+j*Nx+i;
-				id[n] = value;
-				if (value > 0) sum++;
-			}
-		}
-	}
-	PM.close();
-	printf("File porosity = %f\n", double(sum)/N);
-	//.......................................................................
-	//...........device phase ID.................................................
-	char *ID;
-	cudaMalloc((void **) &ID, N);						// Allocate device memory
-	// Copy to the device
-	cudaMemcpy(ID, id, N, cudaMemcpyHostToDevice);
-	//...........................................................................
-	
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	//...........................................................................
-	cudaMalloc((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
-	cudaMalloc((void **) &f_odd, 9*dist_mem_size);		// Allocate device memory
-	//...........................................................................
-
-	//...........................................................................
-//	cudaHostAlloc(&fa,dist_mem_size,cudaHostAllocPortable);
-//	cudaHostAlloc(&fb,dist_mem_size,cudaHostAllocPortable);
-//	cudaHostRegister(fa,dist_mem_size,cudaHostRegisterPortable);
-//	cudaHostRegister(fb,dist_mem_size,cudaHostRegisterPortable);
-//	cudaHostRegister(id,N*sizeof(char),cudaHostAllocPortable);
-
-	printf("Setting the distributions, size = : %i\n", N);
-	//...........................................................................
-//	INITIALIZE <<< grid, nthreads >>>  (ID, f_even, f_odd, Nx, Ny, Nz, S);
-	//...........................................................................
-	dvc_InitD3Q19(ID,f_even,f_odd,Nx,Ny,Nz,nBlocks,nthreads,S);
-	//*************************************************************************
-
-	int timestep = 0;
-	printf("No. of timesteps: %i \n", timestepMax);
-	
-	//.......create a stream for the LB calculation.......
-	cudaStream_t stream;
-	cudaStreamCreate(&stream);
-	
-	//.......create and start timer............
-	cudaEvent_t start, stop; 
-	float time; 
- 
-	cudaEventCreate(&start); 
-	cudaEventCreate(&stop); 
-	cudaEventRecord( start, 0 ); 
-	//.........................................
-	
-	//************ MAIN ITERATION LOOP ***************************************/
-	while (timestep < timestepMax){
-	
-		//...................................................................
-		//........ Execute the swap kernel (device) .........................
-//		SWAP <<< grid, nthreads >>> (ID, f_even, f_odd, Nx, Ny, Nz, S);
-		//...................................................................
-		dvc_SwapD3Q19(ID,f_even,f_odd,Nx,Ny,Nz,nBlocks,nthreads,S);
-
-		//........ Execute the collision kernel (device) ....................
-//		MRT <<< grid, nthreads >>> (ID, f_even, f_odd, Nx, Ny, Nz, S,
-//									rlx_setA, rlx_setB, Fx, Fy, Fz);
-		//............................................................
-		dvc_MRT(ID, f_even, f_odd, rlx_setA, rlx_setB, Fx, Fy, Fz,Nx,Ny,Nz,nBlocks,nthreads,S);
-		// Iteration completed!
-
-		timestep++;
-		//...................................................................
-		
-	}
-	//************************************************************************/
-	
-	cudaThreadSynchronize();
-	//.......... stop and destroy timer.............................
-	cudaEventRecord( stop, stream); 
-	cudaEventSynchronize( stop ); 
- 
-	cudaEventElapsedTime( &time, start, stop ); 
-	printf("CPU time = %f \n", time);
-	
-	float MLUPS = 0.001*float(Nx*Ny*Nz)*timestep/time;
-	printf("MLUPS = %f \n", MLUPS);
-
-	cudaStreamDestroy(stream);
-	cudaEventDestroy( start ); 
-	cudaEventDestroy( stop ); 
-	//..............................................................
-	
-	//..............................................................
-	//.........Compute the velocity and copy result to host ........
-	double *velocity;
-	velocity = new double[3*N];
-	//......................device distributions....................................
-	double *vel;
-	//..............................................................................
-	cudaMalloc((void **) &vel, 3*dist_mem_size);	// Allocate device memory
-	//..............................................................................
-//	Compute_VELOCITY <<< grid, nthreads >>>  (ID, f_even, f_odd, vel, Nx, Ny, Nz, S);
-	//..............................................................................
-	cudaMemcpy(velocity, vel, 3*dist_mem_size, cudaMemcpyDeviceToHost);
-	//..............................................................................
-
-	//............................................................	
-	//....Write the z-velocity to test poiseuille flow............
-	double vz,vz_avg;	
-	vz_avg = 0.0;
-
-	FILE *output;
-	output = fopen("velocity.out","w");
-	for (int k=0; k<1; k++){
-		for (int j=0; j<1; j++){
-			for (int i=0; i<Nx; i++){
-				int n = k*Nx*Ny+j*Nx+i;
-				//.....print value........
-				vz = velocity[2*N+n];
-				vz_avg += vz;
-				fprintf(output, " %e",vz);
-			}
-		}
-	}
-	fclose(output);
-	
-	vz = vz_avg/double(sum);
-	printf("Average Velocity = %e\n", vz);
-
-
-	// cleanup	
-	cudaFree(f_even);	cudaFree(f_odd);	cudaFree(vel);	cudaFree(ID);
-	free (velocity);	free(id);
-	
-}

cuda/exe/lb1_MRT.cu

@@ -1,246 +0,0 @@
-#include <stdio.h>
-#include <iostream>
-#include <fstream>
-#include <cuda.h>
-
-using namespace std;
-
-//*************************************************************************
-extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size);
-//*************************************************************************
-extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size);
-//*************************************************************************
-extern "C" void dvc_Barrier();
-//*************************************************************************
-extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx,
-							  int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-extern "C" void dvc_SwapD3Q19(char *ID, double *f_even, double *f_odd, int Nx,
-							  int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-extern "C" void dvc_MRT(char *ID, double *f_even, double *f_odd, double rlxA, double rlxB, double Fx, double Fy, double Fz,
-		int Nx, int Ny, int Nz, int nblocks, int nthreads, int S);
-//*************************************************************************
-
-void Write_Out(double *array, int Nx, int Ny, int Nz){
-	int value;
-	FILE *output;
-	output = fopen("dist.list","w");
-	for (int k=0; k<Nz; k++){
-		for (int j=0; j<Ny; j++){
-			for (int i=0; i<Nx; i++){
-				int index = k*Nx*Ny+j*Nx+i;
-				value = int(array[index]);
-				fprintf(output, "| %i",value);
-			}
-			fprintf(output, " | \n");
-		}
-		fprintf(output,"************************************** \n");	
-	}
-	fclose(output);
-}
-
-//**************************************************************************
-// MRT implementation of the LBM using CUDA
-//**************************************************************************
-
-int main(void)
-{
-
-	// BGK Model parameters
-	string FILENAME;	
-	unsigned int nBlocks, nthreads;
-	int timestepMax, interval;
-	double tau,Fx,Fy,Fz,tol;
-	// Domain variables
-	int Nx,Ny,Nz;
-
-	ifstream input("MRT.in");
-	input >> FILENAME;		// name of the input file
-	input >> Nz;			// number of nodes (x,y,z)
-	input >> nBlocks;				
-	input >> nthreads;				
-	input >> tau;				// relaxation time 
-	input >> Fx;			// External force components (x,y,z)
-	input >> Fy;
-	input >> Fz;
-	input >> timestepMax;			// max no. of timesteps
-	input >> interval;			// error interval
-	input >> tol;				// error tolerance
-	
-	double rlx_setA = 1.f/tau;
-	double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
-	
-	printf("tau = %f \n", tau);
-	printf("Set A = %f \n", rlx_setA);
-	printf("Set B = %f \n", rlx_setB);
-	printf("Force(x) = %f \n", Fx);
-	printf("Force(y) = %f \n", Fy);
-	printf("Force(z) = %f \n", Fz);
-
-	Nx = Ny = Nz;	// Cubic domain
-	
-	int N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
-	
-//	unsigned int nBlocks = 32;
-//	int nthreads = 128;
-	int S = N/nthreads/nBlocks;
-	
-//	unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1);
-	dim3 grid(nBlocks,1,1);
-		
-	printf("Number of blocks = %i \n", nBlocks);
-	printf("Threads per block = %i \n", nthreads);
-	printf("Sweeps per thread = %i \n", S);
-	printf("Number of nodes per side = %i \n", Nx);
-	printf("Total Number of nodes = %i \n", N);
-	
-	//.......................................................................
-	printf("Read input media... \n");
-	// .......... READ THE INPUT FILE .......................................
-	int n;
-	char value;
-	char *id;
-	id = new char[N];	
-	int sum = 0;
-	double porosity;
-	ifstream PM(FILENAME.c_str(),ios::binary);
-	for (int k=0;k<Nz;k++){
-		for (int j=0;j<Ny;j++){
-			for (int i=0;i<Nx;i++){
-				PM.read((char *) (&value), sizeof(value));
-				n = k*Nx*Ny+j*Nx+i;
-				id[n] = value;
-				if (value > 0) sum++;
-			}
-		}
-	}
-	PM.close();
-	printf("File porosity = %f\n", double(sum)/N);
-	//.......................................................................
-	//...........device phase ID.................................................
-	char *ID;
-	 dvc_AllocateDeviceMemory((void **) &ID, N);						// Allocate device memory
-	// Copy to the device
-	dvc_CopyToDevice(ID, id, N);
-	//...........................................................................
-	
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	//...........................................................................
-	 dvc_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
-	 dvc_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size);		// Allocate device memory
-	//...........................................................................
-
-	//...........................................................................
-//	cudaHostAlloc(&fa,dist_mem_size,cudaHostAllocPortable);
-//	cudaHostAlloc(&fb,dist_mem_size,cudaHostAllocPortable);
-//	cudaHostRegister(fa,dist_mem_size,cudaHostRegisterPortable);
-//	cudaHostRegister(fb,dist_mem_size,cudaHostRegisterPortable);
-//	cudaHostRegister(id,N*sizeof(char),cudaHostAllocPortable);
-
-	printf("Setting the distributions, size = : %i\n", N);
-	//...........................................................................
-//	INITIALIZE <<< grid, nthreads >>>  (ID, f_even, f_odd, Nx, Ny, Nz, S);
-	//...........................................................................
-	dvc_InitD3Q19(ID,f_even,f_odd,Nx,Ny,Nz,nBlocks,nthreads,S);
-	//*************************************************************************
-
-	int timestep = 0;
-	printf("No. of timesteps: %i \n", timestepMax);
-	
-	//.......create a stream for the LB calculation.......
-	cudaStream_t stream;
-	cudaStreamCreate(&stream);
-	
-	//.......create and start timer............
-	cudaEvent_t start, stop; 
-	float time; 
- 
-	cudaEventCreate(&start); 
-	cudaEventCreate(&stop); 
-	cudaEventRecord( start, 0 ); 
-	//.........................................
-	
-	//************ MAIN ITERATION LOOP ***************************************/
-	while (timestep < timestepMax){
-	
-		//...................................................................
-		//........ Execute the swap kernel (device) .........................
-//		SWAP <<< grid, nthreads >>> (ID, f_even, f_odd, Nx, Ny, Nz, S);
-		//...................................................................
-		dvc_SwapD3Q19(ID,f_even,f_odd,Nx,Ny,Nz,nBlocks,nthreads,S);
-
-		//........ Execute the collision kernel (device) ....................
-//		MRT <<< grid, nthreads >>> (ID, f_even, f_odd, Nx, Ny, Nz, S,
-//									rlx_setA, rlx_setB, Fx, Fy, Fz);
-		//............................................................
-		dvc_MRT(ID, f_even, f_odd, rlx_setA, rlx_setB, Fx, Fy, Fz,Nx,Ny,Nz,nBlocks,nthreads,S);
-		// Iteration completed!
-
-		timestep++;
-		//...................................................................
-		
-	}
-	//************************************************************************/
-	
-//	cudaThreadSynchronize();
-	dvc_Barrier();
-	//.......... stop and destroy timer.............................
-	cudaEventRecord( stop, stream); 
-	cudaEventSynchronize( stop ); 
- 
-	cudaEventElapsedTime( &time, start, stop ); 
-	printf("CPU time = %f \n", time);
-	
-	float MLUPS = 0.001*float(Nx*Ny*Nz)*timestep/time;
-	printf("MLUPS = %f \n", MLUPS);
-
-	cudaStreamDestroy(stream);
-	cudaEventDestroy( start ); 
-	cudaEventDestroy( stop ); 
-	//..............................................................
-	
-	//..............................................................
-	/*//.........Compute the velocity and copy result to host ........
-	double *velocity;
-	velocity = new double[3*N];
-	//......................device distributions....................................
-	double *vel;
-	//..............................................................................
-	 dvc_AllocateDeviceMemory((void **) &vel, 3*dist_mem_size);	// Allocate device memory
-	//..............................................................................
-//	Compute_VELOCITY <<< grid, nthreads >>>  (ID, f_even, f_odd, vel, Nx, Ny, Nz, S);
-	//..............................................................................
-//	cudaMemcpy(velocity, vel, 3*dist_mem_size, cudaMemcpyDeviceToHost);
-	//..............................................................................
-
-	//............................................................	
-	//....Write the z-velocity to test poiseuille flow............
-	double vz,vz_avg;	
-	vz_avg = 0.0;
-
-/*	FILE *output;
-	output = fopen("velocity.out","w");
-	for (int k=0; k<1; k++){
-		for (int j=0; j<1; j++){
-			for (int i=0; i<Nx; i++){
-				int n = k*Nx*Ny+j*Nx+i;
-				//.....print value........
-				vz = velocity[2*N+n];
-				vz_avg += vz;
-				fprintf(output, " %e",vz);
-			}
-		}
-	}
-	fclose(output);
-	
-	vz = vz_avg/double(sum);
-	printf("Average Velocity = %e\n", vz);
-*/
-	// cleanup	
-//	cudaFree(f_even);	cudaFree(f_odd);	cudaFree(vel);	cudaFree(ID);
-//	free (velocity);	free(id);
-	
-}

cuda/exe/lb2_Color.cpp

@@ -1,372 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <iostream>
-#include <fstream>
-
-//*************************************************************************
-// Functions defined in Color.cu
-//*************************************************************************
-extern "C" void dvc_InitDenColor( int nblocks, int nthreads, int S,
-		char *ID, double *Den, double *Phi, double das, double dbs, int N);
-//*************************************************************************
-extern "C" void dvc_ComputeColorGradient(int nBlocks, int nthreads, int S,
-		char *ID, double *Phi, double *ColorGrad, int Nx, int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_ColorCollide(int nBlocks, int nthreads, int S,
-		char *ID, double *f_even, double *f_odd, double *ColorGrad, double *Velocity,
-		double rlxA, double rlxB,double alpha, double beta, double Fx, double Fy, double Fz,
-		int Nx, int Ny, int Nz, bool pBC);
-//*************************************************************************
-extern "C" void dvc_DensityStreamD3Q7(int nBlocks, int nthreads, int S,
-		char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity,
-		double beta, int Nx, int Ny, int Nz, bool pBC);
-//*************************************************************************
-extern "C" void dvc_ComputePhi(int nBlocks, int nthreads, int S,
-		char *ID, double *Phi, double *Copy, double *Den, int N);
-//*************************************************************************
-//*************************************************************************
-// Functions defined in D3Q19.cu
-//*************************************************************************
-extern "C" void dvc_InitD3Q19(int nblocks, int nthreads, int S, char *ID, double *f_even, double *f_odd, int Nx,
-							  int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_SwapD3Q19(int nblocks, int nthreads, int S,
-		char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_PackDist(int grid, int threads, int q, int *SendList, int start,
-		int sendCount, double *sendbuf, double *Dist, int N);
-//*************************************************************************
-extern "C" void dvc_UnpackDist(int grid, int threads, int q, int Cqx, int Cqy, int Cqz, int *RecvList, int start,
-		int recvCount, double *recvbuf, double *Dist, int Nx, int Ny, int Nz);
-//*************************************************************************
-//***************************************************************************************
-// Functions defined in D3Q7.cu
-//***************************************************************************************
-extern "C" void dvc_PackDenD3Q7(int grid, int threads, int *list, int count, double *sendbuf,
-		int number, double *Data, int N);
-//***************************************************************************************
-extern "C" void dvc_UnpackDenD3Q7(int grid, int threads, int *list, int count, double *recvbuf,
-		int number, double *Data, int N);
-//***************************************************************************************
-extern "C" void dvc_PackValues(int grid, int threads, int *list, int count, double *sendbuf,
-		double *Data, int N);
-//***************************************************************************************
-extern "C" void dvc_UnpackValues(int grid, int threads, int *list, int count, double *recvbuf,
-		double *Data, int N);
-//***************************************************************************************
-//*************************************************************************
-// Functions defined in CudaExtras.cu
-//*************************************************************************
-extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size);
-//*************************************************************************
-extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size);
-//*************************************************************************
-extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size);
-//*************************************************************************
-extern "C" void dvc_Barrier();
-//*************************************************************************
-
-//*************************************************************************
-// Implementation of Two-Phase Immiscible LBM using CUDA
-//*************************************************************************
-
-using namespace std;
-
-
-inline void PackID(int *list, int count, char *sendbuf, char *ID){
-	// Fill in the phase ID values from neighboring processors
-	// This packs up the values that need to be sent from one processor to another
-	int idx,n;
-
-	for (idx=0; idx<count; idx++){
-		n = list[idx];
-		sendbuf[idx] = ID[n];
-	}
-}
-//***************************************************************************************
-inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
-	// Fill in the phase ID values from neighboring processors
-	// This unpacks the values once they have been recieved from neighbors
-	int idx,n;
-
-	for (idx=0; idx<count; idx++){
-		n = list[idx];
-		ID[n] = recvbuf[idx];
-	}
-}
-//***************************************************************************************
-
-int main(int argc, char **argv)
-{
-
-	int rank = 0;
-	int nprocs =1;
-	int nprocx,nprocy,nprocz;
-	int iproc,jproc,kproc;
-	if (rank == 0){
-		printf("********************************************************\n");
-		printf("Running Hybrid Implementation of Color LBM	\n");
-		printf("********************************************************\n");
-	}
-	// Color Model parameters
-	string FILENAME;
-	unsigned int nBlocks, nthreads;
-	int Nx,Ny,Nz;
-	int timestepMax, interval;
-	double tau,Fx,Fy,Fz,tol;
-	double alpha, beta;
-	double das, dbs;
-	double din,dout;
-	bool pBC;
-	int i,j,k,n;
-
-	if (rank==0){
-		//.............................................................
-		//		READ SIMULATION PARMAETERS FROM INPUT FILE
-		//.............................................................
-		ifstream input("Color.in");
-		// Line 1: Name of the phase indicator file (s=0,w=1,n=2)
-		input >> FILENAME;
-		// Line 2: domain size (Nx, Ny, Nz)
-		input >> Nz;				// number of nodes (x,y,z)
-		input >> nBlocks;
-		input >> nthreads;
-		// Line 3: model parameters (tau, alpha, beta, das, dbs)
-		input >> tau;
-		input >> alpha;
-		input >> beta;
-		input >> das;
-		input >> dbs;
-		// Line 4: External force components (Fx,Fy, Fz)
-		input >> Fx;
-		input >> Fy;
-		input >> Fz;
-		// Line 5: Pressure Boundary conditions
-		input >> pBC;
-		input >> din;
-		input >> dout;
-		// Line 6: time-stepping criteria
-		input >> timestepMax;			// max no. of timesteps
-		input >> interval;			// error interval
-		input >> tol;				// error tolerance
-		//.............................................................
-
-		ifstream domain("Domain.in");
-		domain >> nprocx;
-		domain >> nprocy;
-		domain >> nprocz;
-	}
-
-	double rlxA = 1.f/tau;
-	double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
-
-	if (nprocs != nprocx*nprocy*nprocz){
-		printf("Fatal error in processor number! \n");
-		printf("nprocx =  %i \n",nprocx);
-		printf("nprocy =  %i \n",nprocy);
-		printf("nprocz =  %i \n",nprocz);
-	}
-
-	if (rank==0){
-		printf("********************************************************\n");
-		printf("tau = %f \n", tau);
-		printf("alpha = %f \n", alpha);
-		printf("beta = %f \n", beta);
-		printf("das = %f \n", beta);
-		printf("dbs = %f \n", beta);
-		printf("Force(x) = %f \n", Fx);
-		printf("Force(y) = %f \n", Fy);
-		printf("Force(z) = %f \n", Fz);
-		printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
-		printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
-		printf("********************************************************\n");
-
-	}
-
-	Nz += 2;
-	Nx = Ny = Nz;	// Cubic domain
-
-	int N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
-
-//	unsigned int nBlocks = 32;
-//	int nthreads = 128;
-	int S = N/nthreads/nBlocks;
-
-//	unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1);
-//	dim3 grid(nBlocks,1,1);
-
-	if (rank==0) printf("Number of blocks = %i \n", nBlocks);
-	if (rank==0) printf("Threads per block = %i \n", nthreads);
-	if (rank==0) printf("Sweeps per thread = %i \n", S);
-	if (rank==0) printf("Number of nodes per side = %i \n", Nx);
-	if (rank==0) printf("Total Number of nodes = %i \n", N);
-	if (rank==0) printf("********************************************************\n");
-
-	//.......................................................................
-	if (rank == 0)	printf("Read input media... \n");
-	//.......................................................................
-	char LocalRankString[8];
-	char LocalRankFilename[40];
-	sprintf(LocalRankString,"%05d",rank);
-	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-//	printf("Local File Name =  %s \n",LocalRankFilename);
-	// .......... READ THE INPUT FILE .......................................
-	char value;
-	char *id;
-	id = new char[N];
-	int sum = 0;
-//	double porosity;
-	//.......................................................................
-	ifstream PM(LocalRankFilename,ios::binary);
-	for (k=0;k<Nz;k++){
-		for (j=0;j<Ny;j++){
-			for (i=0;i<Nx;i++){
-				n = k*Nx*Ny+j*Nx+i;
-				id[n] = 0;
-			}
-		}
-	}
-	for ( k=1;k<Nz-1;k++){
-		for ( j=1;j<Ny-1;j++){
-			for ( i=1;i<Nx-1;i++){
-				PM.read((char *) (&value), sizeof(value));
-				n = k*Nx*Ny+j*Nx+i;
-				id[n] = value;
-				if (value > 0) sum++;
-			}
-		}
-	}
-	PM.close();
-//	printf("File porosity = %f\n", double(sum)/N);
-
-	//...........device phase ID.................................................
-	if (rank==0)	printf ("Copying phase ID to device \n");
-	char *ID;
-	dvc_AllocateDeviceMemory((void **) &ID, N);						// Allocate device memory
-	// Copy to the device
-	dvc_CopyToDevice(ID, id, N);
-	//...........................................................................
-
-	if (rank==0)	printf ("Allocating distributions \n");
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	//...........................................................................
-	dvc_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
-	dvc_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size);	// Allocate device memory
-	//...........................................................................
-	//...........................................................................
-	//				MAIN  VARIABLES ALLOCATED HERE
-	//...........................................................................
-	double *Phi,*Den,*Copy;
-	double *ColorGrad, *Velocity;
-	//...........................................................................
-	dvc_AllocateDeviceMemory((void **) &Phi, dist_mem_size);
-	dvc_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
-	dvc_AllocateDeviceMemory((void **) &Copy, 2*dist_mem_size);
-	dvc_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
-	dvc_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
-	//...........................................................................
-	if (rank==0)	printf("Setting the distributions, size = : %i\n", N);
-	//...........................................................................
-	dvc_InitD3Q19(nBlocks, nthreads, S, ID, f_even, f_odd, Nx, Ny, Nz);
-	dvc_InitDenColor(nBlocks, nthreads, S, ID, Den, Phi,  das, dbs, N);
-	//...........................................................................
-	dvc_ComputePhi(nBlocks, nthreads, S,ID, Phi, Copy, Den, N);
-	//...........................................................................
-	
-	//...........................................................................
-	// Grids used to pack faces on the GPU for MPI
-	int faceGrid,edgeGrid,packThreads;
-	packThreads=512;
-	edgeGrid=1;
-	faceGrid=Nx*Ny/packThreads;
-
-
-	int timestep = 0;
-	if (rank==0) printf("********************************************************\n");
-	if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
-
-	//.......create a stream for the LB calculation.......
-//	cudaStream_t stream;
-//	cudaStreamCreate(&stream);
-
-	//.......create and start timer............
-	double start,stop;
-	double walltime;
-	start = clock();
-
-
-	//************ MAIN ITERATION LOOP ***************************************/
-	while (timestep < timestepMax){
-
-		//*************************************************************************
-		// 		Compute the color gradient
-		//*************************************************************************
-		dvc_ComputeColorGradient(nBlocks, nthreads, S,
-				ID, Phi, ColorGrad, Nx, Ny, Nz);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Perform collision step for the momentum transport
-		//*************************************************************************
-		dvc_ColorCollide(nBlocks, nthreads, S,
-				ID, f_even, f_odd, ColorGrad, Velocity,
-				rlxA, rlxB,alpha, beta, Fx, Fy, Fz, Nx, Ny, Nz, pBC);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Carry out the density streaming step for mass transport
-		//*************************************************************************
-		dvc_DensityStreamD3Q7(nBlocks, nthreads, S,
-				ID, Den, Copy, Phi, ColorGrad, Velocity,beta, Nx, Ny, Nz, pBC);
-		//*************************************************************************
-		
-		//*************************************************************************
-		// 		Swap the distributions for momentum transport
-		//*************************************************************************
-		dvc_SwapD3Q19(nBlocks, nthreads, S, ID, f_even, f_odd, Nx, Ny, Nz);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Compute the phase indicator field and reset Copy, Den
-		//*************************************************************************
-		dvc_ComputePhi(nBlocks, nthreads, S,ID, Phi, Copy, Den, N);
-		//*************************************************************************
-
-		// Iteration completed!
-		timestep++;
-		
-		//...................................................................
-	}
-	//************************************************************************/
-	dvc_Barrier();
-	stop = clock();
-
-//	cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl;
-	walltime = (stop - start)/CLOCKS_PER_SEC;
-//	cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 <<  " MLUPS" << endl;
-	double MLUPS = double(Nx*Ny*Nz*timestep)/walltime/1000000;
-	if (rank==0) printf("********************************************************\n");
-	if (rank==0) printf("CPU time = %f \n", walltime);
-	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");
-	
-	//************************************************************************/
-	// Write out the phase indicator field 
-	//************************************************************************/
-	sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
-	//	printf("Local File Name =  %s \n",LocalRankFilename);
-	double *phiOut;
-	phiOut = new double[N];
-	dvc_CopyToHost(phiOut,Phi,N*sizeof(double));
-		
-	FILE *PHASE;
-	PHASE = fopen(LocalRankFilename,"wb");
-	fwrite(phiOut,8,N,PHASE);
-	fclose(PHASE);
-	//************************************************************************/
-
-}

cuda/exe/lb2_Color.cu

@@ -1,425 +0,0 @@
-#ifdef useMPI
-#include <mpi.h>
-#endif
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <iostream>
-#include <fstream>
-#include <math.h>
-#include <cuda.h>
-
-using namespace std;
-//*************************************************************************
-// HokieSpeed
-//nvcc -Xcompiler -fopenmp -lgomp -O3 -arch sm_20 -o hybridATLKR lb2_ATLKR_hybrid.cu
-// -I$VT_MPI_INC -L$VT_MPI_LIB -lmpi
-//*************************************************************************
-
-//*************************************************************************
-// Implementation of Two-Phase Immiscible LBM using CUDA
-//*************************************************************************
-
-//*************************************************************************
-extern "C" void dvc_InitD3Q19(int nblocks, int nthreads, int S,
-		char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_InitDenColor( int nblocks, int nthreads, int S,
-		char *ID, double *Den, double *Phi, double das, double dbs, int N);
-//*************************************************************************
-extern "C" void dvc_ComputeColorGradient(int nBlocks, int nthreads, int S,
-		char *ID, double *Phi, double *ColorGrad, int Nx, int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_ColorCollide(int nBlocks, int nthreads, int S,
-		char *ID, double *f_even, double *f_odd, double *ColorGrad, double *Velocity,
-		double rlxA, double rlxB,double alpha, double beta, double Fx, double Fy, double Fz,
-		int Nx, int Ny, int Nz, bool pBC);
-//*************************************************************************
-extern "C" void dvc_DensityStreamD3Q7(int nBlocks, int nthreads, int S,
-		char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity,
-		double beta, int Nx, int Ny, int Nz, bool pBC);
-//*************************************************************************
-extern "C" void dvc_ComputePhi(int nBlocks, int nthreads, int S,
-		char *ID, double *Phi, double *Copy, double *Den, int N);
-//*************************************************************************
-extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size);
-//*************************************************************************
-extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size);
-//*************************************************************************
-extern "C" void dvc_Barrier();
-//*************************************************************************
-extern "C" void dvc_SwapD3Q19(int nblocks, int nthreads, int S,
-		char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);
-//*************************************************************************
-extern "C" void dvc_PackDist(int grid, int threads, int q, int *SendList, int start,
-		int sendCount, double *sendbuf, double *Dist, int N);
-//*************************************************************************
-extern "C" void dvc_UnpackDist(int grid, int threads, int q, int Cqx, int Cqy, int Cqz, int *RecvList, int start,
-		int recvCount, double *recvbuf, double *Dist, int Nx, int Ny, int Nz);
-//*************************************************************************
-
-int main(int argc, char *argv[])
-{
-	
-	//********** Initialize MPI ****************
-	int numprocs,rank;
-#ifdef useMPI
-	MPI_Status stat;
-	MPI_Init(&argc,&argv);
-    MPI_Comm comm = MPI_COMM_WORLD;
-	MPI_Comm_size(comm,&numprocs);
-	MPI_Comm_rank(comm,&rank);
-#else
-    MPI_Comm comm = MPI_COMM_WORLD;
-	numprocs = 1;
-	rank = 0;
-#endif
-	//******************************************
-	
-	if (rank == 0){
-		printf("********************************************************\n");
-		printf("Running Hybrid Implementation of Color LBM	\n");
-		printf("********************************************************\n");
-	}
-	// Color Model parameters
-	string FILENAME;
-	unsigned int nBlocks, nthreads;
-	int Nx,Ny,Nz;
-	int timestepMax, interval;
-	double tau,Fx,Fy,Fz,tol;
-	double alpha, beta;
-	double das, dbs;
-	double din,dout;
-	bool pBC;
-
-	if (rank==0){
-		//.............................................................
-		//		READ SIMULATION PARMAETERS FROM INPUT FILE 
-		//.............................................................
-		ifstream input("Color.in");
-		// Line 1: Name of the phase indicator file (s=0,w=1,n=2)
-		input >> FILENAME;				
-		// Line 2: domain size (Nx, Ny, Nz)
-		input >> Nz;				// number of nodes (x,y,z)
-		input >> nBlocks;				
-		input >> nthreads;			
-		// Line 3: model parameters (tau, alpha, beta, das, dbs)
-		input >> tau;
-		input >> alpha;
-		input >> beta;
-		input >> das;	
-		input >> dbs;
-		// Line 4: External force components (Fx,Fy, Fz)
-		input >> Fx;				
-		input >> Fy;
-		input >> Fz;
-		// Line 5: Pressure Boundary conditions
-		input >> pBC;
-		input >> din;
-		input >> dout;
-		// Line 6: time-stepping criteria
-		input >> timestepMax;			// max no. of timesteps
-		input >> interval;			// error interval
-		input >> tol;				// error tolerance
-		//.............................................................
-	}
-#ifdef useMPI
-	// **************************************************************
-	// Broadcast simulation parameters from rank 0 to all other procs
-	MPI_Barrier(comm);
-	//.................................................
-	MPI_Bcast(&Nz,1,MPI_INT,0,comm);
-	MPI_Bcast(&nBlocks,1,MPI_INT,0,comm);
-	MPI_Bcast(&nthreads,1,MPI_INT,0,comm);
-	MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm);
-	MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
-	MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
-	MPI_Bcast(&interval,1,MPI_INT,0,comm);
-	MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm);
-	//.................................................
-	MPI_Barrier(comm);
-	// **************************************************************
-#endif
-	
-	double rlxA = 1.f/tau;
-	double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
-
-	if (pBC && rank == 0){
-		printf("Assigning presusre boundary conditions \n");
-		printf("Inlet density = %f \n", din);
-		printf("Outlet density = %f \n", dout);
-	}
-
-	if (rank==0){
-		printf("....Parameters................\n");
-		printf("tau = %f \n", tau);
-		printf("alpha = %f \n", alpha);
-		printf("beta = %f \n", beta);
-		printf("das = %f \n", das);
-		printf("dbs = %f \n", dbs);
-		printf("Force(x) = %f \n", Fx);
-		printf("Force(y) = %f \n", Fy);
-		printf("Force(z) = %f \n", Fz);
-		printf("Nz = %i \n", Nz);
-		printf("timestepMax = %i \n", timestepMax);
-		printf("...............................\n");
-	}
-	
-	// Identical cubic sub-domains
-	Nx = Ny = Nz;// = 16*s;		// Cubic domain
-	int N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
-
-//	unsigned int nBlocks = 32;
-//	int nthreads = 128;
-	int S = N/nthreads/nBlocks;
-	if (nBlocks*nthreads*S < N)	S++;
-//	int S = 1;
-	
-//	unsigned int nBlocks = N/nthreads + (N%nthreads == 0?0:1);
-//	dim3 grid(nBlocks,1,1);
-	if (rank==1){
-		printf("Number of blocks = %i \n", nBlocks);
-		printf("Threads per block = %i \n", nthreads);
-		printf("Sweeps per thread = %i \n", S);
-		printf("Number of nodes per side = %i \n", Nx);
-		printf("Total Number of nodes = %i \n", N);
-		printf("...............................\n");
-	}
-	
-	//.......................................................................
-	// .......... READ THE INPUT FILE .......................................
-	int n;
-	char value;
-	char *id;
-	id = new char[N];	
-	int sum = 0;
-	// RANK 0 READS THE INPUT FILE	
-	if (rank==0){
-		printf("Read input media... \n");
-		ifstream PM(FILENAME.c_str(),ios::binary);
-		for (int k=0;k<Nz;k++){
-			for (int j=0;j<Ny;j++){
-				for (int i=0;i<Nx;i++){
-					PM.read((char *) (&value), sizeof(value));
-					n = k*Nx*Ny+j*Nx+i;
-
-					if (value>0){
-						if (pBC) value=2; 	// Saturate with NWP
-						if (k<8){
-							value=1;
-						}
-					}
-
-					id[n] = value;
-					if (value > 0) sum++;
-				}
-			}
-		}
-		PM.close();
-		printf("File porosity = %f\n", double(sum)/N);
-	}
-	//......... for pressure BC only............................
-	// Void the first / last rows if pressure BC are to be used
-	if (pBC){
-		for (int k=0;k<Nz;k++){
-			for (int j=0;j<Ny;j++){
-				for (int i=0;i<Nx;i++){
-					n = k*Nx*Ny+j*Nx+i;
-					if (k<4) id[n] = 1;
-					if (k>Nz-5) id[n] = 2;
-				}
-			}
-			// Skip the non-boundary values
-			if (k==4)	k=Nz-5;
-		}
-	}
-#ifdef useMPI	//............................................................
-	MPI_Barrier(comm);
-	MPI_Bcast(&id[0],N,MPI_CHAR,0,comm);
-	MPI_Barrier(comm);
-#endif
-	if (rank == 0) printf("Domain set.\n");
-	//...........................................................................
-
-	int SBC;
-	int outlet = N-Nx*Ny;
-	if (pBC){
-		SBC = Nx*Ny/nthreads/nBlocks+1;
-		printf("Number of sweeps for inlet / outlet: %i \n", SBC);
-	}
-	//...........................................................................
-	
-	//...........................................................................
-	//...........device phase ID.................................................
-	char *ID;
-	cudaMalloc((void **) &ID, N);						// Allocate device memory
-	// Copy to the device
-	cudaMemcpy(ID, id, N, cudaMemcpyHostToDevice);
-	//...........................................................................
-	
-	//......................device distributions.................................
-	double *f_even,*f_odd;
-	//...........................................................................
-	cudaMalloc((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
-	cudaMalloc((void **) &f_odd, 9*dist_mem_size);		// Allocate device memory
-//	f_even = new double[10*N];
-//	f_odd = new double[9*N];
-	//...........................................................................
-
-	//...........................................................................
-	//				MAIN  VARIABLES ALLOCATED HERE
-	//...........................................................................
-	double *Phi,*Den,*Copy;
-	double *ColorGrad, *Velocity;
-	//...........................................................................
-	cudaMalloc((void **) &Phi, dist_mem_size);
-	cudaMalloc((void **) &Den, 2*dist_mem_size);
-	cudaMalloc((void **) &Copy, 2*dist_mem_size);
-	cudaMalloc((void **) &Velocity, 3*dist_mem_size);
-	cudaMalloc((void **) &ColorGrad, 3*dist_mem_size);
-	//...........................................................................
-	
-	//...........................................................................
-	if (rank==0)	printf("Setting the distributions, size = : %i\n", N);
-	//...........................................................................
-	dvc_InitD3Q19(nBlocks, nthreads, S, ID, f_even, f_odd, Nx, Ny, Nz);
-	dvc_InitDenColor(nBlocks, nthreads, S, ID, Den, Phi,  das, dbs, N);
-	//...........................................................................
-	dvc_ComputePhi(nBlocks, nthreads, S,ID, Phi, Copy, Den, N);
-	//...........................................................................
-
-	int timestep;
-//	double starttime,stoptime;
-	if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
-	timestep = 0;
-	//.......create and start timer............
-	cudaEvent_t start, stop;
-	float time;
-	//.......create a stream for the LB calculation.......
-	cudaStream_t stream;
-	cudaStreamCreate(&stream);
-
-	cudaEventCreate(&start);
-	cudaEventCreate(&stop);
-	cudaEventRecord( start, 0 );
-	//.........................................
-	//************ MAIN TIMESTEP LOOP ***************************************/
-	while (timestep < timestepMax){
-
-		//*************************************************************************
-		// 		Compute the color gradient
-		//*************************************************************************
-		dvc_ComputeColorGradient(nBlocks, nthreads, S,
-				ID, Phi, ColorGrad, Nx, Ny, Nz);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Perform collision step for the momentum transport
-		//*************************************************************************
-		dvc_ColorCollide(nBlocks, nthreads, S,
-				ID, f_even, f_odd, ColorGrad, Velocity,
-				rlxA, rlxB,alpha, beta, Fx, Fy, Fz, Nx, Ny, Nz, pBC);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Carry out the density streaming step for mass transport
-		//*************************************************************************
-		dvc_DensityStreamD3Q7(nBlocks, nthreads, S,
-				ID, Den, Copy, Phi, ColorGrad, Velocity,beta, Nx, Ny, Nz, pBC);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Swap the distributions for momentum transport
-		//*************************************************************************
-		dvc_SwapD3Q19(nBlocks, nthreads, S, ID, f_even, f_odd, Nx, Ny, Nz);
-		//*************************************************************************
-
-		//*************************************************************************
-		// 		Compute the phase indicator field and reset Copy, Den
-		//*************************************************************************
-		dvc_ComputePhi(nBlocks, nthreads, S,ID, Phi, Copy, Den, N);
-		//*************************************************************************
-
-		dvc_Barrier();
-		timestep++;
-		//.............................................................................
-	}
-	//************************************************************************/
-	dvc_Barrier();
-	//.......... stop and destroy timer.............................
-	cudaEventRecord( stop, stream);
-	cudaEventSynchronize( stop );
-
-	cudaEventElapsedTime( &time, start, stop );
-	printf("CPU time = %f \n", time);
-
-	float MLUPS = 0.001*float(Nx*Ny*Nz)*timestep/time;
-	printf("MLUPS = %f \n", MLUPS);
-
-	cudaEventDestroy( start );
-	cudaEventDestroy( stop );
-
-	double *Data;
-	Data = new double[3*N];
-
-	cudaMemcpy(Data, Phi, dist_mem_size, cudaMemcpyDeviceToHost);
-
-	// Write out the Phase Indicator Field
-	FILE *phase;
-	phase = fopen("Phase.out","wb");
-	fwrite(Data,8,N,phase);
-	fclose(phase);
-
-	//....................................................
-	// Write out the pressure - (reuse Phi arrays since we're done with those)
-//	ComputeDensity<<< grid, nthreads>>> (ID, f_even, f_odd, Phi, Nx, Ny, Nz, S);
-//	cudaMemcpy(Data, Phi, dist_mem_size, cudaMemcpyDeviceToHost);
-//	FILE *PRESSURE;
-//	PRESSURE = fopen("Pressure.out","wb");
-//	fwrite(Phi,8,N,PRESSURE);
-//	fclose(PRESSURE);
-	//....................................................
-
-	// Write out the Color Gradient
-
-	cudaMemcpy(Data, ColorGrad, 3*dist_mem_size, cudaMemcpyDeviceToHost);
-
-	FILE *CG;
-	CG = fopen("ColorGrad.out","wb");
-	fwrite(Data,8,3*N,CG);
-	fclose(CG);
-	
-	// Write out the Velocity
-//	FILE *VEL;
-//	VEL = fopen("Velocity.out","wb");
-//	fwrite(Velocity,8,3*N,VEL);
-//	fclose(VEL);
-
-	// cleanup	
-	cudaFree(ID);
-	cudaFree(f_even);	cudaFree(f_odd);	
-	cudaFree(Velocity);
-	cudaFree(Phi);
-	
-	cudaFree (ColorGrad);
-	cudaFree (Den);		cudaFree(Copy);
-	cudaFree (Phi);
-	free(id);
-	
-	//***********Finish up!*********************************
-#ifdef useMPI
-	MPI_Finalize();
-#endif
-	return 0;
-	
-}