From a267c7b38b37071676dce680b705fec44894a96d Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Wed, 3 Sep 2025 14:51:23 -0300 Subject: [PATCH 01/29] permeability correction - MRT-routine --- models/MRTModel.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index ae3f530e..8fe19172 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -391,8 +391,7 @@ void ScaLBL_MRTModel::Run() { Xs = Dm->Comm.sumReduce(Xs); double h = Dm->voxel_length; - double absperm = - h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag; + double absperm = h * h * mu * Mask->Porosity() * flow_rate / force_mag; absperm *= 1013.0; // Convert to mDarcy if (rank == 0) { From 9b920db0abc458a3204fbbab6e20a5de8ce5ac48 Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Mon, 8 Sep 2025 15:13:28 -0300 Subject: [PATCH 02/29] Reader.cpp: including silo.h --- IO/Reader.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/IO/Reader.cpp b/IO/Reader.cpp index 566c07cb..0634545a 100644 --- a/IO/Reader.cpp +++ b/IO/Reader.cpp @@ -4,6 +4,7 @@ #include "IO/Mesh.h" #include "IO/MeshDatabase.h" #include "common/Utilities.h" +#include "IO/silo.h" #include #include From a47cfbf6f534b5b2ef4fd62a2c21d7d120ef5ea2 Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Tue, 9 Sep 2025 18:42:16 -0300 Subject: [PATCH 03/29] permeabilit.csv --- models/MRTModel.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 8fe19172..2baf2cb3 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -263,7 +263,7 @@ void ScaLBL_MRTModel::Run() { if (WriteHeader) { log_file = fopen("Permeability.csv", "a+"); - fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm\n"); + fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm(mDa) absperm*(mDa)\n"); fclose(log_file); } } @@ -399,9 +399,9 @@ void ScaLBL_MRTModel::Run() { FILE *log_file = fopen("Permeability.csv", "a"); fprintf(log_file, "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g " - "%.8g %.8g\n", + "%.8g %.8g %.8g\n", timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As, - h * Hs, Xs, vax, vay, vaz, absperm); + h * Hs, Xs, vax, vay, vaz, absperm, absperm * Mask->Porosity()); fclose(log_file); } } From 3e1427d10a53946847f28c5c1e902b42ef1f89ef Mon Sep 17 00:00:00 2001 From: Maria Rosa Amorim Date: Wed, 24 Sep 2025 18:36:01 +0200 Subject: [PATCH 04/29] tau correlation cpu --- cpu/Color.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 526bcafa..5e67478f 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1472,7 +1472,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = tauA + 0.5 * (1.0 - phi) * (tauB - tauA); + tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); @@ -2087,7 +2087,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = tauA + 0.5 * (1.0 - phi) * (tauB - tauA); + tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); From 55db011eae4fd28d5ef0624e4210e77a4d7c863a Mon Sep 17 00:00:00 2001 From: Maria Rosa Amorim Date: Wed, 24 Sep 2025 18:39:26 +0200 Subject: [PATCH 05/29] tau correlation gpu --- cuda/Color.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cuda/Color.cu b/cuda/Color.cu index c25e0f08..d32e0799 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -1311,7 +1311,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); + tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); @@ -1898,7 +1898,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); + tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); From 81c959365e5fead2f1fb4eb5fa12d4f640cac532 Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Mon, 6 Oct 2025 14:31:00 -0300 Subject: [PATCH 06/29] missing * --- cpu/Color.cpp | 4 ++-- cuda/Color.cu | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 5e67478f..a4c8cc07 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1472,7 +1472,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); @@ -2087,7 +2087,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); diff --git a/cuda/Color.cu b/cuda/Color.cu index d32e0799..a56d6f72 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -1311,7 +1311,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); @@ -1898,7 +1898,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); From 96524fa153e396566a439187b8aaeda47d268696 Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Mon, 6 Oct 2025 17:23:57 -0300 Subject: [PATCH 07/29] missing others * --- cpu/Color.cpp | 4 ++-- cuda/Color.cu | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/cpu/Color.cpp b/cpu/Color.cpp index a4c8cc07..30a7967d 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1472,7 +1472,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); @@ -2087,7 +2087,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); diff --git a/cuda/Color.cu b/cuda/Color.cu index a56d6f72..e5c92b65 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -1311,7 +1311,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); @@ -1898,7 +1898,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); From da9c52aca809860f395514f7d6b36c2713c19dec Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Tue, 7 Oct 2025 14:28:05 -0300 Subject: [PATCH 08/29] Update Color.cpp Change the constant in the interpolation of the viscosity to make it simpler --- cpu/Color.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 30a7967d..b4356126 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1472,7 +1472,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); From cd43ca77a896f574c99051af436b5ac449877e23 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Tue, 7 Oct 2025 14:29:35 -0300 Subject: [PATCH 09/29] Update Color.cpp Change interpolation function of the viscosity to make it simpler. --- cpu/Color.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpu/Color.cpp b/cpu/Color.cpp index b4356126..9e32147c 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -2087,7 +2087,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( // local density rho0 = rhoA + 0.5 * (1.0 - phi) * (rhoB - rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); From 617d50946e804f590b9ccf6e4616ba17a639fdc6 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Tue, 7 Oct 2025 14:30:46 -0300 Subject: [PATCH 10/29] Update Color.cu Change interpolation function of the viscosity to make it simpler. --- cuda/Color.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cuda/Color.cu b/cuda/Color.cu index e5c92b65..fc35c5c9 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -1311,7 +1311,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); @@ -1898,7 +1898,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // local density rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); // local relaxation time - tau = (tauA+tauB)*0.5 + tanh(phi*180.0/3.14159265359)*(tauA-tauB)*0.5; + tau = (tauA+tauB)*0.5 + tanh(phi*60.)*(tauA-tauB)*0.5; rlx_setA = 1.f/tau; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); From 4386eef5f44a1f941f71ce14eaba0f00b969b4ed Mon Sep 17 00:00:00 2001 From: Ricardoleite Date: Tue, 7 Oct 2025 15:43:11 -0300 Subject: [PATCH 11/29] Second-order velocity field Calculation --- common/ScaLBL.h | 2 ++ cpu/D3Q19.cpp | 45 ++++++++++++++++++++++++++++++++++ cuda/D3Q19.cu | 60 +++++++++++++++++++++++++++++++++++++++++++++ models/MRTModel.cpp | 2 +- 4 files changed, 108 insertions(+), 1 deletion(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 7bdc360a..9addb648 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -188,6 +188,8 @@ extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np); */ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np); +extern "C" void ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int Np, double Fx, double Fy, double Fz); + /** * \brief compute pressure from D3Q19 distribution * @param dist - D3Q19 distributions diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index 2a85dd40..f50acf2c 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -918,6 +918,51 @@ extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np) { } } +extern "C" void ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int Np, double Fx, double Fy, double Fz) { + int n; + int N = Np; + // distributions + double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9; + double f10, f11, f12, f13, f14, f15, f16, f17, f18; + double rho,vx, vy, vz, rho; + + for (n = 0; n < N; n++) { + //........................................................................ + // Registers to store the distributions + //........................................................................ + f2 = dist[2 * N + n]; + f4 = dist[4 * N + n]; + f6 = dist[6 * N + n]; + f8 = dist[8 * N + n]; + f10 = dist[10 * N + n]; + f12 = dist[12 * N + n]; + f14 = dist[14 * N + n]; + f16 = dist[16 * N + n]; + f18 = dist[18 * N + n]; + //........................................................................ + f0 = dist[n]; + f1 = dist[N + n]; + f3 = dist[3 * N + n]; + f5 = dist[5 * N + n]; + f7 = dist[7 * N + n]; + f9 = dist[9 * N + n]; + f11 = dist[11 * N + n]; + f13 = dist[13 * N + n]; + f15 = dist[15 * N + n]; + f17 = dist[17 * N + n]; + //.................Compute the velocity................................... + rho = f0 + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10 + f11 + f12 + f13 + f14 + f15 + f16 + f17 + f18; + vx = f1 - f2 + f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14; + vy = f3 - f4 + f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18; + vz = f5 - f6 + f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18; + //..................Write the velocity..................................... + vel[n] = vx/rho - Fx*0.5; + vel[N + n] = vy/rho - Fy*0.5; + vel[2 * N + n] = vz/rho - Fz*0.5; + //........................................................................ + } +} + extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *Pressure, int N) { for (int n = 0; n < N; n++) { //........................................................................ diff --git a/cuda/D3Q19.cu b/cuda/D3Q19.cu index f895aff5..c7d8effc 100644 --- a/cuda/D3Q19.cu +++ b/cuda/D3Q19.cu @@ -1593,6 +1593,55 @@ __global__ void dvc_ScaLBL_D3Q19_Momentum(double *dist, double *vel, int N) } } +__global__ void dvc_ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int N, double Fx, double Fy, double Fz) +{ + int n; + // distributions + double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9; + double f10,f11,f12,f13,f14,f15,f16,f17,f18; + double rho,vx,vy,vz; + + int S = N/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(dist, vel, Np, Fx, Fy, Fz); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_Velocity_2nd_order: %s \n",cudaGetErrorString(err)); + } +} + + extern "C" void ScaLBL_D3Q19_Pressure(double *fq, double *Pressure, int Np){ dvc_ScaLBL_D3Q19_Pressure<<< NBLOCKS,NTHREADS >>>(fq, Pressure, Np); } diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 2baf2cb3..41d68cc7 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -328,7 +328,7 @@ void ScaLBL_MRTModel::Run() { //************************************************************************/ if (timestep % ANALYSIS_INTERVAL == 0) { - ScaLBL_D3Q19_Momentum(fq, Velocity, Np); + ScaLBL_D3Q19_Momentum_2nd_order(fq, Velocity, Np, Fx, Fy, Fz); ScaLBL_DeviceBarrier(); comm.barrier(); ScaLBL_Comm->RegularLayout(Map, &Velocity[0], Velocity_x); From cb92d4a518b5c484a602b6efd97d1644b8715cc1 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 15 Oct 2025 15:52:11 -0300 Subject: [PATCH 12/29] Fixing IO:initialize call, third optional parameters should be a bool not a string --- analysis/runAnalysis.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 79ff0834..cc1c6675 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -591,7 +591,7 @@ runAnalysis::runAnalysis(std::shared_ptr input_db, //std::string format = "silo"; format = vis_db->getWithDefault("format", "silo"); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct d_meshData.resize(1); @@ -736,7 +736,7 @@ runAnalysis::runAnalysis(ScaLBL_ColorModel &ColorModel) format = vis_db->getWithDefault("format", "silo"); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct d_meshData.resize(1); From c1e5131ea2fe0c98e0464c9938c6bc72f53619d3 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 15 Oct 2025 16:47:02 -0300 Subject: [PATCH 13/29] Variables positions in visData change depending upon the selection on the db file (#Issue 99) --- analysis/runAnalysis.cpp | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index cc1c6675..e62b33ec 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -248,46 +248,53 @@ class IOWorkItem : public ThreadPool::WorkItemRet { fillHalo fillData(comm.comm, rank_info, n, {1, 1, 1}, 0, 1); + int varCounter = 0; if (vis_db->getWithDefault("save_phase_field", true)) { - ASSERT(visData[0].vars[0]->name == "phase"); - Array &PhaseData = visData[0].vars[0]->data; + ASSERT(visData[0].vars[varCounter]->name == "phase"); + Array &PhaseData = visData[0].vars[varCounter]->data; fillData.copy(Averages.Phi, PhaseData); + varCounter+=1; } if (vis_db->getWithDefault("save_pressure", false)) { - ASSERT(visData[0].vars[1]->name == "Pressure"); - Array &PressData = visData[0].vars[1]->data; + ASSERT(visData[0].vars[varCounter]->name == "Pressure"); + Array &PressData = visData[0].vars[varCounter]->data; fillData.copy(Averages.Pressure, PressData); + varCounter+=1; } if (vis_db->getWithDefault("save_velocity", false)) { - ASSERT(visData[0].vars[2]->name == "Velocity_x"); - ASSERT(visData[0].vars[3]->name == "Velocity_y"); - ASSERT(visData[0].vars[4]->name == "Velocity_z"); - Array &VelxData = visData[0].vars[2]->data; - Array &VelyData = visData[0].vars[3]->data; - Array &VelzData = visData[0].vars[4]->data; + ASSERT(visData[0].vars[varCounter]->name == "Velocity_x"); + ASSERT(visData[0].vars[varCounter+1]->name == "Velocity_y"); + ASSERT(visData[0].vars[varCounter+2]->name == "Velocity_z"); + Array &VelxData = visData[0].vars[varCounter]->data; + Array &VelyData = visData[0].vars[varCounter+1]->data; + Array &VelzData = visData[0].vars[varCounter+2]->data; fillData.copy(Averages.Vel_x, VelxData); fillData.copy(Averages.Vel_y, VelyData); fillData.copy(Averages.Vel_z, VelzData); + varCounter+=3; } if (vis_db->getWithDefault("save_dissipation", false)) { - ASSERT(visData[0].vars[5]->name == "ViscousDissipation"); - Array &ViscousDissipation = visData[0].vars[5]->data; + ASSERT(visData[0].vars[varCounter]->name == "ViscousDissipation"); + Array &ViscousDissipation = visData[0].vars[varCounter]->data; fillData.copy(Averages.Dissipation, ViscousDissipation); + varCounter+1; } if (vis_db->getWithDefault("save_distance", false)) { - ASSERT(visData[0].vars[6]->name == "SignDist"); - Array &SignData = visData[0].vars[6]->data; + ASSERT(visData[0].vars[varCounter]->name == "SignDist"); + Array &SignData = visData[0].vars[varCounter]->data; fillData.copy(Averages.SDs, SignData); + varCounter+=1; } if (vis_db->getWithDefault("save_connected_components", false)) { - ASSERT(visData[0].vars[7]->name == "BlobID"); - Array &BlobData = visData[0].vars[7]->data; + ASSERT(visData[0].vars[varCounter]->name == "BlobID"); + Array &BlobData = visData[0].vars[varCounter]->data; fillData.copy(Averages.morph_n->label, BlobData); + varCounter+=1; } if (vis_db->getWithDefault("write_silo", true)) From 7f6c021794b47b52a2e0f362c9731a2e50aae6ae Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 5 Nov 2025 16:13:05 -0300 Subject: [PATCH 14/29] Small Building error in function ScaLBL_D3Q19_Momentum_2nd_order for cpu --- cpu/D3Q19.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp index f50acf2c..dd37bec6 100644 --- a/cpu/D3Q19.cpp +++ b/cpu/D3Q19.cpp @@ -924,7 +924,7 @@ extern "C" void ScaLBL_D3Q19_Momentum_2nd_order(double *dist, double *vel, int N // distributions double f0, f1, f2, f3, f4, f5, f6, f7, f8, f9; double f10, f11, f12, f13, f14, f15, f16, f17, f18; - double rho,vx, vy, vz, rho; + double rho,vx, vy, vz; for (n = 0; n < N; n++) { //........................................................................ From 31a7c186734c9f1f3f51382ae155fb484cf5f049 Mon Sep 17 00:00:00 2001 From: Jeremy First Date: Fri, 17 Oct 2025 10:40:19 -0500 Subject: [PATCH 15/29] bugfix: restart to initialize from previous timestep --- models/ColorModel.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..cb77352e 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -111,6 +111,9 @@ void ScaLBL_ColorModel::ReadParams(string filename) { if (color_db->keyExists("flux")) { flux = color_db->getScalar("flux"); } + if (color_db->keyExists("timestep")) { + timestep = color_db->getScalar("timestep"); + } inletA = 1.f; inletB = 0.f; outletA = 0.f; From 6db13f4fd2d60485a3f9d6321b3cfd54d77ef1f8 Mon Sep 17 00:00:00 2001 From: jeremyfirst22 <9468245+jeremyfirst22@users.noreply.github.com> Date: Wed, 11 Feb 2026 06:48:44 -0600 Subject: [PATCH 16/29] =?UTF-8?q?bugfix:=20kernel=20packing=20of=20send=20?= =?UTF-8?q?buffers=20must=20complete=20before=20initiatin=E2=80=A6=20(#105?= =?UTF-8?q?)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * bugfix: kernel packing of send buffers must complete before initiating message passing * Reorganize packing and sending of buffers Just reorganize the unpacking, syncronization and send for efficiency as point by JamesEMcClure on PR #105 before merging the code --------- Co-authored-by: Diogo Nardelli Siebert --- common/ScaLBL.cpp | 104 ++++++++++++++++++++++++---------------------- 1 file changed, 55 insertions(+), 49 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index c753be18..284066f0 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -2402,64 +2402,60 @@ void ScaLBL_Communicator::BiSendD3Q7AA(double *Aq, double *Bq) { ScaLBL_D3Q19_Pack(2, dvcSendList_x, sendCount_x, sendCount_x, sendbuf_x, Bq, N); - ScaLBL_DeviceBarrier(); - req1[0] = - MPI_COMM_SCALBL.Isend(sendbuf_x, 2 * sendCount_x, rank_x, sendtag + 0); - req2[0] = - MPI_COMM_SCALBL.Irecv(recvbuf_X, 2 * recvCount_X, rank_X, recvtag + 0); - //...Packing for X face(1,7,9,11,13)................................ ScaLBL_D3Q19_Pack(1, dvcSendList_X, 0, sendCount_X, sendbuf_X, Aq, N); ScaLBL_D3Q19_Pack(1, dvcSendList_X, sendCount_X, sendCount_X, sendbuf_X, Bq, N); - ScaLBL_DeviceBarrier(); - req1[1] = - MPI_COMM_SCALBL.Isend(sendbuf_X, 2 * sendCount_X, rank_X, sendtag + 1); - req2[1] = - MPI_COMM_SCALBL.Irecv(recvbuf_x, 2 * recvCount_x, rank_x, recvtag + 1); - //...Packing for y face(4,8,9,16,18)................................. ScaLBL_D3Q19_Pack(4, dvcSendList_y, 0, sendCount_y, sendbuf_y, Aq, N); ScaLBL_D3Q19_Pack(4, dvcSendList_y, sendCount_y, sendCount_y, sendbuf_y, Bq, N); - ScaLBL_DeviceBarrier(); - req1[2] = - MPI_COMM_SCALBL.Isend(sendbuf_y, 2 * sendCount_y, rank_y, sendtag + 2); - req2[2] = - MPI_COMM_SCALBL.Irecv(recvbuf_Y, 2 * recvCount_Y, rank_Y, recvtag + 2); - - //...Packing for Y face(3,7,10,15,17)................................. + //...Packing for Y face(3,7,10,15,17)................................. ScaLBL_D3Q19_Pack(3, dvcSendList_Y, 0, sendCount_Y, sendbuf_Y, Aq, N); ScaLBL_D3Q19_Pack(3, dvcSendList_Y, sendCount_Y, sendCount_Y, sendbuf_Y, Bq, N); - ScaLBL_DeviceBarrier(); - req1[3] = - MPI_COMM_SCALBL.Isend(sendbuf_Y, 2 * sendCount_Y, rank_Y, sendtag + 3); - req2[3] = - MPI_COMM_SCALBL.Irecv(recvbuf_y, 2 * recvCount_y, rank_y, recvtag + 3); - //...Packing for z face(6,12,13,16,17)................................ ScaLBL_D3Q19_Pack(6, dvcSendList_z, 0, sendCount_z, sendbuf_z, Aq, N); ScaLBL_D3Q19_Pack(6, dvcSendList_z, sendCount_z, sendCount_z, sendbuf_z, Bq, N); - ScaLBL_DeviceBarrier(); - req1[4] = - MPI_COMM_SCALBL.Isend(sendbuf_z, 2 * sendCount_z, rank_z, sendtag + 4); - req2[4] = - MPI_COMM_SCALBL.Irecv(recvbuf_Z, 2 * recvCount_Z, rank_Z, recvtag + 4); - //...Packing for Z face(5,11,14,15,18)................................ ScaLBL_D3Q19_Pack(5, dvcSendList_Z, 0, sendCount_Z, sendbuf_Z, Aq, N); ScaLBL_D3Q19_Pack(5, dvcSendList_Z, sendCount_Z, sendCount_Z, sendbuf_Z, Bq, N); - //................................................................................... + //................................................................................... // Send all the distributions ScaLBL_DeviceBarrier(); + + req1[0] = + MPI_COMM_SCALBL.Isend(sendbuf_x, 2 * sendCount_x, rank_x, sendtag + 0); + req2[0] = + MPI_COMM_SCALBL.Irecv(recvbuf_X, 2 * recvCount_X, rank_X, recvtag + 0); + + req1[1] = + MPI_COMM_SCALBL.Isend(sendbuf_X, 2 * sendCount_X, rank_X, sendtag + 1); + req2[1] = + MPI_COMM_SCALBL.Irecv(recvbuf_x, 2 * recvCount_x, rank_x, recvtag + 1); + + req1[2] = + MPI_COMM_SCALBL.Isend(sendbuf_y, 2 * sendCount_y, rank_y, sendtag + 2); + req2[2] = + MPI_COMM_SCALBL.Irecv(recvbuf_Y, 2 * recvCount_Y, rank_Y, recvtag + 2); + + req1[3] = + MPI_COMM_SCALBL.Isend(sendbuf_Y, 2 * sendCount_Y, rank_Y, sendtag + 3); + req2[3] = + MPI_COMM_SCALBL.Irecv(recvbuf_y, 2 * recvCount_y, rank_y, recvtag + 3); + + req1[4] = + MPI_COMM_SCALBL.Isend(sendbuf_z, 2 * sendCount_z, rank_z, sendtag + 4); + req2[4] = + MPI_COMM_SCALBL.Irecv(recvbuf_Z, 2 * recvCount_Z, rank_Z, recvtag + 4); + req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, 2 * sendCount_Z, rank_Z, sendtag + 5); req2[5] = @@ -2544,46 +2540,54 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for x face(2,8,10,12,14)................................ ScaLBL_D3Q19_Pack(2, dvcSendList_x, 0, sendCount_x, sendbuf_x, &Aq[Component * 7 * N], N); - req1[0] = + + //...Packing for X face(1,7,9,11,13)................................ + ScaLBL_D3Q19_Pack(1, dvcSendList_X, 0, sendCount_X, sendbuf_X, + &Aq[Component * 7 * N], N); + + //...Packing for y face(4,8,9,16,18)................................. + ScaLBL_D3Q19_Pack(4, dvcSendList_y, 0, sendCount_y, sendbuf_y, + &Aq[Component * 7 * N], N); + + //...Packing for Y face(3,7,10,15,17)................................. + ScaLBL_D3Q19_Pack(3, dvcSendList_Y, 0, sendCount_Y, sendbuf_Y, + &Aq[Component * 7 * N], N); + + //...Packing for z face(6,12,13,16,17)................................ + ScaLBL_D3Q19_Pack(6, dvcSendList_z, 0, sendCount_z, sendbuf_z, + &Aq[Component * 7 * N], N); + + //...Packing for Z face(5,11,14,15,18)................................ + ScaLBL_D3Q19_Pack(5, dvcSendList_Z, 0, sendCount_Z, sendbuf_Z, + &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); + + req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x, rank_x, sendtag + 0); req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X, rank_X, recvtag + 0); - //...Packing for X face(1,7,9,11,13)................................ - ScaLBL_D3Q19_Pack(1, dvcSendList_X, 0, sendCount_X, sendbuf_X, - &Aq[Component * 7 * N], N); req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X, rank_X, sendtag + 1); req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x, rank_x, recvtag + 1); - //...Packing for y face(4,8,9,16,18)................................. - ScaLBL_D3Q19_Pack(4, dvcSendList_y, 0, sendCount_y, sendbuf_y, - &Aq[Component * 7 * N], N); req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y, rank_y, sendtag + 2); req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y, rank_Y, recvtag + 2); - //...Packing for Y face(3,7,10,15,17)................................. - ScaLBL_D3Q19_Pack(3, dvcSendList_Y, 0, sendCount_Y, sendbuf_Y, - &Aq[Component * 7 * N], N); req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y, rank_Y, sendtag + 3); req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y, rank_y, recvtag + 3); - //...Packing for z face(6,12,13,16,17)................................ - ScaLBL_D3Q19_Pack(6, dvcSendList_z, 0, sendCount_z, sendbuf_z, - &Aq[Component * 7 * N], N); req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z, rank_z, sendtag + 4); req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z, rank_Z, recvtag + 4); - //...Packing for Z face(5,11,14,15,18)................................ - ScaLBL_D3Q19_Pack(5, dvcSendList_Z, 0, sendCount_Z, sendbuf_Z, - &Aq[Component * 7 * N], N); req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z, rank_Z, sendtag + 5); req2[5] = @@ -2697,6 +2701,7 @@ void ScaLBL_Communicator::TriSendD3Q7AA(double *Aq, double *Bq, double *Cq) { //................................................................................... // Send all the distributions + ScaLBL_DeviceBarrier(); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, 3 * sendCount_x, rank_x, sendtag + 0); req2[0] = @@ -2831,10 +2836,11 @@ void ScaLBL_Communicator::SendHalo(double *data) { ScaLBL_Scalar_Pack(dvcSendList_yZ, sendCount_yZ, sendbuf_yZ, data, N); ScaLBL_Scalar_Pack(dvcSendList_Yz, sendCount_Yz, sendbuf_Yz, data, N); ScaLBL_Scalar_Pack(dvcSendList_YZ, sendCount_YZ, sendbuf_YZ, data, N); + //................................................................................... // Send / Recv all the phase indcator field values //................................................................................... - + ScaLBL_DeviceBarrier(); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x, rank_x, sendtag + 0); req2[0] = From 24d42fab269eb0abb0b3923d0304efbe1a25d047 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Mon, 16 Feb 2026 14:00:14 -0300 Subject: [PATCH 17/29] Updating old github Actions scripts --- .github/workflows/build_guide.yml | 39 +++++++ .github/workflows/c-cpp.yml | 187 +++++++++++++----------------- 2 files changed, 122 insertions(+), 104 deletions(-) create mode 100644 .github/workflows/build_guide.yml diff --git a/.github/workflows/build_guide.yml b/.github/workflows/build_guide.yml new file mode 100644 index 00000000..d132cb08 --- /dev/null +++ b/.github/workflows/build_guide.yml @@ -0,0 +1,39 @@ +name: build_guide + +# Generate sphinx docs + +on: + workflow_dispatch: + +jobs: + + build_docs_job: + runs-on: ubuntu-latest + permissions: + contents: write + steps: + - name: Checkout + uses: actions/checkout@v2.3.4 + + - name: Set up Python + uses: actions/setup-python@v2.2.1 + with: + python-version: 3.9 + + - name: Install dependencies + run: | + python -m pip install -U sphinx + python -m pip install sphinx-rtd-theme + python -m pip install breathe + python -m pip install ghp-import + + - name: make the sphinx docs + run: | + make -C docs clean + make -C docs html + + - name: Run ghp-import + run: | + ghp-import -n -p -f ~/local/doc/build/html + + diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml index 42421825..4fd4f4a4 100644 --- a/.github/workflows/c-cpp.yml +++ b/.github/workflows/c-cpp.yml @@ -1,123 +1,102 @@ name: LBPM CI on: + workflow_dispatch: push: - branches: [ master ] + branches: [master, dev] pull_request: - branches: [ master ] + branches: [master, dev] jobs: build-and-test: - runs-on: ubuntu-latest env: - LBPM_ZLIB_DIR: /home/runner/extlib/zlib - LBPM_HDF5_DIR: /home/runner/extlib/hdf5 - LBPM_SILO_DIR: /home/runner/extlib/silo - MPI_DIR: /home/runner/.openmpi - - steps: - - name: download dependencies - run: | - echo $LBPM_ZLIB_DIR - echo $LBPM_HDF5_DIR - echo $LBPM_SILO_DIR - echo $GITHUB_PATH - echo $GITHUB_WORKSPACE - - sudo apt-get update -y - - wget https://bitbucket.org/AdvancedMultiPhysics/tpl-builder/downloads/Silo-4.10.3RC.modified.tar.gz - wget https://www.zlib.net/fossils/zlib-1.2.11.tar.gz - wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.12/src/hdf5-1.8.12.tar.gz - #wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.10/src/hdf5-1.8.10.tar.gz - - tar -xzvf zlib-1.2.11.tar.gz - tar -xzvf hdf5-1.8.12.tar.gz - tar -xzvf Silo-4.10.3RC.modified.tar.gz - - - - name: check out commit - uses: actions/checkout@v2 - with: - path: LBPM - + MPI_VERSION: "4.1.8" + ZLIB_VERSION: "1.3.1" + HDF5_VERSION: "1.14.6" + MPI_DIR: ${{ github.workspace }}/deps/openmpi + LBPM_ZLIB_DIR: ${{ github.workspace }}/deps/zlib + LBPM_HDF5_DIR: ${{ github.workspace }}/hdf5 - - name: install-openmpi - run: | - wget https://download.open-mpi.org/release/open-mpi/v3.1/openmpi-3.1.2.tar.gz - tar -xvf ./openmpi-3.1.2.tar.gz - ./openmpi-3.1.2/configure --prefix="$HOME/.openmpi" - make -j - sudo make install - echo "$HOME/.openmpi/bin" >> $GITHUB_PATH - - + steps: + - uses: actions/checkout@v4 + with: + path: LBPM - - name: install zlib dependencies - run: | - cd zlib-1.2.11 - ./configure --prefix=$LBPM_ZLIB_DIR - make - sudo make install - cd .. + - name: Install build tools + run: | + sudo apt-get update -y + sudo apt-get install -y build-essential cmake wget tar + - name: Cache compiled deps (OpenMPI, zlib, HDF5) + id: cache-deps + uses: actions/cache@v4 + with: + path: | + deps/openmpi + deps/zlib + deps/hdf5 + key: ${{ runner.os }}-lbpm-deps-mpi${{ env.MPI_VERSION }}-zlib${{ env.ZLIB_VERSION }}-hdf5${{ env.HDF5_VERSION }} + restore-keys: | + ${{ runner.os }}-lbpm-deps- - - name: install hdf5 dependencies - run: | - cd hdf5-1.8.12 - CC=/home/runner/.openmpi/bin/mpicc CXX=/home/runner/.openmpi/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ - ./configure --prefix=$LBPM_HDF5_DIR --enable-parallel --enable-shared --with-zlib=$LBPM_ZLIB_DIR - make - sudo make install - cd .. + - name: Build OpenMPI + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://download.open-mpi.org/release/open-mpi/v4.1/openmpi-${MPI_VERSION}.tar.gz + tar -xzf openmpi-${MPI_VERSION}.tar.gz + cd openmpi-${MPI_VERSION} + ./configure --prefix="${MPI_DIR}" + make -j + make install + - name: Build zlib + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://zlib.net/zlib-${ZLIB_VERSION}.tar.gz + tar -xzf zlib-${ZLIB_VERSION}.tar.gz + cd zlib-${ZLIB_VERSION} + ./configure --prefix="${LBPM_ZLIB_DIR}" + make -j + make install - - name: install silo dependencies - run: | - cd Silo-4.10.3RC.modified - CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ - ./configure --prefix=$LBPM_SILO_DIR -with-hdf5="$LBPM_HDF5_DIR/include,$LBPM_HDF5_DIR/lib" --enable-static - make - sudo make install - cd .. - - - - name: configure cmake - run: | - mkdir build - cd build - rm -rf CMake* - cmake \ - -D CMAKE_BUILD_TYPE:STRING=Release \ - -D CMAKE_C_COMPILER:PATH=$MPI_DIR/bin/mpicc \ - -D CMAKE_CXX_COMPILER:PATH=$MPI_DIR/bin/mpicxx \ - -D MPI_CXX_COMPILER=$MPI_DIR/bin/mpicxx \ - -D CMAKE_C_FLAGS="-fPIC" \ - -D CMAKE_CXX_FLAGS="-fPIC" \ - -D CMAKE_CXX_STD=14 \ - -D TEST_MAX_PROCS=1 \ - -D USE_TIMER=0 \ - -D TIMER_DIRECTORY=$LBPM_TIMER_DIR \ - -D USE_NETCDF=0 \ - -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \ - -D USE_SILO=0 \ - -D HDF5_DIRECTORY=$LBPM_HDF5_DIR \ - -D SILO_DIRECTORY=$LBPM_SILO_DIR \ - -D USE_CUDA=0 \ - $GITHUB_WORKSPACE/LBPM + - name: Build HDF5 (parallel) + if: steps.cache-deps.outputs.cache-hit != 'true' + run: | + wget -q https://sourceforge.net/projects/hdf5.mirror/files/hdf5_${HDF5_VERSION}/hdf5-${HDF5_VERSION}.tar.gz/download -O hdf5-${HDF5_VERSION}.tar.gz + tar -xzf hdf5-${HDF5_VERSION}.tar.gz + cd hdf5-${HDF5_VERSION} + CC="${MPI_DIR}/bin/mpicc" CXX="${MPI_DIR}/bin/mpicxx" CXXFLAGS="-fPIC -O3 -std=c++14" \ + ./configure --prefix="${LBPM_HDF5_DIR}" --enable-parallel --enable-shared --with-zlib="${LBPM_ZLIB_DIR}" + make -j + make install + + # Ensure the cached/restored MPI is on PATH + - name: Add MPI to PATH + run: echo "${MPI_DIR}/bin" >> $GITHUB_PATH - + - name: Configure CMake + run: | + cmake -S LBPM -B build \ + -D USE_SILO=0 \ + -D USE_TIMER=0 \ + -D USE_CUDA=0 \ + -D CMAKE_BUILD_TYPE=Release \ + -D CMAKE_C_COMPILER:PATH="${MPI_DIR}/bin/mpicc" \ + -D CMAKE_CXX_COMPILER:PATH="${MPI_DIR}/bin/mpicxx" \ + -D MPI_CXX_COMPILER="${MPI_DIR}/bin/mpicxx" \ + -D CMAKE_C_FLAGS="-fPIC" \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D CMAKE_CXX_STD=14 \ + -D TEST_MAX_PROCS=1 \ + -D HDF5_DIRECTORY="${LBPM_HDF5_DIR}" - - name: build and make - run: | - cd build - make - sudo make install - cd .. - - - name: tests - run: | - cd build - ctest + - name: build and make + run: | + cd build + make install + + - name: tests + run: | + ctest --test-dir build --output-on-failure From 854302ce2979817b6ce75954bedb96b6802c5fdf Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Mon, 16 Feb 2026 14:01:16 -0300 Subject: [PATCH 18/29] Removing old test script for building openmpi --- .github/workflows/test_install_openmpi.yml | 41 ---------------------- 1 file changed, 41 deletions(-) delete mode 100644 .github/workflows/test_install_openmpi.yml diff --git a/.github/workflows/test_install_openmpi.yml b/.github/workflows/test_install_openmpi.yml deleted file mode 100644 index b193a6e6..00000000 --- a/.github/workflows/test_install_openmpi.yml +++ /dev/null @@ -1,41 +0,0 @@ -name: Install OpenMPI test - -on: - push: - branches: [ master ] - pull_request: - branches: [ master ] - -jobs: - install-openmpi: - - runs-on: ubuntu-latest - - steps: - - name: check path - run: | - echo $PATH - echo $GITHUB_PATH - cmake --version - - - name: download-openmpi - run: wget https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-4.0.2.tar.gz - - - name: extract-openmpi - run: tar -xvf ./openmpi-4.0.2.tar.gz - - - name: configure-openmpi - run: ./openmpi-4.0.2/configure --prefix="/home/${USER}/.openmpi" - - - name: install-openmpi - run: | - make -j - sudo make install - - - name: setting path - run: | - echo "/home/${USER}/.openmpi/bin" >> $GITHUB_PATH - #echo "/home/${USER}/.openmpi/bin" >> $PATH - - name: checking version - run: mpirun --version - From 58d0f7b53e1a8dac825f58ad6e1268cc81b6b9e5 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 18 Feb 2026 13:58:36 -0300 Subject: [PATCH 19/29] VTI (VTK XML Imagedata) Support (#110) I proposed this implementation last year but closed my own PR because it included changes unrelated to this new feature. Since then, we have been using this feature in all our simulations, which has allowed us to identify and fix several small bugs. I believe it is now quite mature, so I am merging it into the dev branch so other developers and users can test it and suggest further improvements. * Adding VTI (VTK XML Imagedata) output support --- IO/Mesh.cpp | 2 + IO/Mesh.h | 2 +- IO/VtiWriter.cpp | 117 ++++++++++ IO/Writer.cpp | 49 +++- IO/Writer.h | 6 +- IO/xmlvtk.cpp | 559 ++++++++++++++++++++++++++++++++++++++++++++ IO/xmlvtk.h | 185 +++++++++++++++ models/MRTModel.cpp | 4 +- 8 files changed, 914 insertions(+), 10 deletions(-) create mode 100755 IO/VtiWriter.cpp create mode 100644 IO/xmlvtk.cpp create mode 100644 IO/xmlvtk.h diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp index 5e5eb96a..0eb9283f 100644 --- a/IO/Mesh.cpp +++ b/IO/Mesh.cpp @@ -598,6 +598,8 @@ std::string getString( FileFormat type ) return "silo"; else if ( type == FileFormat::HDF5 ) return "hdf5"; + else if ( type == FileFormat::VTK ) + return "vtk"; else ERROR( "Invalid type" ); return ""; diff --git a/IO/Mesh.h b/IO/Mesh.h index 9e5f32e6..2ada354c 100644 --- a/IO/Mesh.h +++ b/IO/Mesh.h @@ -24,7 +24,7 @@ enum class VariableType { }; enum class DataType { Double, Float, Int, Null }; enum class MeshType { PointMesh, SurfaceMesh, VolumeMesh, Unknown }; -enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 }; +enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 , VTK}; //! Convert enums to/from strings (more future-proof than static_cast) diff --git a/IO/VtiWriter.cpp b/IO/VtiWriter.cpp new file mode 100755 index 00000000..1e6223e0 --- /dev/null +++ b/IO/VtiWriter.cpp @@ -0,0 +1,117 @@ +#include "IO/HDF5_IO.h" +#include "IO/IOHelpers.h" +#include "IO/MeshDatabase.h" +#include "IO/Writer.h" +#include "IO/silo.h" +#include "IO/xmlvtk.h" +#include "common/MPI.h" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include + +static void writeVti( + const std::string &fullpath, const IO::MeshDataStruct &meshData) +{ + const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + + VTIWriter vti = VTIWriter(std::string(fullpath)); + vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, + ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz); + + vti.setSpacing( 1.0 , 1.0 , 1.0 ); + vti.setOrigin(0,0,0); + vti.setCompress(); + + for ( size_t i = 0; i < meshData.vars.size(); i++ ) + { + const auto &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + vti.addCellData( var.name , "Float64" , "binary" , var.dim , (unsigned char*) var.data.begin() ); + } else if ( var.precision == IO::DataType::Float ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + vti.addCellData( var.name , "Float32" , "binary" , var.dim , (unsigned char*) data2.begin()); + } else if ( var.precision == IO::DataType::Int ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + vti.addCellData( var.name , "Int32" , "binary" , var.dim , (unsigned char*) var.data.begin() ); + } else { + ERROR( "Unsupported format" ); + } + } + + vti.write(); +} + +void writeVtiSummary( + const std::vector &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename ) +{ + const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + PVTIWriter pvti = PVTIWriter( filename ); + int rank = 0; + for ( const auto &data : meshes_written ) + { + for ( const auto &tmp : data.domains ) + { + RankInfoStruct info( rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + char filename[100]; + sprintf( filename, "%05i.vti", rank ); + + VTIWriter vti = VTIWriter( filename ); + vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, + ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz); + vti.setSpacing( 1.0 , 1.0 , 1.0 ); + vti.setOrigin(0,0,0); + vti.setCompress(); + + for ( size_t i = 0; i < meshData.vars.size(); i++ ) + { + const auto &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + vti.addCellData( var.name , "Float64" , "binary" , var.dim , nullptr ); + } else if ( var.precision == IO::DataType::Float ) { + vti.addCellData( var.name , "Float32" , "binary" , var.dim , nullptr ); + } else if ( var.precision == IO::DataType::Int ) { + vti.addCellData( var.name , "Int32" , "binary" , var.dim , nullptr ); + } else { + ERROR( "Unsupported format" ); + } + } + + pvti.addVTIWriter(vti); + rank++; + + } + } + pvti.write(); +} + +std::vector writeMeshesVti( const std::vector &meshData, + const std::string &path, int rank ) +{ + std::vector meshes_written; + char filename[100], fullpath[200]; + sprintf( filename, "%05i.vti", rank ); + sprintf( fullpath, "%s/%s", path.c_str(), filename ); + + for ( size_t i = 0; i < meshData.size(); i++ ) { +// auto mesh = meshData[i].mesh; + auto database = getDatabase( fullpath , meshData[i], IO::FileFormat::VTK, rank ); + + if ( database.meshClass == "DomainMesh" ) { + writeVti( fullpath, meshData[i] ); + } else { + ERROR( "Unknown mesh class or not implemented for vtk/vti output" ); + } + + meshes_written.push_back( database ); + } + return meshes_written; +} + diff --git a/IO/Writer.cpp b/IO/Writer.cpp index 0fb8d135..1ed3cdfa 100644 --- a/IO/Writer.cpp +++ b/IO/Writer.cpp @@ -15,7 +15,7 @@ #include -enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN }; +enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN, VTK }; /**************************************************** @@ -26,7 +26,10 @@ std::vector writeMeshesSilo( void writeSiloSummary( const std::vector &, const std::string & ); std::vector writeMeshesHDF5( const std::vector &, const std::string &, IO::FileFormat, int, Xdmf & ); - +std::vector writeMeshesVti( const std::vector &meshData, + const std::string &path, int rank ); +void writeVtiSummary( + const std::vector &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename ); /**************************************************** * Recursively create the subdirectory * @@ -90,6 +93,8 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap global_IO_format = Format::SILO; else if ( format == "hdf5" ) global_IO_format = Format::HDF5; + else if ( format == "vtk" ) + global_IO_format = Format::VTK; else ERROR( "Unknown format" ); int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank(); @@ -100,9 +105,13 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap filename = global_IO_path + "/summary.LBM"; else if ( global_IO_format == Format::SILO || global_IO_format == Format::HDF5 ) filename = global_IO_path + "/LBM.visit"; + else if ( global_IO_format == Format::VTK) + filename = global_IO_path + "/LBM.pvd"; else ERROR( "Unknown format" ); auto fid = fopen( filename.c_str(), "wb" ); + if ( global_IO_format == Format::VTK) + fprintf( fid, "\n\n \n \n\n" ); fclose( fid ); } } @@ -278,13 +287,12 @@ static std::vector writeMeshesNewFormat( return meshes_written; } - /**************************************************** * Write the mesh data * ****************************************************/ void IO::writeData( const std::string &subdir, const std::vector &meshData, - const Utilities::MPI &comm ) -{ + const Utilities::MPI &comm , int timestep ) +{ if ( global_IO_path.empty() ) IO::initialize(); PROFILE_START( "writeData" ); @@ -310,6 +318,9 @@ void IO::writeData( const std::string &subdir, const std::vector\n \n\n", timestep, subdir.c_str()); + + char *pos = strstr(buffer, " \n\n"); + pos[0] = '\0'; + + fid = fopen( filename.c_str(), "w" ); + fputs(buffer,fid); + fputs(newrow,fid); + fclose(fid); + + free(buffer); + } else { ERROR( "Unknown format" ); } diff --git a/IO/Writer.h b/IO/Writer.h index 3844f3b2..e1d2f8f4 100644 --- a/IO/Writer.h +++ b/IO/Writer.h @@ -38,7 +38,7 @@ void initialize( * @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof) */ void writeData( const std::string &subdir, const std::vector &meshData, - const Utilities::MPI &comm ); + const Utilities::MPI &comm , int timestep = -1); /*! @@ -48,12 +48,14 @@ void writeData( const std::string &subdir, const std::vector * @param[in] meshData The data to write * @param[in] comm The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof) */ + + inline void writeData( int timestep, const std::vector &meshData, const Utilities::MPI &comm ) { char subdir[100]; sprintf( subdir, "vis%03i", timestep ); - writeData( subdir, meshData, comm ); + writeData( subdir, meshData, comm , timestep); } diff --git a/IO/xmlvtk.cpp b/IO/xmlvtk.cpp new file mode 100644 index 00000000..aceb4e30 --- /dev/null +++ b/IO/xmlvtk.cpp @@ -0,0 +1,559 @@ +/* + * Copyright (c) 2025 Diogo Nardelli Siebert + * + * Licensed under either of + * - Apache License, Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0) + * - GNU General Public License, Version 3.0 or later (https://www.gnu.org/licenses/gpl-3.0.html) + * + * SPDX-License-Identifier: (Apache-2.0 OR GPL-3.0-or-later) + */ + +#include "xmlvtk.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +/** @brief Indicates whether data compression is used by default **/ +bool Element::compress = false; +/** @brief Default cache size used during data processing. */ +headerType cacheSize = 10000; + +/** + * @brief Formats a string using printf-style syntax. + * @tparam Args Variadic arguments for formatting. + * @param fmt Format string. + * @return Formatted string. + */ +template +std::string format_string(const char* fmt, Args&&... args) +{ + int size = std::snprintf(nullptr, 0, fmt, std::forward(args)...); + if (size < 0) { + throw std::runtime_error("format_string: snprintf error"); + } + std::vector buf(size + 1); + int size2 = std::snprintf(buf.data(), buf.size(), fmt, std::forward(args)...); + if (size2 < 0) { + throw std::runtime_error("format_string: snprintf error"); + } + + return std::string(buf.data(), buf.data() + size2); +} + +/** + * @brief Encodes binary input to Base64 string representation. + * @param input Input binary data. + * @param len Length of input data. + * @return Encoded Base64 string. + */ +std::string spc_base64_encode(const unsigned char* input, size_t len) +{ + static const char table[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"; + std::string output; + output.reserve((len / 3 + (len % 3 != 0)) * 4); + + for (size_t i = 0; i < len; i += 3) { + uint32_t n = (input[i] << 16) | (i + 1 < len ? input[i + 1] << 8 : 0) | (i + 2 < len ? input[i + 2] : 0); // AAAAAABB BBBBCCCC CCDDDDDD + output.push_back(table[(n >> 18) & 0x3F]); + output.push_back(table[(n >> 12) & 0x3F]); + output.push_back(i + 1 < len ? table[(n >> 6) & 0x3F] : '='); + output.push_back(i + 2 < len ? table[n & 0x3F] : '='); + } + return output; +} + +std::string CellData::header() +{ + return "\n"; +} + +std::string CellData::footer() +{ + return "\n"; +} + +std::string PointData::header() +{ + return "\n"; +} + +std::string PointData::footer() +{ + return "\n"; +} + +std::string AppendData::header() +{ + return "\n"; +} + +std::string AppendData::footer() +{ + return "\n"; +} + +AppendData::AppendData() +{ + this -> totalSize = 0; +} + +unsigned int AppendData::addData(unsigned char* pointer, unsigned int size) +{ + int offset = this->totalSize; + this->sizeList.push_back(size); + this->pointerList.push_back(pointer); + this->totalSize += (size+4); + return offset; +} + +DataArray::DataArray(const std::string& name_,const std::string& type_, const std::string& format_, int components_, unsigned char* pointer_, uint64_t points_) + : components(components_), name(name_), type(type_), format(format_), pointer(pointer_), points(points_) + +{ + if ( (this -> type == "Int8") || (this -> type == "UInt8") ) this -> typeSize = 1; + else if ( (this -> type == "Int16") || (this -> type == "UInt16") ) this -> typeSize = 2; + else if ( (this -> type == "Int32") || (this -> type == "UInt32") || (this -> type == "Float32") ) this -> typeSize = 4; + else if ( (this -> type == "Int64") || (this -> type == "UInt64") || (this -> type == "Float64") ) this -> typeSize = 8; + this -> dataSize = components * points; +} + +std::string DataArray::header() +{ + std::ostringstream stringStream; + stringStream << format_string("", offset).c_str() : ">"); + return stringStream.str(); +} + +std::string DataArray::footer() +{ + return "\n\n"; +} + +void VTIWriter::setCompress() +{ + this->compress = true; +} + +PVTIWriter::PVTIWriter(const std::string& filename_) + : filename(filename_), + originX(0), originY(0), originZ(0), + wholeMinX(0), wholeMinY(0), wholeMinZ(0), wholeMaxX(0), wholeMaxY(0), wholeMaxZ(0), pieceCounter(0) +{ +} + +VTIWriter::VTIWriter(const std::string& filename_) + : filename(filename_) , vtkVersion("0.1"), + originX(0), originY(0), originZ(0) +{ + this -> compress = false; + if (sizeof(headerType) == 8) this -> headerTypeName = "UInt64"; + else if (sizeof(headerType) == 4) this -> headerTypeName = "UInt32"; + this -> byteOrder = "LittleEndian"; +} + +void VTIWriter::addPointData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer) +{ + std::unique_ptr data = std::make_unique( name , type , format , components , (unsigned char*) pointer , pd.points ); + data -> compress = this -> compress; + + if (format == "appended") + { + data -> offset = ad.addData( data->pointer , data->dataSize * data->typeSize ); + ad.compress = this->compress; + } + + pd.addChild( std::move(data) ); +} + +void VTIWriter::addCellData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer) +{ + std::unique_ptr data = std::make_unique( name , type , format , components , (unsigned char*) pointer , cd.cells ); + data -> compress = this -> compress; + + if (format == "appended") + { + data -> offset = ad.addData( data->pointer , data->dataSize * data->typeSize ); + ad.compress = this->compress; + } + + cd.addChild( std::move(data) ); +} + +std::string VTIWriter::header() +{ + std::ostringstream stringStream; + stringStream << "" << endl; + stringStream << format_string("", "1.0", byteOrder.c_str() , headerTypeName.c_str() , compress ? "compressor=\"vtkZLibDataCompressor\"" : "" ) << endl; + stringStream << format_string("", spaceX, spaceY, spaceZ) << endl; + stringStream << format_string("" , pieceMinX , pieceMaxX , pieceMinY , pieceMaxY, pieceMinZ , pieceMaxZ) << endl; + return stringStream.str(); +} + +std::string PVTIWriter::header() +{ + std::ostringstream stringStream; + stringStream << "" << endl; + stringStream << "" << endl; + stringStream << format_string("", spaceX, spaceY, spaceZ) << endl; + + return stringStream.str(); +} + +std::string VTIWriter::footer() +{ + return "\n"; +} + +std::string PVTIWriter::footer() +{ + return "\n\n"; +} + +void VTIWriter::setWholeExtent(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_) +{ + wholeMinX = minX_; + wholeMinY = minY_; + wholeMinZ = minZ_; + wholeMaxX = maxX_; + wholeMaxY = maxY_; + wholeMaxZ = maxZ_; + setPiece(minX_,minY_,minZ_,maxX_,maxY_, maxZ_); +} + +std::ostream& operator<<(std::ostream& os, AppendData& obj) +{ + if (obj.totalSize > 0) + { + os << obj.header(); + os << "_"; + for (int k = 0; k < obj.pointerList.size(); k++) + { + if (obj.compress == false) + { + os.write( (char*) &obj.sizeList[k], sizeof(headerType) ); + os.write( (char*) obj.pointerList[k], obj.sizeList[k] ); + } + else + { + streampos beginPos = os.tellp(); + + headerType totalByteSize = obj.sizeList[k]; + headerType numberOfBlocks = totalByteSize / cacheSize + (totalByteSize % cacheSize != 0 ); + + size_t infoSize = sizeof(headerType) * (3 + numberOfBlocks); + std::vectorcompressedInfo(infoSize ); + + compressedInfo[0] = numberOfBlocks; + compressedInfo[1] = (numberOfBlocks > 1) ? cacheSize : totalByteSize; + compressedInfo[2] = (totalByteSize % cacheSize == 0) ? cacheSize : totalByteSize % cacheSize ; + + os.write((char*) compressedInfo.data() , infoSize ); + + for (int n = 1; n <= numberOfBlocks; n++) + { + int numberOfBytesInBlock = (n < numberOfBlocks) ? compressedInfo[1] : compressedInfo[2]; + + uLongf compressedLength = compressBound(numberOfBytesInBlock); + std::vectorcompressedData (compressedLength); + + const Bytef* src = reinterpret_cast( obj.pointerList[k] + (n - 1) * cacheSize ); + int ret = compress( (Bytef *) (compressedData.data() ), &compressedLength, src, numberOfBytesInBlock ); + if (ret != Z_OK) + { + std::cerr << "Error compressing data , code =" << ret << std::endl; + throw std::runtime_error("Failed compressing data using zlib in AppendData"); + } + compressedInfo[n+2] = compressedLength; + + os.write( reinterpret_cast( compressedData.data() ) , compressedLength ); + } + + streampos endPos = os.tellp(); + + os.seekp( beginPos ); + os.write( (char*) compressedInfo.data() , infoSize ); + os.seekp( endPos ); + } + } + os << endl; + os << obj.footer(); + + } + return os; +} + +std::ostream& operator<<(std::ostream& os, DataArray& obj) +{ + os << obj.header() << endl; + + headerType totalByteSize = static_cast( (obj.dataSize) * (obj.typeSize) ); + headerType count = min( 12 - sizeof(headerType) , totalByteSize ); + + if (obj.format == "binary") + { + if (obj.compress == false) + { + os << spc_base64_encode( (unsigned char *) &totalByteSize, sizeof(headerType) ); + count = 0; + + int leftOverSize = 0; + unsigned char leftOverBuffer[3]; + + for (; count < totalByteSize ; count += cacheSize) + { + unsigned char* pointer = obj.pointer + count; + size_t size = min( cacheSize, totalByteSize - count); + if (leftOverSize > 0) + { + for (; leftOverSize < 3; leftOverSize++) + { + if (size > 0) + { + leftOverBuffer[leftOverSize] = *(pointer++); + size--; + } + else break; + } + + os << spc_base64_encode( leftOverBuffer , leftOverSize ); + leftOverSize = 0; + } + + leftOverSize = (size % 3); + size -= leftOverSize; + for (int i = 0; i < leftOverSize; i++) + { + leftOverBuffer[i] = pointer[size + i]; + } + + os << spc_base64_encode( pointer , size ); + + } + + if (leftOverSize > 0) + { + os << spc_base64_encode( leftOverBuffer , leftOverSize ); + leftOverSize = 0; + } + } + else + { + streampos beginPos = os.tellp(); + + std::vectorcompressedData (cacheSize + 3); + int leftOver = 0; + + headerType numberOfBlocks = totalByteSize / cacheSize + (totalByteSize % cacheSize != 0 ); + size_t infoSize = sizeof(headerType) * (3 + numberOfBlocks); + std::vectorcompressedInfo(infoSize ); + + compressedInfo[0] = numberOfBlocks; + compressedInfo[1] = (numberOfBlocks > 1) ? cacheSize : totalByteSize; + compressedInfo[2] = (totalByteSize % cacheSize == 0) ? cacheSize : totalByteSize % cacheSize ; + + os << spc_base64_encode( reinterpret_cast(compressedInfo.data()) , infoSize ); + + for (int n = 1; n <= numberOfBlocks; n++) + { + uLongf compressedLength = cacheSize; + int numberOfBytesInBlock = (n < numberOfBlocks) ? compressedInfo[1] : compressedInfo[2]; + int ret = compress( reinterpret_cast(compressedData.data() + 3) , &compressedLength, obj.pointer + (n-1) * cacheSize, numberOfBytesInBlock ); + if (ret != Z_OK) + { + std::cerr << "Error compressing data , code =" << ret << std::endl; + throw std::runtime_error("Failed compressing data using zlib in AppendData"); + } + compressedInfo[n+2] = compressedLength; + + int encodeSize = 3* ( (compressedLength + leftOver)/3 ); + os << spc_base64_encode( compressedData.data() + 3 - leftOver , encodeSize ); + + int newLeftOver = compressedLength + leftOver - encodeSize ; + if (newLeftOver > 0) + { + memcpy( compressedData.data() + 3 - newLeftOver, compressedData.data() + 3 - leftOver + encodeSize, newLeftOver); + } + leftOver = newLeftOver; + } + + if (leftOver > 0) + { + os << spc_base64_encode(compressedData.data() + 3 - leftOver , leftOver ); + } + + streampos endPos = os.tellp(); + os.seekp( beginPos ); + + os << spc_base64_encode( reinterpret_cast(compressedInfo.data()) , infoSize ); + os.seekp( endPos ); + } + + os << obj.footer(); + } + + return os; +} + +void VTIWriter::setPiece(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_) +{ + pieceMinX = minX_; + pieceMinY = minY_; + pieceMinZ = minZ_; + pieceMaxX = maxX_; + pieceMaxY = maxY_; + pieceMaxZ = maxZ_; + sizeX = maxX_ - minX_ + 1; + sizeY = maxY_ - minY_ + 1; + sizeZ = maxZ_ - minZ_ + 1; + pd.points = sizeX * sizeY * sizeZ; + cd.cells = (sizeX-1) * (sizeY-1) * (sizeZ-1); +} + +void VTIWriter::setOrigin(double x_,double y_, double z_) +{ + originX = x_; + originY = y_; + originZ = z_; +} + +void VTIWriter::setSpacing(double sx_,double sy_, double sz_) +{ + spaceX = sx_; + spaceY = sy_; + spaceZ = sz_; +} + +void PVTIWriter::write() +{ + std::ofstream file; + file.open(this -> filename); + file << header(); + + for (int i = 0; i< cellDataName.size(); i++) + { + if (i==0) file << "" << endl; + file << format_string( "" ,cellDataName[i].c_str() ,cellDataComponents[i] ,cellDataType[i].c_str() ) << endl; + if (i== cellDataName.size()-1) file << "" << endl; + } + + for (int i = 0; i< pointDataName.size(); i++) + { + if (i==0) file << "" << endl; + file << format_string( "" ,pointDataName[i].c_str() ,pointDataComponents[i] ,pointDataType[i].c_str() ) << endl; + if (i== pointDataName.size()-1) file << "" << endl; + } + + + for (int i = 0; i< pieceFilename.size(); i++) + { + file << format_string("", pieceMinX[i], pieceMaxX[i], pieceMinY[i] , pieceMaxY[i] , pieceMinZ[i] , pieceMaxZ[i] , pieceFilename[i].c_str() ) << endl; + } + + file << footer(); + file.close(); +} + +void VTIWriter::write() +{ + file.open(this -> filename); + file << header(); + + if (pd.sizeChild() > 0) + { + file << pd.header(); + for (int n = 0; n < pd.sizeChild() ; n++) + { + DataArray* array = (DataArray*) pd.getChild(n); + file << *array; + } + file << pd.footer(); + } + + if (cd.sizeChild() > 0) + { + file << cd.header(); + for (int n = 0; n < cd.sizeChild() ; n++) + { + DataArray* array = (DataArray*) cd.getChild(n); + file << *array; + } + file << cd.footer(); + } + + file << "" << endl; + file << "" << endl; + + file << ad; + file << footer(); + file.close(); +} + +void PVTIWriter::addVTIWriter(VTIWriter& write) +{ + pieceFilename.push_back( write.filename ); + + pieceMaxX.push_back( write.pieceMaxX ); + pieceMaxY.push_back( write.pieceMaxY ); + pieceMaxZ.push_back( write.pieceMaxZ ); + pieceMinX.push_back( write.pieceMinX ); + pieceMinY.push_back( write.pieceMinY ); + pieceMinZ.push_back( write.pieceMinZ ); + + if (pieceFilename.size() == 1) + { + for (int n = 0; n < write.pd.sizeChild() ; n++) + { + DataArray* data = (DataArray*) write.pd.getChild(n); + pointDataName.push_back( data -> name ); + pointDataType.push_back( data -> type); + pointDataComponents.push_back( data -> components ); + } + + for (int n = 0; n < write.cd.sizeChild() ; n++) + { + DataArray* data = (DataArray*) write.cd.getChild(n); + cellDataName.push_back( data -> name ); + cellDataType.push_back( data -> type); + cellDataComponents.push_back( data -> components ); + } + + originX = write.originX; + originY = write.originX; + originZ = write.originX; + + spaceX = write.spaceX; + spaceY = write.spaceY; + spaceZ = write.spaceZ; + + wholeMaxX = write.wholeMaxX; + wholeMaxY = write.wholeMaxY; + wholeMaxZ = write.wholeMaxZ; + wholeMinX = write.wholeMinX; + wholeMinY = write.wholeMinY; + wholeMinZ = write.wholeMinZ; + } + else + { + wholeMaxX = max(wholeMaxX,write.wholeMaxX); + wholeMaxY = max(wholeMaxY,write.wholeMaxY); + wholeMaxZ = max(wholeMaxZ,write.wholeMaxZ); + wholeMinX = min(wholeMinX,write.wholeMinX); + wholeMinY = min(wholeMinY,write.wholeMinY); + wholeMinZ = min(wholeMinZ,write.wholeMinZ); + } + + +} diff --git a/IO/xmlvtk.h b/IO/xmlvtk.h new file mode 100644 index 00000000..9538eda4 --- /dev/null +++ b/IO/xmlvtk.h @@ -0,0 +1,185 @@ +/* + * Copyright (c) 2025 Diogo Nardelli Siebert + * + * Licensed under either of + * - Apache License, Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0) + * - GNU General Public License, Version 3.0 or later (https://www.gnu.org/licenses/gpl-3.0.html) + * + * SPDX-License-Identifier: (Apache-2.0 OR GPL-3.0-or-later) + */ + +#ifndef __XMLVTK_H_INCLUDED__ +#define __XMLVTK_H_INCLUDED__ + +#include +#include +#include +#include +#include + +#define headerType u_int64_t + + +/** + * @brief Encodes binary data into Base64 format. + * @param input Pointer to the input data. + * @param len Length of the input data. + * @return Base64-encoded string. + */ +std::string spc_base64_encode(const unsigned char* input, size_t len); + +/** + * @class Element + * @brief Represents a hierarchical XML-like element supporting compression and VTK serialization. + */ +class Element +{ + public: + static bool compress; +/** @brief Returns the parent element. */ + Element* getParent() { return parent ;} +/** @brief Returns the nth child element. *//** ... */ + Element* getChild(int n) { return this->child[n].get(); } +/** @brief Adds a new child element and sets its parent. */ + void addChild(std::unique_ptr e) { + e -> parent = this; + child.push_back(std::move( e ) ); + } + +/** @brief Returns the number of child elements. */ + int sizeChild( ) { return child.size(); } +/** @brief Generates a header string for VTK output. */ + virtual std::string header() = 0; +/** @brief Generates a footer string for VTK output. */ + virtual std::string footer() = 0; + private: + Element* parent; + std::vector< std::unique_ptr > child; +}; + +class PointData: public Element +{ + public: + std::string header(); + std::string footer(); + unsigned int points; +}; + +class CellData: public Element +{ + public: + std::string header(); + std::string footer(); + unsigned int cells; +}; + +class DataArray : public Element +{ + public: + DataArray(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer, uint64_t size); + unsigned char* pointer; + std::string header(); + std::string footer(); + void write( std::ofstream& file ); + friend std::ostream& operator<<(std::ostream& os, DataArray& obj); + uint64_t typeSize; + uint64_t dataSize ; + uint64_t offset; + uint64_t points; + std::string name; + std::string format; + std::string type; + int components; + private: + int mode = 0; + + +}; + +class AppendData : public Element +{ + public: + AppendData(); + std::string header(); + std::string footer(); + unsigned int addData(unsigned char* pointer, unsigned int size); + friend std::ostream& operator<<(std::ostream& os, AppendData& obj); + private: + std::vector pointerList; + std::vector sizeList; + headerType totalSize; +}; + +class VTIWriter : public Element +{ + public: + VTIWriter(const std::string& filename); + std::string footer(); + std::string header(); + + void write(); + void setWholeExtent(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_); + void setPiece(int64_t minX_,int64_t minY_, int64_t minZ_,int64_t maxX_,int64_t maxY_, int64_t maxZ_); + void setOrigin(double x_,double y_, double z_); + void setSpacing(double sx_,double sy_, double sz_); + void addPointData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer); + void addCellData(const std::string& name,const std::string& type, const std::string& format, int components, unsigned char* pointer); + void setCompress(); + + std::string filename; + std::string vtkVersion; /*!< The Vtk File Format Version of the file */ + std::string fileTitle; /*!< The title of the file (do not confuse with the name of the file) */ + std::string dataSetType; /*!< The type of geometry (grid) that data is associeted to (STRUCTURED GRID for LBM applications) */ + std::string headerTypeName; + + int64_t wholeMinX, wholeMinY, wholeMinZ; + int64_t wholeMaxX, wholeMaxY, wholeMaxZ; + int64_t pieceMinX, pieceMinY, pieceMinZ; + int64_t pieceMaxX, pieceMaxY, pieceMaxZ; + int64_t sizeX, sizeY, sizeZ; /*!< Lenght in pixels of the image in each axis */ + + /*!< Lenght in pixels of the image in each axis */ + double spaceX, spaceY, spaceZ; /*!< Ratio of the different axis */ + double originX,originY,originZ; /*!< Position of the origin of the image */ + + std::string byteOrder; /*!< Position of the origin of the image */ + std::ofstream file; + CellData cd; + PointData pd; + AppendData ad; +}; + +class PVTIWriter : public Element +{ + public: + PVTIWriter(const std::string& filename); + std::string filename; + std::string footer(); + std::string header(); + + double spaceX, spaceY, spaceZ; + double originX,originY,originZ; + + int pieceCounter; + + int64_t wholeMinX, wholeMinY, wholeMinZ; + int64_t wholeMaxX, wholeMaxY, wholeMaxZ; + + std::vector pieceMinX, pieceMinY, pieceMinZ; + std::vector pieceMaxX, pieceMaxY, pieceMaxZ; + std::vector pieceFilename; + + std::vector cellDataComponents; + std::vector cellDataName; + std::vector cellDataType; + + std::vector pointDataComponents; + std::vector pointDataName; + std::vector pointDataType; + + void addVTIWriter(VTIWriter& write); + void write(); +}; + + +#endif diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 41d68cc7..674132f1 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -483,7 +483,7 @@ void ScaLBL_MRTModel::VelocityField() { auto VzVar = std::make_shared(); auto SignDistVar = std::make_shared(); - IO::initialize("", format, "false"); + IO::initialize("", format, false); // Create the MeshDataStruct visData.resize(1); visData[0].meshName = "domain"; @@ -526,7 +526,7 @@ void ScaLBL_MRTModel::VelocityField() { fillData.copy(Velocity_x, VelxData); fillData.copy(Velocity_y, VelyData); fillData.copy(Velocity_z, VelzData); - + IO::writeData(timestep, visData, Dm->Comm); } } From 1432725c55ce6d7078a266ff5dd08bd8005994e7 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 18 Feb 2026 14:20:52 -0300 Subject: [PATCH 20/29] Fix building error on newer MPI/GCC versions --- tests/DataAggregator.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/DataAggregator.cpp b/tests/DataAggregator.cpp index 137ebd61..466b25a9 100644 --- a/tests/DataAggregator.cpp +++ b/tests/DataAggregator.cpp @@ -5,6 +5,7 @@ #include #include #include +#include using namespace std; From 8f77aed7a71fb68b3f2cf4c204a079bcb9bb2de6 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Wed, 18 Feb 2026 14:21:10 -0300 Subject: [PATCH 21/29] Fix small bugs on the CTest scripts --- tests/CMakeLists.txt | 2 +- tests/test_MPI.cpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ceafff4..64e53260 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -84,7 +84,7 @@ ADD_LBPM_TEST_1_2_4( TestBlobIdentify ) #ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 ) #ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 ) ADD_LBPM_TEST_PARALLEL( TestSegDist 8 ) -ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 ) +ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 test.db) ADD_LBPM_TEST_1_2_4( testCommunication ) ADD_LBPM_TEST( TestWriter ) ADD_LBPM_TEST( TestDatabase ) diff --git a/tests/test_MPI.cpp b/tests/test_MPI.cpp index c6d12011..f14ba219 100644 --- a/tests/test_MPI.cpp +++ b/tests/test_MPI.cpp @@ -1652,7 +1652,9 @@ int main( int argc, char *argv[] ) } // Test GPU aware MPI +#ifdef USE_CUDA test_GPU_aware( &ut ); +#endif } // Limit the scope so objects are destroyed From 6711a7c4b035ee5f58fa734003179b73acbf65d4 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Tue, 21 Apr 2026 10:46:37 -0300 Subject: [PATCH 22/29] Update ZLIB_VERSION to 1.3.2 in workflow --- .github/workflows/c-cpp.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml index 4fd4f4a4..098decaf 100644 --- a/.github/workflows/c-cpp.yml +++ b/.github/workflows/c-cpp.yml @@ -12,7 +12,7 @@ jobs: runs-on: ubuntu-latest env: MPI_VERSION: "4.1.8" - ZLIB_VERSION: "1.3.1" + ZLIB_VERSION: "1.3.2" HDF5_VERSION: "1.14.6" MPI_DIR: ${{ github.workspace }}/deps/openmpi From 5c2fbfef6c45423413c97591a828eb91f26669db Mon Sep 17 00:00:00 2001 From: Ricardo Leite Martin Bazarin Date: Tue, 21 Apr 2026 12:06:04 -0300 Subject: [PATCH 23/29] spurious diffusive effect correction (#112) * spurious diffusive effect correction * missing ID_solid * Trying to implement new wettability imposition using bitwise operations for optimization * Fix building errors * Fix building errors * Fix building error * Fix small building errors * Try to remove unnecessary variables * Fix unexcepted error when computing normal vectors * Another attempt to fix Subphase analysis error * Optimized Implementation of the wettability fix * Adding wettability fix (optimized) to CPU implementation * Fix sign * reading communication as a solid * plus/minus addition for contact angle * cleaning * Reverting changes in CMakeLists.txt * Updating CMakeLists.txt with dev branch to avoid merge issues * Update c-cpp.yml --------- Co-authored-by: Diogo Nardelli Siebert Co-authored-by: Diogo Nardelli Siebert --- common/ScaLBL.h | 4 +- cpu/Color.cpp | 549 ++++++++++++++++++++++++++++++-------- cuda/Color.cu | 601 ++++++++++++++++++++++++++++++++---------- models/ColorModel.cpp | 57 +++- models/ColorModel.h | 1 + 5 files changed, 948 insertions(+), 264 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 9addb648..521bda87 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -628,7 +628,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *neighborList, double *dist, * @param Np - size of local sub-domain (derived from Domain structure) */ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); @@ -660,7 +660,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( */ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *NeighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 9e32147c..15d84f78 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -17,6 +17,7 @@ #include #define STOKES +#define f64_eps 1E-12 extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, @@ -1429,7 +1430,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, // double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int* NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -1446,6 +1447,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + double nspx, nspy, nspz; // Rock Fluid interface normal vector const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -1476,79 +1478,239 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); - // Get the 1D index based on regular data layout - ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk - 1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk + 1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk - strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk + strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk - strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk + strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk - strideY - 1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk + strideY + 1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk + strideY - 1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk - strideY + 1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk - strideZ - 1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk + strideZ + 1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk + strideZ - 1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk - strideZ + 1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk - strideZ - strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk + strideZ + strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk + strideZ - strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk - strideZ + strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 - //............Compute the Color Gradient................................... - nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); - ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); - nz = -(m5 - m6 + 0.5 * (m11 - m12 - m13 + m14 + m15 - m16 - m17 + m18)); - - //...........Normalize the Color Gradient................................. - C = sqrt(nx * nx + ny * ny + nz * nz); - double ColorMag = C; - if (C == 0.0) - ColorMag = 1.0; - nx = nx / ColorMag; - ny = ny / ColorMag; - nz = nz / ColorMag; + ijk = Map[n]; + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + + //............Compute the Color Gradient................................... + nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); + + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } + // q=0 fq = dist[n]; rho = fq; @@ -1989,7 +2151,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2005,7 +2167,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 2 // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2020,7 +2182,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; @@ -2041,7 +2203,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *neighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, unsigned int* NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -2061,6 +2223,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + double nspx, nspy, nspz; const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -2091,78 +2254,238 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( rlx_setA = 1.f / tau; rlx_setB = 8.f * (2.f - rlx_setA) / (8.f - rlx_setA); - // Get the 1D index based on regular data layout + // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ + nn = ijk - 1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ + nn = ijk + 1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ + nn = ijk - strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ + nn = ijk + strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ + nn = ijk - strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ + nn = ijk + strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ + nn = ijk - strideY - 1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ + nn = ijk + strideY + 1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ + nn = ijk + strideY - 1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ + nn = ijk - strideY + 1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ + nn = ijk - strideZ - 1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ + nn = ijk + strideZ + 1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ + nn = ijk + strideZ - 1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ + nn = ijk - strideZ + 1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ + nn = ijk - strideZ - strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ + nn = ijk + strideZ + strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ + nn = ijk + strideZ - strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ + nn = ijk - strideZ + strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... - nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); - ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); - nz = -(m5 - m6 + 0.5 * (m11 - m12 - m13 + m14 + m15 - m16 - m17 + m18)); + nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); //...........Normalize the Color Gradient................................. - C = sqrt(nx * nx + ny * ny + nz * nz); - double ColorMag = C; - if (C == 0.0) - ColorMag = 1.0; - nx = nx / ColorMag; - ny = ny / ColorMag; - nz = nz / ColorMag; + C = sqrt(nx*nx+ny*ny+nz*nz); + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -2666,7 +2989,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2685,7 +3008,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2705,7 +3028,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * nspz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; diff --git a/cuda/Color.cu b/cuda/Color.cu index fc35c5c9..d4e8b90f 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -20,6 +20,7 @@ #define NBLOCKS 1024 #define NTHREADS 256 +#define f64_eps 1E-12 __global__ void dvc_ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz) { @@ -1266,7 +1267,7 @@ __global__ void dvc_ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int } -__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -1281,6 +1282,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + double nspx, nspy, nspz; // const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1317,63 +1319,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk-1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk+1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk-strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk+strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk-strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk+strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk-strideY-1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk+strideY+1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk+strideY-1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk-strideY+1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk-strideZ-1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk+strideZ+1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk+strideZ-1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk-strideZ+1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk-strideZ-strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk+strideZ+strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk+strideZ-strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk-strideZ+strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1381,11 +1381,174 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //...........Normalize the Color Gradient................................. C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -1660,7 +1823,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -1802,7 +1964,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*nspx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -1817,7 +1979,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 2 // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*nspy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -1831,7 +1993,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*nspz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -1850,7 +2012,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, unsigned int *NeighborSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int n,nn,ijk,nread; @@ -1869,6 +2031,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + double nspx, nspy, nspz; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1904,63 +2067,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double // Get the 1D index based on regular data layout ijk = Map[n]; - // COMPUTE THE COLOR GRADIENT - //........................................................................ - //.................Read Phase Indicator Values............................ - //........................................................................ - nn = ijk-1; // neighbor index (get convention) - m1 = Phi[nn]; // get neighbor for phi - 1 - //........................................................................ - nn = ijk+1; // neighbor index (get convention) - m2 = Phi[nn]; // get neighbor for phi - 2 - //........................................................................ - nn = ijk-strideY; // neighbor index (get convention) - m3 = Phi[nn]; // get neighbor for phi - 3 - //........................................................................ - nn = ijk+strideY; // neighbor index (get convention) - m4 = Phi[nn]; // get neighbor for phi - 4 - //........................................................................ - nn = ijk-strideZ; // neighbor index (get convention) - m5 = Phi[nn]; // get neighbor for phi - 5 - //........................................................................ - nn = ijk+strideZ; // neighbor index (get convention) - m6 = Phi[nn]; // get neighbor for phi - 6 - //........................................................................ - nn = ijk-strideY-1; // neighbor index (get convention) - m7 = Phi[nn]; // get neighbor for phi - 7 - //........................................................................ - nn = ijk+strideY+1; // neighbor index (get convention) - m8 = Phi[nn]; // get neighbor for phi - 8 - //........................................................................ - nn = ijk+strideY-1; // neighbor index (get convention) - m9 = Phi[nn]; // get neighbor for phi - 9 - //........................................................................ - nn = ijk-strideY+1; // neighbor index (get convention) - m10 = Phi[nn]; // get neighbor for phi - 10 - //........................................................................ - nn = ijk-strideZ-1; // neighbor index (get convention) - m11 = Phi[nn]; // get neighbor for phi - 11 - //........................................................................ - nn = ijk+strideZ+1; // neighbor index (get convention) - m12 = Phi[nn]; // get neighbor for phi - 12 - //........................................................................ - nn = ijk+strideZ-1; // neighbor index (get convention) - m13 = Phi[nn]; // get neighbor for phi - 13 - //........................................................................ - nn = ijk-strideZ+1; // neighbor index (get convention) - m14 = Phi[nn]; // get neighbor for phi - 14 - //........................................................................ - nn = ijk-strideZ-strideY; // neighbor index (get convention) - m15 = Phi[nn]; // get neighbor for phi - 15 - //........................................................................ - nn = ijk+strideZ+strideY; // neighbor index (get convention) - m16 = Phi[nn]; // get neighbor for phi - 16 - //........................................................................ - nn = ijk+strideZ-strideY; // neighbor index (get convention) - m17 = Phi[nn]; // get neighbor for phi - 17 - //........................................................................ - nn = ijk-strideZ+strideY; // neighbor index (get convention) - m18 = Phi[nn]; // get neighbor for phi - 18 + + nn = ijk - 1; // neighbor index (get convention) + m1 = Phi[nn]; // get neighbor for phi - 1 + + nn = ijk + 1; // neighbor index (get convention) + m2 = Phi[nn]; // get neighbor for phi - 2 + + nn = ijk - strideY; // neighbor index (get convention) + m3 = Phi[nn]; // get neighbor for phi - 3 + + nn = ijk + strideY; // neighbor index (get convention) + m4 = Phi[nn]; // get neighbor for phi - 4 + + nn = ijk - strideZ; // neighbor index (get convention) + m5 = Phi[nn]; // get neighbor for phi - 5 + + nn = ijk + strideZ; // neighbor index (get convention) + m6 = Phi[nn]; // get neighbor for phi - 6 + + nn = ijk - strideY - 1; // neighbor index (get convention) + m7 = Phi[nn]; // get neighbor for phi - 7 + + nn = ijk + strideY + 1; // neighbor index (get convention) + m8 = Phi[nn]; // get neighbor for phi - 8 + + nn = ijk + strideY - 1; // neighbor index (get convention) + m9 = Phi[nn]; // get neighbor for phi - 9 + + nn = ijk - strideY + 1; // neighbor index (get convention) + m10 = Phi[nn]; // get neighbor for phi - 10 + + nn = ijk - strideZ - 1; // neighbor index (get convention) + m11 = Phi[nn]; // get neighbor for phi - 11 + + nn = ijk + strideZ + 1; // neighbor index (get convention) + m12 = Phi[nn]; // get neighbor for phi - 12 + + nn = ijk + strideZ - 1; // neighbor index (get convention) + m13 = Phi[nn]; // get neighbor for phi - 13 + + nn = ijk - strideZ + 1; // neighbor index (get convention) + m14 = Phi[nn]; // get neighbor for phi - 14 + + nn = ijk - strideZ - strideY; // neighbor index (get convention) + m15 = Phi[nn]; // get neighbor for phi - 15 + + nn = ijk + strideZ + strideY; // neighbor index (get convention) + m16 = Phi[nn]; // get neighbor for phi - 16 + + nn = ijk + strideZ - strideY; // neighbor index (get convention) + m17 = Phi[nn]; // get neighbor for phi - 17 + + nn = ijk - strideZ + strideY; // neighbor index (get convention) + m18 = Phi[nn]; // get neighbor for phi - 18 + //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1968,11 +2129,174 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //...........Normalize the Color Gradient................................. C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + if (C > f64_eps) + { + nx = nx/C; + ny = ny/C; + nz = nz/C; + } + //...........Correct wettability vector for Mass Balance................................. + + unsigned int data = NeighborSolid[n]; + if ( (data != 0) ) + { + char isNeighborSolid; + char countSolid = 0; + + nspx = 0; + nspy = 0; + nspz = 0; + + isNeighborSolid = (data >> 1) & 1u; + countSolid += isNeighborSolid; + m1 = isNeighborSolid * m1; + nspx += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 2) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m2; + nspx -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 3) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m3; + nspy += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 4) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m4; + nspy -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 5) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m5; + nspz += 2 * isNeighborSolid; + + isNeighborSolid = (data >> 6) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m6; + nspz -= 2 * isNeighborSolid; + + isNeighborSolid = (data >> 7) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m7; + nspx += isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 8) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m8; + nspx -= isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 9) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m9; + nspx += isNeighborSolid; + nspy -= isNeighborSolid; + + isNeighborSolid = (data >> 10) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m10; + nspx -= isNeighborSolid; + nspy += isNeighborSolid; + + isNeighborSolid = (data >> 11) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m11; + nspx += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 12) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m12; + nspx -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 13) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m13; + nspx += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 14) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m14; + nspx -= isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 15) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m15; + nspy += isNeighborSolid; + nspz += isNeighborSolid; + + isNeighborSolid = (data >> 16) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m16; + nspy -= isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 17) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m17; + nspy += isNeighborSolid; + nspz -= isNeighborSolid; + + isNeighborSolid = (data >> 18) & 1u; + countSolid += isNeighborSolid; + m1 += isNeighborSolid * m18; + nspy -= isNeighborSolid; + nspz += isNeighborSolid; + + m3 = sqrt( nspx * nspx + nspy * nspy + nspz * nspz ); + + if (m3 == 0.0) + m3 = 1.0f; + + nspx = -nspx / m3; + nspy = -nspy / m3; + nspz = -nspz / m3; + + m1 = m1 / countSolid; + m2 = (nx*nspx + ny*nspy + nz*nspz); + + m3 = 1.0f-m2*m2; + m3 = (m3 > 0.0f) ? sqrtf(m3) : 1.0f; + + double nspxp = (nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspyp = (ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzp = (nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double nspxm = -(nx - nspx*m2)*sqrt(1.0f-m1*m1)/m3 + nspx*m1; + double nspym = -(ny - nspy*m2)*sqrt(1.0f-m1*m1)/m3 + nspy*m1; + double nspzm = -(nz - nspz*m2)*sqrt(1.0f-m1*m1)/m3 + nspz*m1; + + double dotp = nx*nspxp + ny*nspyp + nz*nspzp; + double dotm = nx*nspxm + ny*nspym + nz*nspzm; + + if (dotp > dotm) { + nspx = nspxp; + nspy = nspyp; + nspz = nspzp; + } else { + nspx = nspxm; + nspy = nspym; + nspz = nspzm; + } + } + else + { + nspx = nx; + nspy = ny; + nspz = nz; + } + + if (C < f64_eps) + { + nspx = nspy = nspz = 0.0; + nx = ny = nz = 0; + } // q=0 fq = dist[n]; @@ -2298,7 +2622,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -19*alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -2451,7 +2774,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*nspx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2469,7 +2792,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*nspy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2488,7 +2811,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*nspz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2877,7 +3200,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_ColorMomentum(int *neighborList, double * //........................................................................ //..............carry out relaxation process.............................. //..........Toelke, Fruediger et. al. 2006................................ - if (C == 0.0) nx = ny = nz = 0.0; m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) -alpha*C - m1); m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4); @@ -4007,14 +4329,14 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A } // Pressure Boundary Conditions Functions -extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB, + dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, NeighborSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -4025,13 +4347,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do } extern "C" void ScaLBL_D3Q19_AAodd_Color(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, unsigned int *NeighborSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Vel, + dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, NeighborSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); @@ -4164,4 +4486,3 @@ extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Sour dvc_ScaLBL_CopySlice_z<<>>(Phi,Nx,Ny,Nz,Source,Dest); } - diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..ed643512 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -32,7 +32,7 @@ ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), inletA(0), inletB(0), outletA(0), outletB(0), Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), - BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), + BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), NeighborSolid(nullptr), NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), Bq(nullptr), Den(nullptr), Phi(nullptr), ColorGrad(nullptr), Velocity(nullptr), Pressure(nullptr), comm(COMM) { @@ -448,6 +448,8 @@ void ScaLBL_ColorModel::Create() { ScaLBL_AllocateDeviceMemory((void **)&Pressure, sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&Velocity, 3 * sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&ColorGrad, 3 * sizeof(double) * Np); + ScaLBL_AllocateDeviceMemory((void **)&NeighborSolid, sizeof(unsigned int) * Np); + //........................................................................... // Update GPU data structures if (rank == 0) @@ -503,6 +505,43 @@ void ScaLBL_ColorModel::Create() { if (rank == 0) printf("Model created \n"); delete[] PhaseLabel; + + unsigned int *TmpSolid = new unsigned int[Np]; + + for (int k = 1; k < Nz - 1; k++) { + for (int j = 1; j < Ny - 1; j++) { + for (int i = 1; i < Nx - 1; i++) { + int idx = Map(i, j, k); + unsigned int data = 0; + if (!(idx < 0)) { + if ((Map(i - 1 , j, k) == -1) ) data |= (1u << 1); + if ((Map(i + 1 , j, k) == -1) ) data |= (1u << 2); + if ((Map(i , j - 1, k) == -1) ) data |= (1u << 3); + if ((Map(i , j + 1, k) == -1) ) data |= (1u << 4); + if ((Map(i , j, k - 1) == -1) ) data |= (1u << 5); + if ((Map(i , j, k + 1) == -1) ) data |= (1u << 6); + if ((Map(i - 1 , j - 1, k) == -1) ) data |= (1u << 7); + if ((Map(i + 1 , j + 1, k) == -1) ) data |= (1u << 8); + if ((Map(i - 1 , j + 1, k) == -1) ) data |= (1u << 9); + if ((Map(i + 1, j - 1, k) == -1) ) data |= (1u << 10); + if ((Map(i - 1 , j, k - 1) == -1) ) data |= (1u << 11); + if ((Map(i + 1 , j, k + 1) == -1) ) data |= (1u << 12); + if ((Map(i - 1 , j, k + 1) == -1) ) data |= (1u << 13); + if ((Map(i + 1 , j, k - 1) == -1) ) data |= (1u << 14); + if ((Map(i , j - 1, k - 1) == -1) ) data |= (1u << 15); + if ((Map(i , j + 1, k + 1) == -1) ) data |= (1u << 16); + if ((Map(i , j - 1, k + 1) == -1) ) data |= (1u << 17); + if ((Map(i , j + 1, k - 1) == -1) ) data |= (1u << 18); + TmpSolid[idx] = data; + } + } + } + } + + ScaLBL_CopyToDevice(NeighborSolid, TmpSolid, sizeof(unsigned int) * Np); + ScaLBL_Comm->Barrier(); + delete[] TmpSolid; + } /******************************************************** @@ -690,7 +729,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -709,7 +748,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -735,7 +774,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -754,7 +793,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); @@ -1142,7 +1181,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -1161,7 +1200,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -1187,7 +1226,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1206,7 +1245,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, NeighborSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); diff --git a/models/ColorModel.h b/models/ColorModel.h index 888f4fc2..b8a1dea5 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -138,6 +138,7 @@ class ScaLBL_ColorModel { IntArray Map; signed char *id; + unsigned int *NeighborSolid; int *NeighborList; int *dvcMap; double *fq, *Aq, *Bq; From bcf9cf89f9931484d931a573976b7f01baa8033f Mon Sep 17 00:00:00 2001 From: jeremyfirst22 <9468245+jeremyfirst22@users.noreply.github.com> Date: Tue, 21 Apr 2026 10:24:10 -0500 Subject: [PATCH 24/29] bugfix: min_steady_timesteps in lbpm_color_simulator (#107) I am merging this PR since the bug is clear, and the additional changes made by @jeremyfirst22 guarantee the same behavior for the default parameters. I also agree that dividing by "analysis_interval" makes more sense in terms of the relative error per time step. However, the key point here is that this should be well documented so that the user understands what the tolerance means. --- models/ColorModel.cpp | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 338fdc7e..ff37ea11 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -672,7 +672,6 @@ double ScaLBL_ColorModel::Run(int returntime) { bool RESCALE_FORCE = false; bool SET_CAPILLARY_NUMBER = false; bool TRIGGER_FORCE_RESCALE = false; - double tolerance = 0.01; auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); auto current_db = db->cloneDatabase(); auto flow_db = db->getDatabase("FlowAdaptor"); @@ -698,9 +697,8 @@ double ScaLBL_ColorModel::Run(int returntime) { color_db->getScalar("rescale_force_after_timestep"); RESCALE_FORCE = true; } - if (analysis_db->keyExists("tolerance")) { - tolerance = analysis_db->getScalar("tolerance"); - } + double tolerance = analysis_db->getWithDefault("tolerance", 1e-5); + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map); @@ -806,7 +804,7 @@ double ScaLBL_ColorModel::Run(int returntime) { Den); // allow initial ramp-up to get closer to steady state CURRENT_TIMESTEP += 2; - if (CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { + if (CURRENT_TIMESTEP % analysis_interval == 0 && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { analysis.finish(); double volB = Averages->gwb.V; @@ -843,7 +841,7 @@ double ScaLBL_ColorModel::Run(int returntime) { fabs(muA * flow_rate_A + muB * flow_rate_B) / (5.796 * alpha); bool isSteady = false; - if ((fabs((Ca - Ca_previous) / Ca) < tolerance && + if ((fabs((Ca - Ca_previous) / analysis_interval / Ca) < tolerance && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS)) isSteady = true; if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS) @@ -1113,7 +1111,12 @@ double ScaLBL_ColorModel::Run(int returntime) { printf("Ca = %f, (previous = %f) \n", Ca, Ca_previous); } } + + break; // steady-state achieved, exit. } + + // save for convergence checks + Ca_previous = Ca; } } analysis.finish(); @@ -1142,11 +1145,9 @@ double ScaLBL_ColorModel::Run(int returntime) { void ScaLBL_ColorModel::Run() { int nprocs = nprocx * nprocy * nprocz; const RankInfoStruct rank_info(rank, nprocx, nprocy, nprocz); - int analysis_interval = - 1000; // number of timesteps in between in situ analysis - if (analysis_db->keyExists("analysis_interval")) { - analysis_interval = analysis_db->getScalar("analysis_interval"); - } + + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); + double tolerance = analysis_db->getWithDefault("tolerance", 0.0); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); From 48db887d653b3a09b2ccad1c5a4baa26b9b9d035 Mon Sep 17 00:00:00 2001 From: jeremyfirst22 <9468245+jeremyfirst22@users.noreply.github.com> Date: Tue, 21 Apr 2026 12:07:39 -0500 Subject: [PATCH 25/29] 'sw_steady' protocol for color model (#108) * feature: 'sw_steady' protocol for color model to exit explicitly on sw convergence * Separate handling for centrifuge and sw_steady Refactor protocol handling to separate the centrifuge and sw_steady cases, making them consistent with the existing parameter reading structure and enabling a custom error message. --------- Co-authored-by: Diogo Nardelli Siebert --- .../models/color/protocols/sw_steady.rst | 72 +++++++++++++++++++ models/ColorModel.cpp | 34 +++++++-- tests/lbpm_color_simulator.cpp | 4 +- 3 files changed, 102 insertions(+), 8 deletions(-) create mode 100644 docs/source/userGuide/models/color/protocols/sw_steady.rst diff --git a/docs/source/userGuide/models/color/protocols/sw_steady.rst b/docs/source/userGuide/models/color/protocols/sw_steady.rst new file mode 100644 index 00000000..d63d066a --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/sw_steady.rst @@ -0,0 +1,72 @@ +====================================== +Color model -- Sw Steady +====================================== + +The water saturation steady state protocol is identical to the centrifuge protocol, with +the exception that the simulation explicity converges on the saturation state of fluid +A. + +That is, the simulation exits when + +.. math:: + :nowrap: + + $$ + \frac{\left | S_{w, i+1} - S_{w, i} \right |}{S_{w, i}} \le \epsilon + $$ + +averaged over the ``analysis_interval`` where :math:`S_{w,i}` is the saturation of fluid +A at step :math:`i` and :math:`\epsilon` is an allowed convergence threshold, or when +the ``timestepMax`` is reached, whichever occurs first. + +By default, :math:`\epsilon` is set to zero (i.e., the simulation continues through the +``timestepMax``), but can be set via ``tolerance`` within the ``Analysis`` section of +the input file database as shown below. + + +.. code-block:: c + + Color { + protocol = "sw_steady" + timestepMax = 1000000 // maximum timtestep + alpha = 0.005 // controls interfacial tension + rhoA = 1.0 // controls the density of fluid A + rhoB = 1.0 // controls the density of fluid B + tauA = 0.7 // controls the viscosity of fluid A + tauB = 0.7 // controls the viscosity of fluid B + F = 0, 0, -1.0e-5 // body force + din = 1.0 // inlet density (controls pressure) + dout = 1.0 // outlet density (controls pressure) + WettingConvention = "SCAL" // convention for sign of wetting affinity + ComponentLabels = 0, -1, -2 // image labels for solid voxels + ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label + Restart = false + } + Domain { + Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw" + ReadType = "8bit" // data type + N = 900, 900, 1600 // size of original image + nproc = 2, 2, 2 // process grid + n = 200, 200, 200 // sub-domain size + offset = 300, 300, 300 // offset to read sub-domain + voxel_length = 1.66 // voxel length (in microns) + ReadValues = -2, -1, 0, 1, 2 // labels within the original image + WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM + BC = 3 // boundary condition type (0 for periodic) + } + Analysis { + analysis_interval = 1000 // logging interval for timelog.csv + subphase_analysis_interval = 5000 // loggging interval for subphase.csv + visualization_interval = 100000 // interval to write visualization files + N_threads = 4 // number of analysis threads (GPU version only) + restart_interval = 1000000 // interval to write restart file + restart_file = "Restart" // base name of restart file + tolerance = 1e-9 // Sw convergence tolerance + } + Visualization { + write_silo = true // write SILO databases with assigned variables + save_8bit_raw = true // write labeled 8-bit binary files with phase assignments + save_phase_field = true // save phase field within SILO database + save_pressure = false // save pressure field within SILO database + save_velocity = false // save velocity field within SILO database + } diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index ff37ea11..aedc4795 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -183,7 +183,15 @@ void ScaLBL_ColorModel::ReadParams(string filename) { "pressure boundary condition \n"); } domain_db->putScalar("BC", BoundaryCondition); - } else if (protocol == "core flooding") { + } else if (protocol == "sw_steady") { + if (BoundaryCondition != 3) { + BoundaryCondition = 3; + if (rank == 0) + printf("WARNING: protocol (sw_steady) supports only constant " + "pressure boundary condition \n"); + } + domain_db->putScalar("BC", BoundaryCondition); + } else if (protocol == "core flooding") { if (rank == 0) printf("Using core flooding protocol \n"); if (BoundaryCondition != 4) { @@ -1159,7 +1167,11 @@ void ScaLBL_ColorModel::Run() { Map); //analysis.createThreads( analysis_method, 4 ); auto t1 = std::chrono::system_clock::now(); - while (timestep < timestepMax) { + + double delta_sw = 1.0; + double sw_prev = -1.0; + + while (timestep < timestepMax && delta_sw > tolerance) { PROFILE_START("Update"); // *************ODD TIMESTEP************* @@ -1256,13 +1268,23 @@ void ScaLBL_ColorModel::Run() { //************************************************************************ PROFILE_STOP("Update"); - if (rank == 0 && timestep % analysis_interval == 0 && - BoundaryCondition == 4) { - printf("%i %f \n", timestep, din); - } // Run the analysis analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den); + + if (timestep % analysis_interval == 0){ + analysis.finish(); + + double volA = Averages->gnb.V / Dm->Volume; + double volB = Averages->gwb.V / Dm->Volume; + double sw = volB / (volA + volB); + + delta_sw = fabs(sw - sw_prev) / analysis_interval / sw; + if (rank == 0) + printf("t: %d sw: %0.5e dSw/dt: %.5e\n", timestep, sw, delta_sw); + + sw_prev = sw; + } } analysis.finish(); PROFILE_STOP("Loop"); diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index cafa7a52..977d304b 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -66,7 +66,8 @@ int main( int argc, char **argv ) // structure and allocate variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables - if (SimulationMode == "legacy"){ + auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); + if (PROTOCOL == "sw_steady") { ColorModel.Run(); } else { @@ -75,7 +76,6 @@ int main( int argc, char **argv ) bool ContinueSimulation = true; /* Variables for simulation protocols */ - auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); /* image sequence protocol */ int IMAGE_INDEX = 0; int IMAGE_COUNT = 0; From 2999c65abd4650b45ee9e3a5eb386e272054121b Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Tue, 21 Apr 2026 15:57:30 -0300 Subject: [PATCH 26/29] Small bug where the maximum distance was computed using sumReduce instead of maxReduce, leading to the wrong maximum pore size in the sample which was also dependent on the number of tasks used. (#113) Co-authored-by: Diogo Nardelli Siebert --- analysis/morphology.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 728bd6ed..9667740c 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -508,7 +508,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, // total Global is the number of nodes in the pore-space totalGlobal = Dm->Comm.sumReduce(count); - maxdistGlobal = Dm->Comm.sumReduce(maxdist); + maxdistGlobal = Dm->Comm.maxReduce(maxdist); double volume = double(nprocx * nprocy * nprocz) * double(nx - 2) * double(ny - 2) * double(nz - 2); double volume_fraction = totalGlobal / volume; From 326c0d84389fcde622637619fcf497bf20ca9e5e Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Mon, 4 May 2026 16:12:42 -0300 Subject: [PATCH 27/29] Updating the documentation to include the VTK output option --- docs/source/userGuide/IO/fileformat.rst | 4 +++- docs/source/userGuide/visualization/vtk.rst | 19 +++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) create mode 100644 docs/source/userGuide/visualization/vtk.rst diff --git a/docs/source/userGuide/IO/fileformat.rst b/docs/source/userGuide/IO/fileformat.rst index c1da21cf..012a8088 100644 --- a/docs/source/userGuide/IO/fileformat.rst +++ b/docs/source/userGuide/IO/fileformat.rst @@ -2,12 +2,14 @@ I/O conventions for LBPM ======================== -There are three main kinds of output file that are supported by LBPM. +There are four main kinds of output file that are supported by LBPM. * CSV files -- space-delimited CSV files are used by the internal analysis framework * formatted binary files -- SILO and HDF5 formats are supported for visualization data + +* formatted binary files -- XML VTK format is also supported for visualization data * unformatted binary files -- ``.raw`` extension diff --git a/docs/source/userGuide/visualization/vtk.rst b/docs/source/userGuide/visualization/vtk.rst new file mode 100644 index 00000000..69a519f9 --- /dev/null +++ b/docs/source/userGuide/visualization/vtk.rst @@ -0,0 +1,19 @@ +====================================== +VTK XML output format +====================================== + +Due to several limitations in the support for VisIt files in ParaView, we added support for VTK XML (.vti) files, which are natively supported by ParaView. + +This output format can be enabled by setting format = "vtk" in the Visualization section of the database file. The fields to be written to output can be chosen using the same parameters as in the VisIt output. + +.. code:: c + + Visualization { + format = "vtk" + save_phase_field = true + save_pressure = false + } + +The interval for writing output files is defined by the ``visualization_interval`` parameter in the ``Analysis`` section. + +LBPM also generates and updates a ``LBM.pvd`` file, which allows ParaView to open all .vti files generated during the current simulation. From c0415d939141e8c5a69e986863cf9aec73b731f3 Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Mon, 4 May 2026 21:56:30 -0300 Subject: [PATCH 28/29] A new morphology simulator, **`fmm_micp_drainage`**, has been added. It introduces two protocols: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * **MICP**, where the non-wetting phase invades from all faces of the sample; * **Drainage**, supporting both compressible and incompressible cases. Additionally, new implementations of the **Euclidean Distance Transform (EDT)** and **component labeling** have been developed. This is an initial version and still requires a parallel implementation. The goal is to eventually replace the current morphology functions in LBPM. To achieve this, the following developments are needed: 1. Full parallel MPI implementation of EDT and component labeling; 2. Adaptation of the EDT algorithm to compute distances from solid–fluid interfaces (walls); 3. Addition of a signed version of the EDT algorithm; 4. An option to compute EDT and component labeling with periodic boundary conditions. In the future, the **morphDrain** pre-processor could be unified, allowing users to select the desired method through a protocol option. For now, these functions and the simulator are provided as separate files to avoid confusion with existing implementations until all features listed above are fully implemented and tested. --- docs/source/examples/morphology/micp.rst | 29 + tests/CMakeLists.txt | 2 +- tests/fmm_micp_drainage.cpp | 727 +++++++++++++++++++++++ 3 files changed, 757 insertions(+), 1 deletion(-) create mode 100644 docs/source/examples/morphology/micp.rst create mode 100644 tests/fmm_micp_drainage.cpp diff --git a/docs/source/examples/morphology/micp.rst b/docs/source/examples/morphology/micp.rst new file mode 100644 index 00000000..70ef0f36 --- /dev/null +++ b/docs/source/examples/morphology/micp.rst @@ -0,0 +1,29 @@ +************************************ +MICP Morphology Simulator +************************************ + +To expand and improve the current morphological tools in LBPM, new +morphology pre-processors and simulators were added. These changes +encompass the implementation of state-of-the-art morphology algorithms, +as well as an additional MICP simulator in which the invading non-wetting +phase is injected through all faces of the sample. + +.. code:: c + + Domain { + Filename = "crop_bt_101x112x88_uint8.raw" + N = 101, 112, 88 // domain size + n = 101, 112, 88 + nproc = 1, 1, 1 + ReadValues = 0, 1, 2 + WriteValues = 0, 1, 2 + voxel_length = 1 + BC = 0 + } + + FM { + Diameters = 1, 30, 1 // start, end, step + SaveImage = true + direction = "surround" // axis of intrusion + protocol = "micp" + } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 64e53260..f8a31b60 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -48,7 +48,7 @@ ADD_LBPM_EXECUTABLE( TestIonModel ) ADD_LBPM_EXECUTABLE( TestNernstPlanck ) ADD_LBPM_EXECUTABLE( TestPNP_Stokes ) ADD_LBPM_EXECUTABLE( TestMixedGrad ) - +ADD_LBPM_EXECUTABLE( fmm_micp_drainage ) CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY ) diff --git a/tests/fmm_micp_drainage.cpp b/tests/fmm_micp_drainage.cpp new file mode 100644 index 00000000..ff37b057 --- /dev/null +++ b/tests/fmm_micp_drainage.cpp @@ -0,0 +1,727 @@ +/* +<<<<<<< HEAD + Copyright 2026 + Diogo Nardelli Siebert, Universidade Federal de Santa Catarina + Bernardo Gehlen, Universidade Federal de Santa Catarina + Alexandre Miers Zabot, Universidade Federal de Santa Catarina + + 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 :contentReference[oaicite:0]{index=0}, 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 . +*/ + +/* + ---------------------------------------------------------------------- + Methodological Notes + + This file implements the method described in: + + A. M. Zabot et al., + "A Unified Algorithm for the Young–Laplace Method Applied to Porous Media," + Brazilian Journal of Physics, vol. 54, no. 3, 2024, p. 63. + + The distance transform is computed using the method from: + + P. F. Felzenszwalb and D. P. Huttenlocher, + "Distance transforms of sampled functions," + Theory of Computing, vol. 8, no. 1, pp. 415–428, 2012. + + Connected component labeling is based on: + + L. He, Y. Chao, and K. Suzuki, + "A run-based two-scan labeling algorithm," + IEEE Transactions on Image Processing, vol. 17, no. 5, pp. 749–756, 2008. + ---------------------------------------------------------------------- +*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include "../common/Array.h" +#include "../common/Domain.h" +#include "../common/UtilityMacros.h" + +#define SOLID ((unsigned char)0) +#define DISPLACED ((unsigned char)1) +#define INJECTED ((unsigned char)2) + +#define BACKGROUND 0 +#define FOREGROUND 1 + +using namespace std; + +void merge(const int &u, const int &v, vector &next, vector &tail, + vector &rtable) { + for (int i = v; i != -1;) { + rtable[i] = u; + i = next[i]; + } + next[tail[u]] = v; + tail[u] = tail[v]; +} + +void resolve(const int &x, const int &y, vector &next, vector &tail, + vector &rtable) { + const int u = rtable[x]; + const int v = rtable[y]; + if (u < v) + merge(u, v, next, tail, rtable); + else if (v < u) + merge(v, u, next, tail, rtable); +} + +void component_labeling(IntArray &IMG, const int &F, const int &B) { + + size_t maxNumberOfLabels = (IMG.length() + 1) / 2; + + std::vector next(maxNumberOfLabels); + std::vector tail(maxNumberOfLabels); + std::vector rtable(maxNumberOfLabels); + + const int nx = IMG.size(0); + const int ny = IMG.size(1); + const int nz = IMG.size(2); + + int lx = 0, nl = 1; + vector uniq_labels(3); + int nuniq; + + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + + if (IMG(x, y, z) == F) { + const int lq = (x > 0) ? IMG(x - 1, y, z) : B; + const int lp = (y > 0) ? IMG(x, y - 1, z) : B; + const int lz = (z > 0) ? IMG(x, y, z - 1) : B; + + nuniq = 0; + if (lp != B) { + uniq_labels[nuniq] = lp; + nuniq++; + } + if (lq != B && lq != lp) { + uniq_labels[nuniq] = lq; + nuniq++; + } + if (lz != B && lz != lp && lz != lq) { + uniq_labels[nuniq] = lz; + nuniq++; + } + + // Deals with non unique labels in neighboaring points by unifying them + switch (nuniq) { + + case 0: + nl++; + lx = nl; + + rtable[nl] = nl; + next[nl] = -1; + tail[nl] = nl; + break; + + case 1: + lx = uniq_labels[0]; + break; + + case 2: + + resolve(uniq_labels[0], uniq_labels[1], next, tail, + rtable); + + lx = min(uniq_labels[0], uniq_labels[1]); + break; + + case 3: + + resolve(uniq_labels[0], uniq_labels[1], next, tail, + rtable); + resolve(uniq_labels[0], uniq_labels[2], next, tail, + rtable); + resolve(uniq_labels[1], uniq_labels[2], next, tail, + rtable); + + lx = min(uniq_labels[0], + min(uniq_labels[1], uniq_labels[2])); + break; + } + + IMG(x, y, z) = lx; + } + } + } + } + + int *img = IMG.data(); + + for (size_t n = 0; n < IMG.length(); n++) { + if (*img != B) + *img = rtable[*img]; + img++; + } +} + +// Compute the intersection between two lower-envelope parabolas. +static inline float intersection(int q, int vk, float fq, float fvk) { + float qq = (float)q * (float)q; + float vv = (float)vk * (float)vk; + return ((fq + qq) - (fvk + vv)) / (2.0 * ((float)q - (float)vk)); +} + +// Compute the exact squared Euclidean distance transform for a 1D line. +void edt_1d(const int *f, int *g, int n) { + int *v = (int *)malloc((size_t)n * sizeof(int)); + float *z = (float *)malloc((size_t)(n + 1) * sizeof(float)); + + int k = 0; + v[0] = 0; + z[0] = -INFINITY; + z[1] = INFINITY; + + for (int q = 1; q < n; q++) { + float s; + + while (1) { + int vk = v[k]; + s = intersection(q, vk, f[q], f[vk]); + + if (s <= z[k]) { + k--; + if (k < 0) { + k = 0; + break; + } + } else { + break; + } + } + + if (k == 0) { + int vk = v[k]; + s = intersection(q, vk, f[q], f[vk]); + + if (s <= z[k]) { + v[0] = q; + z[0] = -INFINITY; + z[1] = INFINITY; + continue; + } + } + + k++; + v[k] = q; + z[k] = s; + z[k + 1] = INFINITY; + } + + k = 0; + for (int x = 0; x < n; x++) { + while (z[k + 1] < (float)x) { + k++; + } + + int vk = v[k]; + int dx = x - vk; + g[x] = dx * dx + f[vk]; + } + + free(v); + free(z); +} + +// Copy one strided line, transform it in 1D, and write it back. +void process_line(int *edt2, const size_t first, const size_t stride, int n) { + + IntArray f(n); + IntArray g(n); + + size_t pos = first; + for (int i = 0; i < n; i++) { + f(i) = edt2[pos]; + pos += stride; + } + + edt_1d(f.data(), g.data(), n); + + pos = first; + for (int i = 0; i < n; i++) { + edt2[pos] = g(i); + pos += stride; + } +} + +template +void edt_3d(unsigned char target, Array &image, IntArray &distance2) { + const int nx = image.size(0); + const int ny = image.size(1); + const int nz = image.size(2); + + size_t nvox = image.length(); + int BIG = static_cast(nx * nx + ny * ny + nz * nz) + 1; + + TYPE *img = image.data(); + int *edt2 = distance2.data(); + + // Initialize target voxels with zero distance and all others with a large value. + for (int i = 0; i < (int)nvox; i++) { + edt2[i] = (img[i] == target) ? 0 : BIG; + } + + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + process_line(edt2, nx * (y + z * ny), 1, nx); + } + } + + for (int z = 0; z < nz; z++) { + for (int x = 0; x < nx; x++) { + process_line(edt2, z * nx * ny + x, nx, ny); + } + } + + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + process_line(edt2, y * nx + x, nx * ny, nz); + } + } +} + +// Validate database options and report the accepted values on failure. +void checkOption(std::string a, std::vector s, std::string keyName) { + std::string message = "Error: Invalid option '" + a + "' for " + keyName + + ". Valid options are: "; + for (string b : s) { + if (a == b) + return; + message += "'" + b + "', "; + } + message.pop_back(); + ERROR(message); +} + +template +void setRegion(Array &A, TYPE value, int x0, int x1, int y0, int y1, + int z0, int z1) { + for (int z = z0; z < z1; z++) + for (int y = y0; y < y1; y++) + for (int x = x0; x < x1; x++) { + A(x, y, z) = value; + } +} + +class Full_Morphology { +public: + Full_Morphology(int, char *[]); + int calc(const int &); + +public: + std::vector r_ini = {0, 0, 0}; + std::vector r_end = {0, 0, 0}; + + std::vector flowAxis = { + false, false, false}; // Injection occurs along the selected axis + bool flowPos = false; // Sense of invasion + + int ny, nx, nz; // Work dimensions (with reservoirs, if any) + int dimy, dimx, dimz; // Original dimensions + int NP; // Porous pixels + int inletPos, outletPos; // Reservoir regions + + int rChamberI[3]; + int rChamberO[3]; + + double resolution; + + vector diameter; // Invasion diameters + + bool compressible = false; // Compressibility, enabled for MICP + + bool allFaces; // Add reservoirs around all active faces for MICP + bool saveImg; + + Array originalState; // Original image, used for reference + Array currentState; // Work image, modified during the simulation + Array finalMap; + IntArray originalEDT; + Array trapped; +}; + +Full_Morphology::Full_Morphology(int argc, char *argv[]) { + + string filename; + + filename = argv[1]; + + auto db = std::make_shared(filename); + + auto domain_db = db->getDatabase("Domain"); + auto fm_db = db->getDatabase("FM"); + + auto size = domain_db->getVector("N"); + nx = size[0]; + ny = size[1]; + nz = size[2]; + + finalMap.resize(size[0], size[1], size[2]); + finalMap.fill(-1); + + auto ReadValues = domain_db->getVector("ReadValues"); + auto WriteValues = domain_db->getVector("WriteValues"); + + resolution = domain_db->getScalar("voxel_length"); + + auto READFILE = domain_db->getScalar("Filename"); + const string mmfile(READFILE); + + saveImg = fm_db->getScalar("SaveImage"); + + auto protocol = fm_db->getScalar("protocol"); + checkOption(protocol, {"micp", "drainage"}, "protocol"); + + if (protocol == "micp") { + compressible = true; + allFaces = true; + flowPos = true; + if (size[2] > 1) + flowAxis[2] = true; + else if (size[1] > 1) + flowAxis[1] = true; + else if (size[0] > 1) + flowAxis[0] = true; + } else if (protocol == "drainage") { + allFaces = false; + auto direction = fm_db->getScalar("direction"); + checkOption(direction, {"+x", "-x", "+y", "-y", "+z", "-z"}, + "direction"); + flowAxis[0] = (direction[1] == 'x'); + flowAxis[1] = (direction[1] == 'y'); + flowAxis[2] = (direction[1] == 'z'); + flowPos = (direction[0] == '+'); + compressible = fm_db->getWithDefault("compressible", false); + } + + auto diameterRange = fm_db->getVector("Diameters"); + int numberOfDiameters = + (diameterRange[1] - diameterRange[0]) / diameterRange[2] + 1; + + if (numberOfDiameters <= 0) { + ERROR("Error: It was impossible to create diameters array. "); + } + + diameter.resize(numberOfDiameters); + + for (int i = 0; i < numberOfDiameters; i++) { + diameter[i] = diameterRange[1] - i * diameterRange[2]; + } + + r_end = size; + dimx = nx; + dimy = ny; + dimz = nz; + + // Add additional layers for input/output reservoirs. + for (int i = 0; i < 3; i++) { + if (((flowAxis[i]) && (!allFaces)) || ((size[i] > 1) && (allFaces))) { + size[i] += 2; + r_ini[i] = 1; + r_end[i] = size[i] - 1; + } + } + + nx = size[0]; + ny = size[1]; + nz = size[2]; + + originalState.resize(nx, ny, nz); + currentState.resize(originalState.size()); + originalEDT.resize(originalState.size()); + trapped.resize(originalState.size()); + + currentState.fill(INJECTED); + trapped.fill(false); + + // For directional injection, set the outlet face to DISPLACED fluid. + for (int i = 0; i < 3; i++) { + int rMin[3] = {0, 0, 0}; + int rMax[3] = {nx, ny, nz}; + + if (!allFaces) { + if (flowAxis[i]) { + inletPos = flowPos ? 0 : size[i] - 1; + outletPos = flowPos ? size[i] - 1 : 0; + rMin[i] = outletPos; + rMax[i] = rMin[i] + 1; + setRegion(currentState, DISPLACED, rMin[0], rMax[0], rMin[1], + rMax[1], rMin[2], rMax[2]); + } + rChamberI[i] = flowAxis[i] ? inletPos : rMax[i] / 2; + rChamberO[i] = flowAxis[i] ? outletPos : rMax[i] / 2; + } + + else + rChamberI[i] = 0; + } + + int mapValue[255] = {-1}; + for (size_t idx = 0; idx < ReadValues.size(); idx++) { + + if ((ReadValues[idx] < 0) || (ReadValues[idx] > 255)) { + ERROR("Only values between 0 - 255 can be used as labels in " + "ReadValues"); + cout << ReadValues[idx] << endl; + } + if ((WriteValues[idx] < 0) || (WriteValues[idx] > 2)) { + ERROR("Only values between 0 (SOLID), 1 and 2 (INJECT/DISPLACED " + "FLUIDS) can be used as labels in WriteValues"); + } + mapValue[ReadValues[idx]] = (int)WriteValues[idx]; + } + + FILE *rawFile = fopen(mmfile.c_str(), "r"); + if (rawFile == NULL) { + ERROR("Error openning the file " + mmfile); + } + + long SEEK_BEGIN = ftell(rawFile); + long expectedSize = (long)(r_end[2] - r_ini[2]) * + (long)(r_end[1] - r_ini[1]) * + (long)(r_end[0] - r_ini[0]); + + fseek(rawFile, 0, SEEK_END); + + if (ftell(rawFile) != expectedSize) { + ERROR("File '" + mmfile + "' size is different from the expected (" + + to_string(expectedSize) + " bytes)."); + } + + fseek(rawFile, 0, + SEEK_BEGIN); // Move to the beginning of the file before reading. + + unsigned char readValue; + + NP = 0; + for (int z = r_ini[2]; z < r_end[2]; z++) { + for (int y = r_ini[1]; y < r_end[1]; y++) { + for (int x = r_ini[0]; x < r_end[0]; x++) { + + fread(&readValue, sizeof(unsigned char), 1, rawFile); + if (mapValue[readValue] == -1) { + ERROR(std::string("Not specified value in '" + filename + + "' at (" + to_string(x) + ", " + + to_string(y) + ", " + to_string(z) + + ").")); + } + + currentState(x, y, z) = (unsigned char)mapValue[readValue]; + if (currentState(x, y, z) != SOLID) + NP++; + } + } + } + + fclose(rawFile); + originalState = currentState; + + // Calculate the distance transform of the original solid phase. + edt_3d(SOLID, originalState, originalEDT); + + // Create the output CSV file and header when needed. + bool WriteHeader = false; + FILE *log_file = fopen("injection_output.csv", "r"); + if (log_file != NULL) + fclose(log_file); + else + WriteHeader = true; + + if (WriteHeader) { + log_file = fopen("injection_output.csv", "a+"); + fprintf(log_file, "step diameter_px diameter_um num_px_in frac_in " + "num_px_out frac_out\n"); + fclose(log_file); + } +} + +int Full_Morphology::calc(const int &step) { + + const int D = diameter[step]; + const double D24 = D * D / 4.0; + + IntArray auxMatrix(nx, ny, nz); + + // Mark the candidate invaded region as the region where the center of a + // D diameter sphere can be placed (erosion of the solid region) united with the reservoir. + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + auxMatrix(x, y, z) = (originalEDT(x, y, z) >= D24) || (originalState(x, y, z) == INJECTED) ? FOREGROUND : BACKGROUND; + } + } + } + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + + // Extract the label associated with the injection layer. + int chamber_label = auxMatrix(rChamberI[0], rChamberI[1], rChamberI[2]); + + // Filter only the voxels connected to the input where the center of a a sphere with the specified diameter + // can be placed. + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + auxMatrix(x, y, z) = (originalEDT(x, y, z) >= D24 && auxMatrix(x,y,z) == chamber_label) ? FOREGROUND : BACKGROUND; + } + } + } + + // Compute the distance from the filtered region + edt_3d(FOREGROUND, auxMatrix, auxMatrix); + + // Perform the dilation of the filtered region and performs the union of this result + // with the previous step stored at the currentState + for (int z = 0; z < nz; z++) { + for (int y = 0; y < ny; y++) { + for (int x = 0; x < nx; x++) { + if (auxMatrix(x, y, z) < D24) { + currentState(x, y, z) = INJECTED; + } + auxMatrix(x, y, z) = (currentState(x, y, z) == INJECTED) + ? FOREGROUND + : BACKGROUND; + } + } + } + + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + chamber_label = auxMatrix(rChamberI[0], rChamberI[1], rChamberI[2]); + + // Readds the DISPLACED fluid output layer + if (!allFaces) { + + if (flowAxis[0]) { + setRegion(currentState, DISPLACED, outletPos, outletPos + 1, 0, ny, + 0, nz); + } else if (flowAxis[1]) { + setRegion(currentState, DISPLACED, 0, nx, outletPos, outletPos + 1, + 0, nz); + } else if (flowAxis[2]) { + setRegion(currentState, DISPLACED, 0, nx, 0, ny, outletPos, + outletPos + 1); + } + } + + // For incompressible flow, disconnected displaced regions remain trapped + // in the final state. + if (!compressible) { + + for (int x = 0; x < nx; x++) { + for (int y = 0; y < ny; y++) { + for (int z = 0; z < nz; z++) { + + if (trapped(x, y, z)) + currentState(x, y, z) = DISPLACED; + + auxMatrix(x, y, z) = (currentState(x, y, z) == DISPLACED) + ? FOREGROUND + : BACKGROUND; + } + } + } + + component_labeling(auxMatrix, FOREGROUND, BACKGROUND); + chamber_label = auxMatrix(rChamberO[0], rChamberO[1], rChamberO[2]); + + for (int x = 0; x < nx; x++) { + for (int y = 0; y < ny; y++) { + for (int z = 0; z < nz; z++) { + if (currentState(x, y, z) == DISPLACED && + auxMatrix(x, y, z) != chamber_label) + trapped(x, y, z) = true; + } + } + } + } + + int injectedVolume = 0, displacedVolume = 0; + for (int z = r_ini[2]; z < r_end[2]; z++) { + for (int y = r_ini[1]; y < r_end[1]; y++) { + for (int x = r_ini[0]; x < r_end[0]; x++) { + + if (currentState(x, y, z) == INJECTED) + injectedVolume++; + else if (currentState(x, y, z) == DISPLACED) + displacedVolume++; + + int16_t *value = + &finalMap(x - r_ini[0], y - r_ini[1], z - r_ini[2]); + if (*value == -1 && currentState(x, y, z) == INJECTED) + *value = (int16_t)D; + if (*value == -1 && originalState(x, y, z) == SOLID) + *value = 0; + } + } + } + + if (saveImg && (D == diameter.back())) { + FILE *FRAW; + + FRAW = fopen("invasion_diameters.raw", "wb"); + fwrite(finalMap.data(), sizeof(int16_t), finalMap.length(), FRAW); + fclose(FRAW); + + FILE *FMHD = fopen("invasion_diameters.mhd", "w"); + + fprintf(FMHD, "ObjectType = Image\n"); + fprintf(FMHD, "NDims = 3\n"); + fprintf(FMHD, "DimSize = %d %d %d\n", dimx, dimy, dimz); + fprintf(FMHD, "ElementType = MET_SHORT\n"); + fprintf(FMHD, "ElementSpacing = %.1f %.1f %.1f\n", resolution, + resolution, resolution); + fprintf(FMHD, "ElementByteOrderMSB = False\n"); + fprintf(FMHD, "ElementDataFile = %s\n", "invasion_diameters.raw"); + fprintf(FMHD, "HeaderSize = 0\n"); + fclose(FMHD); + } + + FILE *log_file = fopen("injection_output.csv", "a"); + fprintf(log_file, "%d %d %f %d %f %d %f\n", step, D, D * resolution, + injectedVolume, injectedVolume / (1.0 * NP), displacedVolume, + displacedVolume / (1.0 * NP)); + fclose(log_file); + + return D; +} + +int main(int argc, char *argv[]) { + + if (argc != 2) + ERROR("Wrong number of parameters."); + + Full_Morphology fm(argc, argv); + + const int nsteps = fm.diameter.size(); + for (int step = 0; step < nsteps; step++) { + + int d = fm.calc(step); + cout << "Step " << step << ", D = " << d << " px." << endl; + } + + return 0; +} From bf4261234a59b8ab35ade539aae3fc70416b1a4d Mon Sep 17 00:00:00 2001 From: Diogo Nardelli Siebert Date: Mon, 4 May 2026 23:09:53 -0300 Subject: [PATCH 29/29] Fix an error in the parallel computation of the maximum pore size, where sumReduce was incorrectly used instead of maxReduce. --- analysis/morphology.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 9667740c..6a01812c 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -683,7 +683,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, // total Global is the number of nodes in the pore-space totalGlobal = Dm->Comm.sumReduce(count); - maxdistGlobal = Dm->Comm.sumReduce(maxdist); + maxdistGlobal = Dm->Comm.maxReduce(maxdist); double volume = double(nprocx * nprocy * nprocz) * double(nx - 2) * double(ny - 2) * double(nz - 2); double volume_fraction = totalGlobal / volume;