From 370a3c7669a4d175884aed13b096bf842e145e04 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshk@Dmytros-MacBook-Air.local>
Date: Fri, 4 Aug 2023 16:12:32 +1000
Subject: [PATCH 01/22] Add fixes for building on macOS

---
 src/makefile.main.Rt | 2 +-
 src/models.R         | 2 +-
 tools/dep.R          | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)
 mode change 100644 => 100755 tools/dep.R

diff --git a/src/makefile.main.Rt b/src/makefile.main.Rt
index b0af2d3e4..aa7297e52 100644
--- a/src/makefile.main.Rt
+++ b/src/makefile.main.Rt
@@ -293,7 +293,7 @@ wiki/schema/<?%s m ?>.xsd:$(SRC)/schema.xsd.Rt $(SRC)/conf.R <?%s d ?>/<?%s m ?>
 %_sp.c:%.c
 	@echo "  SP-CONST   $@"
 	@$(MKPATH) $@
-	@cat $< | sed -r 's/([^a-zA-Z0-9][0-9]+\.[0-9]*([eE][-+]?[0-9]+)?)[fL]+/\1f/g' > $@
+	@cat $< | sed -E 's/([^a-zA-Z0-9][0-9]+\.[0-9]*([eE][-+]?[0-9]+)?)[fL]+/\1f/g' > $@
 
 
 
diff --git a/src/models.R b/src/models.R
index a9a1add46..b4db6571f 100644
--- a/src/models.R
+++ b/src/models.R
@@ -19,7 +19,7 @@ get.model.names = function(path) {
 
 get.models = function() {
 	M1 = get.model.dirs("git ls-files | grep 'conf.mk$'")
-	M2 = get.model.dirs("find models/ -name 'conf.mk'")
+	M2 = get.model.dirs("find models -name 'conf.mk'")
 	M3 = union(M1,M2)
 	Models = do.call(rbind, lapply(M3,function (m) {
 		ret = get.model.names(m)
diff --git a/tools/dep.R b/tools/dep.R
old mode 100644
new mode 100755
index 7492616b8..9fc102c0e
--- a/tools/dep.R
+++ b/tools/dep.R
@@ -9,7 +9,7 @@ opt <- parse_args(OptionParser(usage="Usage: ADmod -f inputfile [-o outputfile]"
 if (opt$path != "") setwd(opt$path)
 
 # 'grep' all the includes
-f = pipe("grep '# *include' `find -regex '.*[.]\\(c\\|cu\\|cpp\\|h\\|hpp\\)'` | sed -n 's/^\\([^:]*\\):[ \\t]*#[ \\t]*include[ \\t]*[\"<]\\(.*\\)[>\"]/\\1,\\2/gp'")
+f = pipe("grep '# *include' `find . -type f \\( -name '*.c' -o -name '*.cu' -o -name '*.cpp' -o -name '*.h' -o -name '*.hpp' \\)` | sed -n -e 's/^\\([^:]*\\):[ \\t]*#include[ \\t]*[\"<]\\(.*\\)[>\"]/\\1,\\2/gp'")
 w = read.csv(f,col.names=c("file","dep"), stringsAsFactor=F, header=FALSE);
 
 # make the relative paths in include statement relative to *here*

From c3ac8f416a2760449612007773d61284425cfa7a Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Thu, 8 Feb 2024 14:48:15 +1000
Subject: [PATCH 02/22] Adjust boundary computation for angle 90 (or other
 similar cases)

Compute pf_f in the normal direction as soon as possible, because
we probably want to use it later in the computations
---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 55 +++++++++----------
 1 file changed, 26 insertions(+), 29 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 024f94c20..0ceb6b12b 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -689,37 +689,17 @@ CudaDeviceFunction void calcWallPhase(){
 	PhaseF = PhaseF(0,0,0); //For fluid nodes.
 	if ( IamWall || IamSolid ) {
 		real_t a, h, pf_f;
-
-		// This is needed, because otherwise geometric_staircaseimp performance will drop
+        // First compute pf_f without staircase improvement
+		//
+		// Weird selection here is needed, because otherwise geometric_staircaseimp performance will drop
 		// (presumably because of the dynamic access that it has)
         <?R phase_field_field_name <- if (Options$geometric) 'gradPhi_PhaseF' else 'PhaseF' ?>
-
         pf_f = <?%s phase_field_field_name ?>_dyn(nw_x, nw_y, nw_z);
-
 		h = 0.5 * sqrt(nw_x*nw_x + nw_y*nw_y + nw_z*nw_z);
 
-        // handling special cases
-		if (h < 0.001) {
-			// If I am a wall/solid node and I am surrounded by solid nodes
-			PhaseF = 1;
-		} else if (fabs(radAngle - PI/2.0) < 1e-4) {
-			// If I am not surrounded, but contact angle is pi/2 (90d)
-			PhaseF = pf_f;
-		} else if (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?>) {
-			// Pass for now. Dont calculate anything, just set it to some huge number, purely to make sure
-			// that the value is corrected in calcWall_correction, and otherwise the error would be visible
-			// and hopefulyl break the simulation
-			PhaseF = <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>;
-		} else if  (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_FURTHER_NEXT_NODE ?>) {
-			// Eventhough I am geometric boundary condition, still apply surface energy
-			// here because otherwise we cant really apply anything else
-			a = -h * (4.0/IntWidth) * cos( radAngle );
-			PhaseF = (1 + a - sqrt( (1+a)*(1+a) - 4*a*pf_f))/(a+1e-12) - pf_f;
-		} else {
-			// normal calculation with picking correct form depending on the boundary condition
-
+        // if enabled compute staircase improved version of pf_f,
+        // and wallback to pf_f
 #ifdef OPTIONS_staircaseimp
-
 			int face_index = int(triangle_index) / 8;
 			int face_triangle_index = int(triangle_index) % 8;
 			int vertex_coords[3] = {0,0,0};
@@ -729,9 +709,6 @@ CudaDeviceFunction void calcWallPhase(){
 				int v<?%s v ?>_x = vertex_coords[0], v<?%s v ?>_y = vertex_coords[1], v<?%s v ?>_z = vertex_coords[2];
 			<?R } ?>
 
-
-			<?R phase_field_field_name <- if (Options$geometric) 'gradPhi_PhaseF' else 'PhaseF' ?>
-
 			<?R for (v in c(1, 2, 3)) { ?>
 				real_t pf_v<?%s v?> = <?%s phase_field_field_name ?>_dyn(v<?%s v?>_x, v<?%s v?>_y, v<?%s v?>_z);
 			<?R } ?>
@@ -744,7 +721,27 @@ CudaDeviceFunction void calcWallPhase(){
 				pf_f = pf_interpolated;
 			}
 #endif
-			
+
+        // handling special cases
+		if (h < 0.001) {
+			// If I am a wall/solid node and I am surrounded by solid nodes
+			PhaseF = 1;
+		} else if (fabs(radAngle - PI/2.0) < 1e-4) {
+			// If I am not surrounded, but contact angle is pi/2 (90d)
+			PhaseF = pf_f;
+		} else if (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?>) {
+			// Pass for now. Dont calculate anything, just set it to some huge number, purely to make sure
+			// that the value is corrected in calcWall_correction, and otherwise the error would be visible
+			// and hopefulyl break the simulation
+			PhaseF = <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>;
+		} else if  (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_FURTHER_NEXT_NODE ?>) {
+			// Eventhough I am geometric boundary condition, still apply surface energy
+			// here because otherwise we cant really apply anything else
+			a = -h * (4.0/IntWidth) * cos( radAngle );
+			PhaseF = (1 + a - sqrt( (1+a)*(1+a) - 4*a*pf_f))/(a+1e-12) - pf_f;
+		} else {
+			// normal calculation with picking correct form depending on the boundary condition
+
 #ifndef OPTIONS_geometric
 			// Case 1: Apply surface energy BC (with calculated pf_f with standard or staircase improvement)
 			a = -h * (4.0/IntWidth) * cos( radAngle );

From c3b563f01b00cc6de6d645c41b9220909ed38af7 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 13 Feb 2024 10:43:50 +1000
Subject: [PATCH 03/22] Split normal initialisation in two part

This way we can both initialise external
and internally based on some approximations
---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 125 +++++++++++++-----
 .../multiphase/d3q27_pf_velocity/Dynamics.R   |  26 +++-
 .../d3q27_pf_velocity/Dynamics.c.Rt           |  16 +++
 3 files changed, 125 insertions(+), 42 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 0ceb6b12b..1db215670 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -673,11 +673,32 @@ CudaDeviceFunction void calcPhaseGradCloseToBoundary(){
  */
 CudaDeviceFunction void calcWallPhase_correction() {
 	PhaseF = PhaseF(0,0,0);
-
 	if (IsSpecialBoundaryPoint == <?%s NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?> ) {
 		// take the phase field calculated already from the node in front.
 		// Might as well calculate using the neighbors, e.g. averaging
-		PhaseF = PhaseF_dyn(nw_x, nw_y, nw_z);
+		real_t point_in_the_normal_direction = PhaseF_dyn(nw_x, nw_y, nw_z);
+        // in some cases we might (TODO: investigate) still point into the solid node
+        // that has normal pointing into another, we need to avoid this
+        if (point_in_the_normal_direction == <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>) {
+            // ignore first node
+            int count = 0;
+            real_t fluid_average = 0.0;
+            real_t x, y, z;
+
+			for (int i=1; i<27 ; i++){
+				x = d3q27_ex[i];
+				y = d3q27_ey[i];
+				z = d3q27_ez[i];
+                if (!IsBoundary_dyn(int(x), int(y), int(z)) && PhaseF_dyn(int(x), int(y), int(z)) > -0.5) {
+                    fluid_average += PhaseF_dyn(int(x), int(y), int(z));
+                    count++;
+                }
+			}
+            // print sum and average and count
+            PhaseF = count == 0 ? 1 : fluid_average / count;
+        } else {
+            PhaseF = point_in_the_normal_direction;
+        }
 	}
 }
 
@@ -848,6 +869,57 @@ CudaDeviceFunction real_t getSpecialBoundaryPoint() {
 	return IsSpecialBoundaryPoint;
 }
 
+/* Calculate the actual wall norm based on the solid
+ *field */
+CudaDeviceFunction real_t Init_real_wallNorm()
+{
+    // set initial values
+    nw_x = 0.0;
+    nw_y = 0.0;
+    nw_z = 0.0;
+    if ( IamWall || IamSolid ) {
+
+		int i,j,k;
+	  	real_t tmp = 0.0;
+	  	for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
+			tmp += PhaseF_dyn(i,j,k);
+	  	}}}
+
+	  	if (abs(tmp) > 26000) {
+			// I am surrounded by all solid nodes (sum(pf) = 27*-999 = -26973 if surrounded):
+			nw_x = 0.0; nw_y = 0.0; nw_z = 0.0;
+	  	} else {
+			// no I am not surrounded, so calc normal:
+			int solidFlag[27];
+			// Calculate the normal direction based converting
+			// negative PhaseF into actual solid flags
+			<?R
+			    myN   = PV(paste0("myNorm[",1:3-1,"]"))
+			    pf    = PV(paste0("PhaseF(",U[,1],",",U[,2],",",U[,3],")/-998"))
+			    solid = PV(paste0("solidFlag[",1:27-1,"]"))
+
+			    C(solid, pf)
+			?>
+			for (i=0;i<27;i++){
+				nw_x += wg[i] * solidFlag[i] * d3q27_ex[i];
+				nw_y += wg[i] * solidFlag[i] * d3q27_ey[i];
+				nw_z += wg[i] * solidFlag[i] * d3q27_ez[i];
+			}
+			nw_x *= -1.0/3.0; nw_y *= -1.0/3.0; nw_z *= -1.0/3.0;
+			tmp = nw_x*nw_x + nw_y*nw_y + nw_z*nw_z;
+
+            nw_x = nw_x / sqrt(tmp);
+            nw_y = nw_y / sqrt(tmp);
+            nw_z = nw_z / sqrt(tmp);
+
+        }
+    }
+    else {
+        // I am a fluid node, I dont need no solid normal,
+        // should be set to zero before. Branch kept
+        // for clarity
+    }
+}
 
 /*
  * Initialise and set: wall normals, get interpolating triangles and their coefficients,
@@ -863,7 +935,8 @@ CudaDeviceFunction void Init_wallNorm(){
 	coeff_v2 = 0;
 	coeff_v3 = 0;
 #endif
-	if ( IamWall || IamSolid ) { 
+    //printf("PhaseF value is %f\n", PhaseF);
+	if ( IamWall || IamSolid ) {
         IsBoundary = 1.0;
 		int i,j,k;
 	  	real_t tmp = 0.0;
@@ -880,42 +953,29 @@ CudaDeviceFunction void Init_wallNorm(){
 			nw_actual_z = 0.0;
 #endif
 	  	} else {
-			// no I am not surrounded, so calc normal:
-			int solidFlag[27];
+			// Calculate the closest discrete direction for normal:
+			real_t tmp = 0.0;
 			int maxi = 0;
-			real_t myNorm[3] = {0.0,0.0,0.0};
 			real_t maxn=0.0, dot;
 
-			// Calculate the normal direction based converting
-			// negative PhaseF into actual solid flags
-			<?R
-			    myN   = PV(paste0("myNorm[",1:3-1,"]"))
-			    pf    = PV(paste0("PhaseF(",U[,1],",",U[,2],",",U[,3],")/-998"))
-			    solid = PV(paste0("solidFlag[",1:27-1,"]"))
-
-			    C(solid, pf)
-			?>
-			for (i=0;i<27;i++){
-				myNorm[0] += wg[i] * solidFlag[i] * d3q27_ex[i];
-				myNorm[1] += wg[i] * solidFlag[i] * d3q27_ey[i];
-				myNorm[2] += wg[i] * solidFlag[i] * d3q27_ez[i];
-
-			}
-			myNorm[0] *= -1.0/3.0;myNorm[1] *= -1.0/3.0;myNorm[2] *= -1.0/3.0;
-			tmp = myNorm[0]*myNorm[0] + myNorm[1]*myNorm[1] + myNorm[2]*myNorm[2];
-
-			// Calculate the closest discrete direction for normal:
 			for (i = 0; i<27; i++) {
-				dot = (myNorm[0]*d3q27_ex[i] + myNorm[1]*d3q27_ey[i] + myNorm[2]*d3q27_ez[i]) /
-			      		sqrt( tmp*(d3q27_ex[i]*d3q27_ex[i] + d3q27_ey[i]*d3q27_ey[i] +
-					 	   d3q27_ez[i]*d3q27_ez[i]) + 1e-12);
+			tmp = nw_x*nw_x + nw_y*nw_y + nw_z*nw_z;
+            dot = (nw_x*d3q27_ex[i] + nw_y*d3q27_ey[i] + nw_z*d3q27_ez[i]) /
+                    sqrt( tmp*(d3q27_ex[i]*d3q27_ex[i] + d3q27_ey[i]*d3q27_ey[i] +
+                        d3q27_ez[i]*d3q27_ez[i]) + 1e-12);
 				if (dot > maxn) {
 					maxn = dot; maxi = i;
 				}
 			}
+#ifdef OPTIONS_staircaseimp
+            // save original values before they get truncated
+			real_t original_normal[3] = {nw_x, nw_y, nw_z};
+#endif
 			if (maxi < 0) {
 				// This should not happen ?
-				nw_x = 0.0;nw_y = 0.0;nw_z = 0.0;
+				nw_x = 0.0;
+                nw_y = 0.0;
+                nw_z = 0.0;
 			} else {
 				nw_x = d3q27_ex[maxi];
 				nw_y = d3q27_ey[maxi];
@@ -930,7 +990,6 @@ CudaDeviceFunction void Init_wallNorm(){
 			}
 
 #ifdef OPTIONS_staircaseimp
-			real_t truncated_normals[3] = {nw_x, nw_y, nw_z};
 			real_t intersection_point[3];
 			int face_index;
 			real_t coeff[3];
@@ -938,7 +997,7 @@ CudaDeviceFunction void Init_wallNorm(){
 			int t_index;
 			int t_index2;
 
-			calcIntersectionTriangleId(myNorm, t_index, face_index, coeff, intersection_point, t_index2, coeff2);
+			calcIntersectionTriangleId(original_normal, t_index, face_index, coeff, intersection_point, t_index2, coeff2);
 
 			// make sure normal vector is extended to the surface (to use it during interpolation)
 			nw_actual_x = intersection_point[0];
@@ -1027,10 +1086,6 @@ CudaDeviceFunction void Init_wallNorm(){
 
 	  	}
 	} else {
-	// I am a fluid node, I dont need no solid normal.
-		nw_x = 0.0;
-		nw_y = 0.0;
-		nw_z = 0.0;
 #ifdef OPTIONS_staircaseimp
 		nw_actual_x = 0.0;
 		nw_actual_y = 0.0;
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 2023981b4..5118c085e 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -14,6 +14,11 @@ AddDensity(name="Init_UY_External", group="init", comment="free stream velocity"
 AddDensity(name="Init_UZ_External", group="init", comment="free stream velocity", parameter=TRUE)
 AddDensity(name="Init_PhaseField_External", group="init", dx=0,dy=0,dz=0, parameter=TRUE)
 
+# for initialising the normals
+AddDensity(name="Init_nwx_external", group="init_normals", dx=0,dy=0,dz=0, parameter=TRUE)
+AddDensity(name="Init_nwy_external", group="init_normals", dx=0,dy=0,dz=0, parameter=TRUE)
+AddDensity(name="Init_nwz_external", group="init_normals", dx=0,dy=0,dz=0, parameter=TRUE)
+
 # macroscopic params
 # - consider migrating to fields
 AddDensity(name="pnorm", dx=0, dy=0, dz=0, group="Vel")
@@ -125,16 +130,21 @@ if (Options$thermo){
 	AddStage("calcPhase", "calcPhaseF", save=Fields$name=="PhaseF", load=DensityAll$group %in% load_phase)
 	AddStage("BaseIter" , "Run", save=Fields$group %in% save_iteration, load=DensityAll$group %in% load_iteration )
 	AddStage(name="InitFromFieldsStage", load=DensityAll$group %in% "init",read=FALSE, save=Fields$group %in% save_initial_PF)
+    AddStage(name="InitFromFieldsStageWithNormals", load=DensityAll$group %in% c("init_normals", "init"),read=FALSE, save=Fields$group %in% c("nw", save_initial_PF))
+
 	# STAGES FOR VARIOUS OPTIONS
 	if (Options$geometric){
-		AddStage("WallInit_CA"  , "Init_wallNorm", save=Fields$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
+
+		AddStage("WallInit_Real"  , "Init_real_wallNorm", save=Fields$group %in% c("nw"))
+		AddStage("WallInit_CA"  , "Init_wallNorm", load=DensityAll$group %in% c("nw"), save=Fields$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
 		AddStage("calcWall_CA"  , "calcWallPhase", save=Fields$name %in% c("PhaseF"), load=DensityAll$group %in% c("nw", "gradPhi", "PF", "solid_boundary", extra_fields_to_load_for_bc))
 
 		AddStage('calcPhaseGrad', "calcPhaseGrad", load=DensityAll$group %in% c("nw", "PF", "solid_boundary"), save=Fields$group=="gradPhi")
 		AddStage('calcPhaseGrad_init', "calcPhaseGrad_init", load=DensityAll$group %in% c("nw", "PF", "solid_boundary"), save=Fields$group=="gradPhi")
 		AddStage("calcWallPhase_correction", "calcWallPhase_correction", save=Fields$name=="PhaseF", load=DensityAll$group %in% c("nw", "solid_boundary"))
 	} else {
-		AddStage("WallInit" , "Init_wallNorm", save=Fields$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
+		AddStage("WallInit_Real"  , "Init_real_wallNorm", save=Fields$group %in% c("nw"))
+		AddStage("WallInit" , "Init_wallNorm", load=DensityAll$group %in% c("nw"), save=Fields$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
 		AddStage("calcWall" , "calcWallPhase", save=Fields$name=="PhaseF", load=DensityAll$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
 		AddStage("calcWallPhase_correction", "calcWallPhase_correction", save=Fields$name=="PhaseF", load=DensityAll$group %in% c("nw", "solid_boundary"))
 	}
@@ -156,16 +166,18 @@ if (Options$thermo){
 		AddAction("TempToSteadyState", c("CopyDistributions","RK_1", "RK_2", "RK_3", "RK_4","NonLocalTemp"))
 		AddAction("Iteration", c("BaseIter", "calcPhase", "calcWall","RK_1", "RK_2", "RK_3", "RK_4","NonLocalTemp"))
 		AddAction("IterationConstantTemp", c("BaseIter", "calcPhase", "calcWall","CopyThermal"))
-		AddAction("Init"     , c("PhaseInit","WallInit" , "calcWall","BaseInit"))
+		AddAction("Init"     , c("PhaseInit","WallInit_Real","WallInit" , "calcWall","BaseInit"))
 	} else if (Options$geometric) {
         calcGrad <- if (Options$isograd)  "calcPhaseGrad" else "calcPhaseGrad_init"
         AddAction("Iteration", c("BaseIter", "calcPhase",  calcGrad, "calcWall_CA", "calcWallPhase_correction"))
-	    AddAction("Init"     , c("PhaseInit","WallInit_CA" , "calcPhaseGrad_init"  , "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
-	    AddAction("InitFields"     , c("InitFromFieldsStage","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
+	    AddAction("Init"     , c("PhaseInit","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init"  , "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
+	    AddAction("InitFields"     , c("InitFromFieldsStage","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
+	    AddAction("InitFieldsWithNormals"     , c("InitFromFieldsStageWithNormals","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
     } else {
 		AddAction("Iteration", c("BaseIter", "calcPhase", "calcWall", "calcWallPhase_correction"))
-		AddAction("Init"     , c("PhaseInit","WallInit" , "calcWall","calcWallPhase_correction", "BaseInit"))
-		AddAction("InitFields", c("InitFromFieldsStage","WallInit" , "calcWall", "calcWallPhase_correction", "BaseInit"))
+		AddAction("Init"     , c("PhaseInit", "WallInit_Real", "WallInit" , "calcWall","calcWallPhase_correction", "BaseInit"))
+		AddAction("InitFields", c("InitFromFieldsStage","WallInit_Real","WallInit" , "calcWall", "calcWallPhase_correction", "BaseInit"))
+		AddAction("InitFieldsWithNormals", c("InitFromFieldsStageWithNormals", "WallInit" , "calcWall", "calcWallPhase_correction", "BaseInit"))
 	}
 #######################
 ########OUTPUTS########
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
index 02cf4ce43..ea1059442 100755
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
@@ -287,6 +287,22 @@ CudaDeviceFunction void InitFromFieldsStage()
 	pnorm = 0.0; // initialise as zero and fill in later stage
 }
 
+
+CudaDeviceFunction void InitFromFieldsStageWithNormals()
+{
+	PhaseF = Init_PhaseField_External;
+	U = Init_UX_External;
+	V = Init_UY_External;
+	W = Init_UZ_External;
+	if ( IamWall || IamSolid ) PhaseF = -999;
+	pnorm = 0.0; // initialise as zero and fill in later stage
+    nw_x = Init_nwx_external;
+    nw_y = Init_nwy_external;
+    nw_z = Init_nwz_external;
+}
+
+
+
 CudaDeviceFunction void specialCases_Init()
 {
 	#ifdef OPTIONS_thermo

From 5bcc8d4cc6def5d71afc3b65ca56e4d29774f34a Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Wed, 28 Feb 2024 16:45:41 +1000
Subject: [PATCH 04/22] Fix geometric model

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 1db215670..1fd7c74cc 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -735,8 +735,9 @@ CudaDeviceFunction void calcWallPhase(){
 			<?R } ?>
 
 
-			// dont do staircase improvement if any of the interpolating nodes are solid
-			if (IsSpecialBoundaryPoint != <?%d NEXT_INTERPOLATING_POINTS_ARE_SOLID ?>) {
+			// dont do staircase improvement if any of the interpolating nodes are solid, make sure to cover
+			// both cases......
+			if (IsSpecialBoundaryPoint != <?%d NEXT_INTERPOLATING_POINTS_ARE_SOLID ?> && IsSpecialBoundaryPoint != <?%d BOTH_NEXT_AND_NNEXT_INTERPOLATING_POINTS_ARE_SOLID ?>) {
 				real_t pf_interpolated = coeff_v1 * pf_v1 + coeff_v2 * pf_v2 + coeff_v3 * pf_v3;
 				h = 0.5 * sqrt(nw_actual_x*nw_actual_x + nw_actual_y*nw_actual_y + nw_actual_z*nw_actual_z);
 				pf_f = pf_interpolated;

From fd057e624f2224e65c664fbaf1a9e9365214f1d9 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Wed, 13 Mar 2024 10:15:29 +1000
Subject: [PATCH 05/22] Fix undified behaviours

---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 35 +++++++++++--------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 1fd7c74cc..f7cec0549 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -718,8 +718,21 @@ CudaDeviceFunction void calcWallPhase(){
         pf_f = <?%s phase_field_field_name ?>_dyn(nw_x, nw_y, nw_z);
 		h = 0.5 * sqrt(nw_x*nw_x + nw_y*nw_y + nw_z*nw_z);
 
-        // if enabled compute staircase improved version of pf_f,
-        // and wallback to pf_f
+        // /now handle special cacses, that don't even
+        // need to calculate the phase field now (or have
+        // access to different fields)
+		if (h < 0.001) {
+			// If I am a wall/solid node and I am surrounded by solid nodes
+			PhaseF = 1;
+        } else if (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?>) {
+			// Pass for now. Dont calculate anything, just set it to some huge number, purely to make sure
+			// that the value is corrected in calcWall_correction, and otherwise the error would be visible
+			// and hopefulyl break the simulation
+			PhaseF = <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>;
+        } else {
+
+            // if enabled compute staircase improved version of pf_f,
+            // and wallback to pf_f, which is required is many places
 #ifdef OPTIONS_staircaseimp
 			int face_index = int(triangle_index) / 8;
 			int face_triangle_index = int(triangle_index) % 8;
@@ -743,22 +756,15 @@ CudaDeviceFunction void calcWallPhase(){
 				pf_f = pf_interpolated;
 			}
 #endif
-
         // handling special cases
-		if (h < 0.001) {
-			// If I am a wall/solid node and I am surrounded by solid nodes
-			PhaseF = 1;
-		} else if (fabs(radAngle - PI/2.0) < 1e-4) {
-			// If I am not surrounded, but contact angle is pi/2 (90d)
+		if (fabs(radAngle - PI/2.0) < 1e-4) {
+			// If I am not surrounded, but contact angle is pi/2 (90d). Staircase improved
+			// already handled before
 			PhaseF = pf_f;
-		} else if (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?>) {
-			// Pass for now. Dont calculate anything, just set it to some huge number, purely to make sure
-			// that the value is corrected in calcWall_correction, and otherwise the error would be visible
-			// and hopefulyl break the simulation
-			PhaseF = <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>;
 		} else if  (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_FURTHER_NEXT_NODE ?>) {
 			// Eventhough I am geometric boundary condition, still apply surface energy
-			// here because otherwise we cant really apply anything else
+			// here because otherwise we cant really apply anything else. Staircase improved
+			// already handled before
 			a = -h * (4.0/IntWidth) * cos( radAngle );
 			PhaseF = (1 + a - sqrt( (1+a)*(1+a) - 4*a*pf_f))/(a+1e-12) - pf_f;
 		} else {
@@ -864,6 +870,7 @@ CudaDeviceFunction void calcWallPhase(){
 #endif // end boundary condition pick
 		}
 	}
+  }
 }
 
 CudaDeviceFunction real_t getSpecialBoundaryPoint() {

From 1b4eca66b5e2f014ff0c407ebed18dfaf585951d Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Thu, 22 Jun 2023 09:40:09 +1000
Subject: [PATCH 06/22] Add utilities

---
 models/multiphase/d3q27_pf_velocity/Dynamics.R    | 2 ++
 models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 5118c085e..4d0be761d 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -320,3 +320,5 @@ if (Options$thermo){
 	AddGlobal(name="FluxX",comment='flux in x direction for flux_nodes', unit="1")
 	AddGlobal(name="FluxY",comment='flux in y direction for flux_nodes', unit="1")
 	AddGlobal(name="FluxZ",comment='flux in z direction for flux_nodes', unit="1")
+    AddGlobal(name="LiquidSaturation", comment="Liquid saturation(number of liquid nodes)", unit="1")
+    AddGlobal(name="GasSaturation", comment="Gas saturation(number of gas nodes)", unit="1")
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
index ea1059442..87bf0d626 100755
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
@@ -749,12 +749,14 @@ CudaDeviceFunction void updateMyGlobals(real_t pf){
 	    AddToGasTotalVelocityZ( tmpPF*W );
 	    AddToGasTotalPhase( tmpPF );
 	    AddToXLocation( tmpPF*X );
+        AddToGasSaturation(1);
 	} else {
 	    AddToLiqTotalVelocity( pf*sqrt(u2mag));
 	    AddToLiqTotalVelocityX( U*pf );
 	    AddToLiqTotalVelocityY( V*pf );
 	    AddToLiqTotalVelocityZ( W*pf );
 	    AddToLiqTotalPhase( pf );
+        AddToLiquidSaturation(1);
 	}
 
 	if ((NodeType & NODE_ADDITIONALS) == NODE_flux_nodes) {

From 86d7854c02ea419614b26e811eaede0f80075ef8 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Fri, 12 Apr 2024 10:04:28 +1000
Subject: [PATCH 07/22] Fix weighting

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index f7cec0549..9baab6104 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -684,13 +684,12 @@ CudaDeviceFunction void calcWallPhase_correction() {
             int count = 0;
             real_t fluid_average = 0.0;
             real_t x, y, z;
-
 			for (int i=1; i<27 ; i++){
 				x = d3q27_ex[i];
 				y = d3q27_ey[i];
 				z = d3q27_ez[i];
                 if (!IsBoundary_dyn(int(x), int(y), int(z)) && PhaseF_dyn(int(x), int(y), int(z)) > -0.5) {
-                    fluid_average += PhaseF_dyn(int(x), int(y), int(z));
+                    fluid_average += wg[i]*PhaseF_dyn(int(x), int(y), int(z));
                     count++;
                 }
 			}

From 731ac874d5ef354052e01c1b9eed9e2bc4ff01be Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 21 May 2024 10:57:57 +1000
Subject: [PATCH 08/22] WIP Pressure boundary condition fix

---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 39 ++++++++++++++++---
 .../d3q27_pf_velocity/Dynamics.c.Rt           | 23 ++++++++++-
 2 files changed, 55 insertions(+), 7 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 9baab6104..d2ed4f364 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -73,6 +73,19 @@ CudaDeviceFunction void <?R C(my_velocity_boundaries[ii])?>(){
 
 CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
     real_t d = getRho();
+
+    real_t myPhaseField = 0.0;
+    // TODO: Need a proper fix soon
+    if ((NodeType & NODE_BOUNDARY) == NODE_WPressure){
+        d = Density_l;
+        myPhaseField = 0.0;
+    }
+
+    if ((NodeType & NODE_BOUNDARY) == NODE_EPressure){
+        d = Density_h;
+        myPhaseField = 1.0;
+    }
+
 	real_t pstar = Pressure / (d*cs2);
 	<?R
         EQ = MRT_eq(U,PV(1),u,mat=t(M))
@@ -90,7 +103,7 @@ CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
 	?>
 	<?R
 		U_PF = U[1:PF_velocities,]
-       	pf = PV("PhaseField")
+       	pf = PV("myPhaseField")
 		n = my_normal_directions[ii,]
 		heq = pf * EQ_h$feq
 		bounce = Bounce(U_PF)
@@ -723,6 +736,7 @@ CudaDeviceFunction void calcWallPhase(){
 		if (h < 0.001) {
 			// If I am a wall/solid node and I am surrounded by solid nodes
 			PhaseF = 1;
+
         } else if (IsSpecialBoundaryPoint == <?%d NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?>) {
 			// Pass for now. Dont calculate anything, just set it to some huge number, purely to make sure
 			// that the value is corrected in calcWall_correction, and otherwise the error would be visible
@@ -888,8 +902,23 @@ CudaDeviceFunction real_t Init_real_wallNorm()
 
 		int i,j,k;
 	  	real_t tmp = 0.0;
+        real_t PhaseF_tmp[3][3][3];
+        // mirror in the beginning, in case of non-periodic boundaries
+        int nx = constContainer.nx;
+        for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
+            if ((i == -1) && (X == 0.5)) {
+                PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
+            }
+            else if ((i == 1) && (X == nx - 0.5)) {
+                PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
+            }
+            else {
+                PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(i,j,k);
+            }
+        }}
+        }
 	  	for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
-			tmp += PhaseF_dyn(i,j,k);
+			tmp += PhaseF_tmp[i+1][j+1][k + 1];
 	  	}}}
 
 	  	if (abs(tmp) > 26000) {
@@ -901,17 +930,18 @@ CudaDeviceFunction real_t Init_real_wallNorm()
 			// Calculate the normal direction based converting
 			// negative PhaseF into actual solid flags
 			<?R
-			    myN   = PV(paste0("myNorm[",1:3-1,"]"))
-			    pf    = PV(paste0("PhaseF(",U[,1],",",U[,2],",",U[,3],")/-998"))
+			    pf    = PV(paste0("PhaseF_tmp[",U[,1] + 1,"][",U[,2] + 1,"][",U[,3] + 1,"]/-998"))
 			    solid = PV(paste0("solidFlag[",1:27-1,"]"))
 
 			    C(solid, pf)
 			?>
 			for (i=0;i<27;i++){
+                // no need to calculate the normal direction if in the beginning of the boundary
 				nw_x += wg[i] * solidFlag[i] * d3q27_ex[i];
 				nw_y += wg[i] * solidFlag[i] * d3q27_ey[i];
 				nw_z += wg[i] * solidFlag[i] * d3q27_ez[i];
 			}
+
 			nw_x *= -1.0/3.0; nw_y *= -1.0/3.0; nw_z *= -1.0/3.0;
 			tmp = nw_x*nw_x + nw_y*nw_y + nw_z*nw_z;
 
@@ -942,7 +972,6 @@ CudaDeviceFunction void Init_wallNorm(){
 	coeff_v2 = 0;
 	coeff_v3 = 0;
 #endif
-    //printf("PhaseF value is %f\n", PhaseF);
 	if ( IamWall || IamSolid ) {
         IsBoundary = 1.0;
 		int i,j,k;
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
index 87bf0d626..f2b632c7d 100755
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
@@ -81,12 +81,16 @@ CudaDeviceFunction vector_t getU(){
 }
 
 CudaDeviceFunction real_t getPstar(){
+    updateBoundary();
 	return <?R C(sum(g)) ?>;
 }
 
 CudaDeviceFunction real_t getP(){
+    // note: vtk export should happen after the boundary is updated,
+    // otherwise there would be some strange values
+    updateBoundary();
 	real_t d = getRho();
-	real_t pstar = getPstar();
+	real_t pstar = <?R C(sum(g)) ?>;
 	return pstar*d*cs2;
 }
 
@@ -143,7 +147,14 @@ CudaDeviceFunction vector_t calcGradPhi()
 		    <?R
 			IsotropicGrad('gradPhi', 'PhaseF')
 		    ?>
-	#endif	
+
+        auto boundary_marker = NodeType & NODE_BOUNDARY;
+        if (boundary_marker == NODE_EPressure || boundary_marker == NODE_WPressure) {
+            gradPhi.x = 0.0;
+            gradPhi.y = 0.0;
+            gradPhi.z = 0.0;
+        }
+	#endif
 	return gradPhi;
 }
 
@@ -218,11 +229,19 @@ CudaDeviceFunction real_t calcMu(real_t C)
 			?>
 	#endif
 	#ifdef OPTIONS_thermo
+
 		mu = 4.0*(12.0*SurfaceTension(0,0,0)/IntWidth)*(C-PhaseField_l)*(C-PhaseField_h)*(C-pfavg)
 			- (1.5 *SurfaceTension(0,0,0)*IntWidth) * lpPhi;
 	#else
+
+        auto boundary_marker = NodeType & NODE_BOUNDARY;
+        if (boundary_marker == NODE_EPressure || boundary_marker == NODE_WPressure) {
+            lpPhi = 0.0;
+        }
+
 		mu = 4.0*(12.0*sigma/IntWidth)*(C-PhaseField_l)*(C-PhaseField_h)*(C-pfavg)
 			- (1.5 *sigma*IntWidth) * lpPhi;
+
 	#endif
 	return mu;
 }

From c45c99737a4845a44393b01ae1cec23a34beb6e1 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Fri, 24 May 2024 08:50:07 +1000
Subject: [PATCH 09/22] Add Lukaszes boundary conditions

---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 102 +++++++++++++-----
 1 file changed, 76 insertions(+), 26 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index d2ed4f364..e580d3fc6 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -87,34 +87,76 @@ CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
     }
 
 	real_t pstar = Pressure / (d*cs2);
-	<?R
-        EQ = MRT_eq(U,PV(1),u,mat=t(M))
-       	pstar = PV("pstar")
-		n = my_normal_directions[ii,]
-		geq = pstar*w_g + (EQ$feq - w_g)
-		bounce = Bounce(U)
-		sel = as.vector( (U %*% n) < 0)
-
-		sel2 = as.vector( ( U %*% n) == 0)	
-		exM = (g[sel2] - geq[sel2]) %*% rep(1,sum(sel2)) 
-		Unknowns = 1.0/sum(sel)
 
-		C(g[sel], geq[sel] + geq[bounce][sel] - g[bounce][sel] - 0.5 * Unknowns * (exM))
-	?>
-	<?R
-		U_PF = U[1:PF_velocities,]
-       	pf = PV("myPhaseField")
-		n = my_normal_directions[ii,]
-		heq = pf * EQ_h$feq
-		bounce = Bounce(U_PF)
-		sel = as.vector( (U_PF %*% n) < 0)
-
-		sel2 = as.vector( ( U_PF %*% n) == 0)	
-		exM = (h[sel2] - heq[sel2]) %*% rep(1,sum(sel2)) 
-		Unknowns = 1.0/sum(sel)
+    <?R
+    type_h = "fix_flux"
+    type = "pressure"
+    n = my_normal_directions[ii,]
+
+    EQ = MRT_eq(U,PV(1),u,mat=t(M))
+    w_g = subst(EQ$feq,U=0,V=0,W=0) # calculating weights as equilibrium with zero velocity
+    pstar = PV("pstar")
+    geq = pstar*w_g + (EQ$feq - w_g)
+    bounce = Bounce(U)
+    sel = as.vector( (U %*% n) < 0)
+
+    ng = g
+    ng[sel] = geq[sel] + (g-geq)[bounce][sel]
+
+    equations = V(sum(ng),ng %*% U)
+    rhs = V(pstar,u)
+
+    what_to_pull = c("pstar","U","V","W")
+    which_to_pull = 1:4
+    direction = which(n != 0)
+    if (type == "pressure") {
+        which_to_pull = which_to_pull[-(1+direction)] # throw out the equation related to this velocity direction
+        what_to_pull[1] = what_to_pull[1+direction] # pull out velocity from first equation
+    } else if (type == "velocity") {
+        which_to_pull = which_to_pull[-(1+direction)] # throw out the equation related to this velocity direction
+    } else stop("Unknown type of BC")
+
+    for (i in which_to_pull) {
+        C_pull(equations[i] - rhs[i],what_to_pull[i])
+    }
+    C(g[sel],ng[sel])
+
+    cat("real_t PF;\n")
+    pf = PV("PF")
+    desired_pf = PV("myPhaseField")
+    U_PF = U[1:PF_velocities,]
+    EQ_h = MRT_eq(U_PF,PV(1),u)
+
+    heq = pf * EQ_h$feq
+    bounce = Bounce(U_PF)
+    sel = as.vector( (U_PF %*% n) < 0)
+
+    nh = h
+    nh[sel] = heq[sel] + (h-heq)[bounce][sel]
+
+    equations = V(sum(nh),nh %*% U_PF)
+    rhs = V(desired_pf,desired_pf*u)
+
+    if (type_h == "fix_pf") {
+        C_pull(equations[1] - rhs[1], "PF" )
+    } else if (type_h == "open") {
+        C_pull(equations[1] - pf, "PF" )
+    } else {
+        nh = subst(nh, PF=desired_pf)
+        equations = subst(equations, PF=desired_pf)
+        if (type_h == "fix_flux") {
+        }  else if (type_h == "weird") {
+            what_to_pull = c("","U","V","W")
+            what_to_pull[1] = what_to_pull[1+direction] # pull out velocity from first equation
+            which_to_pull = (1:4)[-(1+direction)] # throw out equation
+            for (i in which_to_pull) {
+                C_pull(equations[i] - rhs[i],what_to_pull[i])
+            }
+        } else stop("Unknown type of BC")
+    }
 
-		C(h[sel], heq[sel] + heq[bounce][sel] - h[bounce][sel] - 0.5 * Unknowns * (exM))
-	?>
+    C(h[sel],nh[sel])
+    ?>
 }
 
 <?R
@@ -582,6 +624,14 @@ CudaDeviceFunction void calcPhaseGrad_init(){
 		}
 	<?R } ?>
 
+    // no need to do anything on the boundary
+    auto boundary_marker = NodeType & NODE_BOUNDARY;
+    if (boundary_marker == NODE_EPressure || boundary_marker == NODE_WPressure) {
+        gradPhiVal_x = 0.0;
+        gradPhiVal_y = 0.0;
+        gradPhiVal_z = 0.0;
+    }
+
     // simply copy the values
 	gradPhi_PhaseF = PhaseF(0,0,0);
 }

From e20e665eca974fd280c0f1e4dc6001bc0a115dfa Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 28 May 2024 08:42:08 +1000
Subject: [PATCH 10/22] Rewrite fixes generically

---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 92 +++++++++++++------
 .../multiphase/d3q27_pf_velocity/Dynamics.R   | 10 +-
 .../d3q27_pf_velocity/Dynamics.c.Rt           | 31 +++++--
 3 files changed, 93 insertions(+), 40 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index e580d3fc6..f7e784f3f 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -70,55 +70,72 @@ CudaDeviceFunction void <?R C(my_velocity_boundaries[ii])?>(){
 		C(h[sel], heq[sel] + heq[bounce][sel] - h[bounce][sel] - 0.5 * Unknowns * (exM))
 	?>
 }
+<?R
+}
+?>
 
-CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
+<?R
+for (ii in 1:cube_faces){
+for (j in 1:length(pressure_boundary_types)) {  ?>
+CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boundary_types[j]) ?>(){
     real_t d = getRho();
-
     real_t myPhaseField = 0.0;
-    // TODO: Need a proper fix soon
-    if ((NodeType & NODE_BOUNDARY) == NODE_WPressure){
-        d = Density_l;
-        myPhaseField = 0.0;
+
+    if ((NodeType & NODE_BOUNDARY) == <?%s paste0("NODE_WPressure", pressure_boundary_types[j]) ?>){
+        d = InvasionDrainage == 2? Density_l : Density_h;
+        myPhaseField = InvasionDrainage == 2? 0.0 : 1.0;
     }
 
-    if ((NodeType & NODE_BOUNDARY) == NODE_EPressure){
-        d = Density_h;
-        myPhaseField = 1.0;
+    if ((NodeType & NODE_BOUNDARY) == <?%s paste0("NODE_EPressure", pressure_boundary_types[j]) ?>){
+       d = InvasionDrainage == 2? Density_h : Density_l;
+       myPhaseField = InvasionDrainage == 2? 1.0 : 0.0;
     }
 
 	real_t pstar = Pressure / (d*cs2);
-
     <?R
-    type_h = "fix_flux"
+    type_h = (if (j == 1) "fpf" else pressure_boundary_types[j])
     type = "pressure"
     n = my_normal_directions[ii,]
 
     EQ = MRT_eq(U,PV(1),u,mat=t(M))
-    w_g = subst(EQ$feq,U=0,V=0,W=0) # calculating weights as equilibrium with zero velocity
+    # calculating weights as equilibrium with zero velocity
+    w_g = subst(EQ$feq,U=0,V=0,W=0)
     pstar = PV("pstar")
     geq = pstar*w_g + (EQ$feq - w_g)
     bounce = Bounce(U)
     sel = as.vector( (U %*% n) < 0)
-
+    # n prefix for newly expressed values
     ng = g
     ng[sel] = geq[sel] + (g-geq)[bounce][sel]
 
-    equations = V(sum(ng),ng %*% U)
-    rhs = V(pstar,u)
+    # write down the closure equation
+    # sum of LBM densities equal to pressure
+    # the first moment equal to the velocity
+    # (assuming no forcing - no phase-field gradient)!
+    equations_lhs = V(sum(ng),ng %*% U)
+    equations_rhs = V(pstar,u)
 
     what_to_pull = c("pstar","U","V","W")
     which_to_pull = 1:4
     direction = which(n != 0)
     if (type == "pressure") {
-        which_to_pull = which_to_pull[-(1+direction)] # throw out the equation related to this velocity direction
+        # pop the pressure out of the list of equation
+        # and don't duplicate
+        which_to_pull = which_to_pull[-(1+direction)]
         what_to_pull[1] = what_to_pull[1+direction] # pull out velocity from first equation
     } else if (type == "velocity") {
+        # pop out the equation related to this velocity direction
+        # from the list of equations, as it is already know
         which_to_pull = which_to_pull[-(1+direction)] # throw out the equation related to this velocity direction
     } else stop("Unknown type of BC")
 
+    # express all the unknowns in terms of knowns
     for (i in which_to_pull) {
-        C_pull(equations[i] - rhs[i],what_to_pull[i])
+        C_pull(equations_lhs[i] - equations_rhs[i],what_to_pull[i])
     }
+
+    # now that the quantities are known, can properly calculate the distributution
+    # functions and reflect
     C(g[sel],ng[sel])
 
     cat("real_t PF;\n")
@@ -127,6 +144,7 @@ CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
     U_PF = U[1:PF_velocities,]
     EQ_h = MRT_eq(U_PF,PV(1),u)
 
+    # assuming no force shift here
     heq = pf * EQ_h$feq
     bounce = Bounce(U_PF)
     sel = as.vector( (U_PF %*% n) < 0)
@@ -134,23 +152,35 @@ CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
     nh = h
     nh[sel] = heq[sel] + (h-heq)[bounce][sel]
 
-    equations = V(sum(nh),nh %*% U_PF)
-    rhs = V(desired_pf,desired_pf*u)
+    # write down the closure equation
+    # sum of phase-field densities equal to phase-field value (desired or imposed)
+    # the first moment equal to the desired phase-field momentum
+    equations_lhs = V(sum(nh),nh %*% U_PF)
+    equations_rhs = V(desired_pf,desired_pf*u)
 
-    if (type_h == "fix_pf") {
-        C_pull(equations[1] - rhs[1], "PF" )
+    if (type_h == "fpf") {
+        # Make sure the sum of h is equal to the desired phase-field
+        # - "interface" will stop on the boundary
+        C_pull(equations_lhs[1] - equations_rhs[1], "PF" )
     } else if (type_h == "open") {
-        C_pull(equations[1] - pf, "PF" )
+        # The sum of phase-field is not enforced - the same is kept
+        # - "interface" is passing through the boundary
+        C_pull(equations_lhs[1] - pf, "PF" )
     } else {
+        # desired phase-field is enforced via substitution in
+        # the equilibrium distribution function directly
         nh = subst(nh, PF=desired_pf)
-        equations = subst(equations, PF=desired_pf)
-        if (type_h == "fix_flux") {
+        #
+        lhs_equations = subst(equations_lhs, PF=desired_pf)
+        if (type_h == "fflux") {
+        # we are pretty much done
         }  else if (type_h == "weird") {
             what_to_pull = c("","U","V","W")
             what_to_pull[1] = what_to_pull[1+direction] # pull out velocity from first equation
             which_to_pull = (1:4)[-(1+direction)] # throw out equation
+            # express velocities through the phase-field?
             for (i in which_to_pull) {
-                C_pull(equations[i] - rhs[i],what_to_pull[i])
+                C_pull(equations_lhs[i] - equations_rhs[i],what_to_pull[i])
             }
         } else stop("Unknown type of BC")
     }
@@ -160,7 +190,7 @@ CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
 }
 
 <?R
-}
+    }}
 ?>
 
 #ifdef OPTIONS_OutFlow
@@ -953,13 +983,17 @@ CudaDeviceFunction real_t Init_real_wallNorm()
 		int i,j,k;
 	  	real_t tmp = 0.0;
         real_t PhaseF_tmp[3][3][3];
-        // mirror in the beginning, in case of non-periodic boundaries
+
+        // find in which direction to do reflection
+        // mirror in the beginning, in case of non-periodic boundaries,
+        // so we don't accidentally take the value from the other phase
+        /// Can keep if the setting set to true
         int nx = constContainer.nx;
         for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
-            if ((i == -1) && (X == 0.5)) {
+            if (InvasionDrainage && (i == -1) && (X == 0.5)) {
                 PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
             }
-            else if ((i == 1) && (X == nx - 0.5)) {
+            else if (InvasionDrainage && (i == 1) && (X == nx - 0.5)) {
                 PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
             }
             else {
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 4d0be761d..6b95173fa 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -172,7 +172,7 @@ if (Options$thermo){
         AddAction("Iteration", c("BaseIter", "calcPhase",  calcGrad, "calcWall_CA", "calcWallPhase_correction"))
 	    AddAction("Init"     , c("PhaseInit","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init"  , "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
 	    AddAction("InitFields"     , c("InitFromFieldsStage","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
-	    AddAction("InitFieldsWithNormals"     , c("InitFromFieldsStageWithNormals","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
+	    AddAction("InitFieldsWithNormals"     , c("InitFromFieldsStageWithNormals",,"WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
     } else {
 		AddAction("Iteration", c("BaseIter", "calcPhase", "calcWall", "calcWallPhase_correction"))
 		AddAction("Init"     , c("PhaseInit", "WallInit_Real", "WallInit" , "calcWall","calcWallPhase_correction", "BaseInit"))
@@ -248,6 +248,7 @@ if (Options$thermo){
 	AddSetting(name="VelocityY", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="VelocityZ", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="Pressure" , default=0.0, comment='inlet/outlet/init density', zonal=T)
+    AddSetting(name='InvasionDrainage', default=0, comment='0 nothing, 1 flooding, -1 drainage')
 	AddSetting(name="GravitationX", default=0.0, comment='applied (rho)*GravitationX')
 	AddSetting(name="GravitationY", default=0.0, comment='applied (rho)*GravitationY')
 	AddSetting(name="GravitationZ", default=0.0, comment='applied (rho)*GravitationZ')
@@ -276,12 +277,17 @@ if (Options$thermo){
 ##########################
 	AddNodeType("Smoothing",group="ADDITIONALS")
 	AddNodeType(name="flux_nodes", group="ADDITIONALS")
+    # the first one are "fix_pf"
+    pressure_boundary_types = c("", "open", "fflux")
 	dotR_my_velocity_boundaries = paste0(c("N","E","S","W","F","B"),"Velocity")
-    dotR_my_pressure_boundaries = paste0(c("N","E","S","W","F","B"),"Pressure")
+    dotR_my_pressure_boundaries = outer(paste0(c("N","E","S","W","F","B"),"Pressure"), pressure_boundary_types, FUN  = paste, sep="")
     for (ii in 1:6){
         AddNodeType(name=dotR_my_velocity_boundaries[ii], group="BOUNDARY")
+    }
+    for (ii in 1:length(dotR_my_pressure_boundaries)) {
         AddNodeType(name=dotR_my_pressure_boundaries[ii], group="BOUNDARY")
     }
+
 	AddNodeType(name="MovingWall_N", group="BOUNDARY")
 	AddNodeType(name="MovingWall_S", group="BOUNDARY")
 	AddNodeType(name="Solid", group="BOUNDARY")
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
index f2b632c7d..3a2ae382d 100755
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
@@ -149,7 +149,8 @@ CudaDeviceFunction vector_t calcGradPhi()
 		    ?>
 
         auto boundary_marker = NodeType & NODE_BOUNDARY;
-        if (boundary_marker == NODE_EPressure || boundary_marker == NODE_WPressure) {
+        <?R condition = paste(paste0("boundary_marker == NODE_", c(dotR_my_velocity_boundaries, dotR_my_pressure_boundaries), collapse = ' || ')) ?>
+        if (<?%s condition ?>) {
             gradPhi.x = 0.0;
             gradPhi.y = 0.0;
             gradPhi.z = 0.0;
@@ -235,7 +236,8 @@ CudaDeviceFunction real_t calcMu(real_t C)
 	#else
 
         auto boundary_marker = NodeType & NODE_BOUNDARY;
-        if (boundary_marker == NODE_EPressure || boundary_marker == NODE_WPressure) {
+        <?R condition = paste(paste0("boundary_marker == NODE_", c(dotR_my_velocity_boundaries, dotR_my_pressure_boundaries), collapse = ' || ')) ?>
+        if (<?%s condition ?>) {
             lpPhi = 0.0;
         }
 
@@ -274,12 +276,17 @@ CudaDeviceFunction real_t calcF_phi(int i, real_t tmp1, real_t nx, real_t ny, re
 	return f_phi;
 }
 
+CudaDeviceFunction void InitPhaseFforWallNormals()
+{
+	if ( IamWall || IamSolid ) PhaseF = -999;
+}
+
 /*	INITIALISATION:	*/
 CudaDeviceFunction void Init() 
 {
 	PhaseF = PhaseField;
 	specialCases_Init();
-	if ( IamWall || IamSolid ) PhaseF = -999;
+    InitPhaseFforWallNormals();
 
 	if (developedFlow > 0.1) {
 		U = 6.0 * Uavg * Y*(HEIGHT - Y)/(HEIGHT*HEIGHT);
@@ -302,7 +309,7 @@ CudaDeviceFunction void InitFromFieldsStage()
 	U = Init_UX_External;
 	V = Init_UY_External;
 	W = Init_UZ_External;
-	if ( IamWall || IamSolid ) PhaseF = -999;
+    InitPhaseFforWallNormals();
 	pnorm = 0.0; // initialise as zero and fill in later stage
 }
 
@@ -313,7 +320,7 @@ CudaDeviceFunction void InitFromFieldsStageWithNormals()
 	U = Init_UX_External;
 	V = Init_UY_External;
 	W = Init_UZ_External;
-	if ( IamWall || IamSolid ) PhaseF = -999;
+    InitPhaseFforWallNormals();
 	pnorm = 0.0; // initialise as zero and fill in later stage
     nw_x = Init_nwx_external;
     nw_y = Init_nwy_external;
@@ -719,16 +726,22 @@ CudaDeviceFunction void updateBoundary(){
 			BounceBack();
 			break;
         <?R
-        for (ii in 1:6){
+        for (ii in 1:cube_faces){
         ?>
         case NODE_<?R C(my_velocity_boundaries[ii])?>:
 			<?R C(my_velocity_boundaries[ii])?>();
 			break;
-		case NODE_<?R C(my_pressure_boundaries[ii])?>:
-			<?R C(my_pressure_boundaries[ii])?>();
-			break;
 		<?R
         }
+        ?>
+        <?R
+        for (my_pressure_boundary in dotR_my_pressure_boundaries) {
+        ?>
+		case NODE_<?%s my_pressure_boundary ?>:
+			<?%s my_pressure_boundary ?>();
+			break;
+		<?R
+              }
         ?>
 		case NODE_MovingWall_N:
 			MovingNWall();

From fe4f08b04f98fa707bf57101546776fd5422c478 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 28 May 2024 09:01:35 +1000
Subject: [PATCH 11/22] Add comments

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index f7e784f3f..6af24b6d6 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -987,7 +987,12 @@ CudaDeviceFunction real_t Init_real_wallNorm()
         // find in which direction to do reflection
         // mirror in the beginning, in case of non-periodic boundaries,
         // so we don't accidentally take the value from the other phase
-        /// Can keep if the setting set to true
+        // ===========================================================
+        // NOTE: Works only in X direction. To use pressure/velocity
+        // bondaries in other directions need to mark the boundaries
+        // with symmetry market, otherwise wall gradients might point
+        // to the other side. (I.e. SymmetryX_plus at the same location
+        // as EPressure, etc)
         int nx = constContainer.nx;
         for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
             if (InvasionDrainage && (i == -1) && (X == 0.5)) {

From 8c48c8df35dcf82af24b9fb16cc1d6f439073bb3 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 28 May 2024 12:47:15 +1000
Subject: [PATCH 12/22] Make more explicit

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 6af24b6d6..fcdfa5c85 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -995,10 +995,10 @@ CudaDeviceFunction real_t Init_real_wallNorm()
         // as EPressure, etc)
         int nx = constContainer.nx;
         for (i=-1;i<2;i++){for (j=-1;j<2;j++){for (k=-1;k<2;k++){
-            if (InvasionDrainage && (i == -1) && (X == 0.5)) {
+            if ((InvasionDrainage > 0) && (i == -1) && (X == 0.5)) {
                 PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
             }
-            else if (InvasionDrainage && (i == 1) && (X == nx - 0.5)) {
+            else if ((InvasionDrainage > 0) && (i == 1) && (X == nx - 0.5)) {
                 PhaseF_tmp[i+1][j+1][k+1] = PhaseF_dyn(-i,j,k);
             }
             else {

From 34ea4257420c2dc38101b8f39292aed17c696131 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 28 May 2024 13:05:28 +1000
Subject: [PATCH 13/22] Fix

---
 models/multiphase/d3q27_pf_velocity/Dynamics.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 6b95173fa..8629d5de7 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -172,7 +172,7 @@ if (Options$thermo){
         AddAction("Iteration", c("BaseIter", "calcPhase",  calcGrad, "calcWall_CA", "calcWallPhase_correction"))
 	    AddAction("Init"     , c("PhaseInit","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init"  , "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
 	    AddAction("InitFields"     , c("InitFromFieldsStage","WallInit_Real","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
-	    AddAction("InitFieldsWithNormals"     , c("InitFromFieldsStageWithNormals",,"WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
+	    AddAction("InitFieldsWithNormals"     , c("InitFromFieldsStageWithNormals","WallInit_CA" , "calcPhaseGrad_init", "calcWall_CA", "calcWallPhase_correction", "BaseInit"))
     } else {
 		AddAction("Iteration", c("BaseIter", "calcPhase", "calcWall", "calcWallPhase_correction"))
 		AddAction("Init"     , c("PhaseInit", "WallInit_Real", "WallInit" , "calcWall","calcWallPhase_correction", "BaseInit"))

From 98accf45f96e0cde59bf8347faf643aab76fe716 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 4 Jun 2024 08:55:16 +1000
Subject: [PATCH 14/22] Fix

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index fcdfa5c85..a672e94e4 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -972,7 +972,7 @@ CudaDeviceFunction real_t getSpecialBoundaryPoint() {
 
 /* Calculate the actual wall norm based on the solid
  *field */
-CudaDeviceFunction real_t Init_real_wallNorm()
+CudaDeviceFunction void Init_real_wallNorm()
 {
     // set initial values
     nw_x = 0.0;

From 660ce39361c9c0d7f6a31dea9e42202953e03a25 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Wed, 5 Jun 2024 09:07:27 +1000
Subject: [PATCH 15/22] Change the default type

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 2 +-
 models/multiphase/d3q27_pf_velocity/Dynamics.R    | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index a672e94e4..65942baa1 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -93,7 +93,7 @@ CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boun
 
 	real_t pstar = Pressure / (d*cs2);
     <?R
-    type_h = (if (j == 1) "fpf" else pressure_boundary_types[j])
+    type_h = (if (j == 1) "fflux" else pressure_boundary_types[j])
     type = "pressure"
     n = my_normal_directions[ii,]
 
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 8629d5de7..45a80d2bb 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -277,8 +277,8 @@ if (Options$thermo){
 ##########################
 	AddNodeType("Smoothing",group="ADDITIONALS")
 	AddNodeType(name="flux_nodes", group="ADDITIONALS")
-    # the first one are "fix_pf"
-    pressure_boundary_types = c("", "open", "fflux")
+    # the first one are "fflux"
+    pressure_boundary_types = c("", "open", "fpf")
 	dotR_my_velocity_boundaries = paste0(c("N","E","S","W","F","B"),"Velocity")
     dotR_my_pressure_boundaries = outer(paste0(c("N","E","S","W","F","B"),"Pressure"), pressure_boundary_types, FUN  = paste, sep="")
     for (ii in 1:6){

From b4696c375728d11db6893b9e801db81bb2239bcb Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Fri, 14 Jun 2024 11:19:21 +1000
Subject: [PATCH 16/22] Add setting pressure on the boundaries using RunR

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 4 ++++
 models/multiphase/d3q27_pf_velocity/Dynamics.R    | 4 ++++
 2 files changed, 8 insertions(+)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 65942baa1..ef80fcd11 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -92,6 +92,10 @@ CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boun
     }
 
 	real_t pstar = Pressure / (d*cs2);
+    if (UseExternalPressure > 0){
+        /// use pressure set externally from RunR
+        p_star = Pressure_external / (d*cs2);
+    }
     <?R
     type_h = (if (j == 1) "fflux" else pressure_boundary_types[j])
     type = "pressure"
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 45a80d2bb..314c1ab05 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -19,6 +19,9 @@ AddDensity(name="Init_nwx_external", group="init_normals", dx=0,dy=0,dz=0, param
 AddDensity(name="Init_nwy_external", group="init_normals", dx=0,dy=0,dz=0, parameter=TRUE)
 AddDensity(name="Init_nwz_external", group="init_normals", dx=0,dy=0,dz=0, parameter=TRUE)
 
+# Add extra density for setting the pressure on boundaries
+AddDensity(name="Pressure_external", group="Vel", dx=0, dy=0, dz=0, parameter=TRUE)
+
 # macroscopic params
 # - consider migrating to fields
 AddDensity(name="pnorm", dx=0, dy=0, dz=0, group="Vel")
@@ -207,6 +210,7 @@ if (Options$thermo){
 	AddSetting(name="omega_phi", comment='one over relaxation time (phase field)')
 	AddSetting(name="M", omega_phi='1.0/(3*M+0.5)', default=0.02, comment='Mobility')
 	AddSetting(name="sigma", comment='surface tension')
+    AddSetting(name="UseExternalPressure", default=0, comment='Use external pressure (set from RunR)')
 	AddSetting(name="force_fixed_iterator", default=2, comment='to resolve implicit relation of viscous force')
   	AddSetting(name="Washburn_start", default="0", comment='Start of washburn gas phase')
   	AddSetting(name="Washburn_end", default="0", comment='End of washburn gas phase')

From 9a0f8d600279f8aeee67540b43f425497c81f05f Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Fri, 14 Jun 2024 12:53:20 +1000
Subject: [PATCH 17/22] Fix typo

---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index ef80fcd11..0a9aae192 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -94,7 +94,7 @@ CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boun
 	real_t pstar = Pressure / (d*cs2);
     if (UseExternalPressure > 0){
         /// use pressure set externally from RunR
-        p_star = Pressure_external / (d*cs2);
+        pstar = Pressure_external / (d*cs2);
     }
     <?R
     type_h = (if (j == 1) "fflux" else pressure_boundary_types[j])

From 5ca385bc8fd96ef02327371a9996e06e42b5cdb6 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Mon, 26 Aug 2024 14:09:22 +1000
Subject: [PATCH 18/22] Fix normalization for the special boundary points

Need to normalise by the sum of weights obviously
(otherwise would leads to crazy low values of the phase field
in random points of the domain)
---
 models/multiphase/d3q27_pf_velocity/Boundary.c.Rt | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 0a9aae192..7cfef71a9 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -781,17 +781,19 @@ CudaDeviceFunction void calcWallPhase_correction() {
             int count = 0;
             real_t fluid_average = 0.0;
             real_t x, y, z;
+            real_t denom = 0.0;
 			for (int i=1; i<27 ; i++){
 				x = d3q27_ex[i];
 				y = d3q27_ey[i];
 				z = d3q27_ez[i];
                 if (!IsBoundary_dyn(int(x), int(y), int(z)) && PhaseF_dyn(int(x), int(y), int(z)) > -0.5) {
                     fluid_average += wg[i]*PhaseF_dyn(int(x), int(y), int(z));
+                    denom += wg[i];
                     count++;
                 }
 			}
             // print sum and average and count
-            PhaseF = count == 0 ? 1 : fluid_average / count;
+            PhaseF = count == 0 ? 1 : fluid_average / denom;
         } else {
             PhaseF = point_in_the_normal_direction;
         }

From 33a2b7f2f8b2cfd7fb793d842e88f87354a0d056 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Thu, 29 Aug 2024 11:18:43 +1000
Subject: [PATCH 19/22] Add handling of errors arising in RunPython

---
 src/Handlers/cbRunR.cpp | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/src/Handlers/cbRunR.cpp b/src/Handlers/cbRunR.cpp
index f09b7277e..5c70dfa30 100644
--- a/src/Handlers/cbRunR.cpp
+++ b/src/Handlers/cbRunR.cpp
@@ -692,10 +692,16 @@ namespace RunPython {
 	void parseEval(const std::string& source) {
 		if (! py_initialised) initializePy();
 		Rcpp::Function py_run_string("py_run_string");
-		py_run_string(source);
+        py_run_string(source);
 		return;
 	}	
 
+    std::tuple<std::string, std::string> getLastError() {
+        Rcpp::Function py_last_error("py_last_error");
+        Rcpp::List error_details = py_last_error();
+        return std::make_tuple(Rcpp::as<std::string>(error_details["type"]), Rcpp::as<std::string>(error_details["value"]));
+    } 
+
 	void initializePy() {
 		if (py_initialised) return;
 		RInside& R = RunR::GetR();
@@ -836,8 +842,17 @@ int cbRunR::DoIt() {
 		ERROR("Caught std exception: %s", ex.what());
 		return -1;
 	} catch (...) {
-		ERROR("Caught uknown exception");
-		return -1;
+        if (python) { 
+            std::string error_type, error_message;
+            std::tie(error_type, error_message) = RunPython::getLastError();
+            if (error_type == "SimulationStopSignal") {
+		        WARNING("From RunPython received signal to stop simulation with message: %s. Stopping the enclosing handler.", 
+                    error_message.c_str());
+                return 1;  
+            } 
+        }
+        ERROR("Caught uknown exception");
+        return -1;
 	}
 	return 0;
 }

From d69ce1fa58071befb67f53c3f6e4e63b09f553c8 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Thu, 5 Sep 2024 15:54:07 +1000
Subject: [PATCH 20/22] Improve the phase field and density on Boundaries

---
 .../multiphase/d3q27_pf_velocity/Boundary.c.Rt   | 16 ++++++----------
 models/multiphase/d3q27_pf_velocity/Dynamics.R   |  2 +-
 2 files changed, 7 insertions(+), 11 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 7cfef71a9..15467d732 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -78,18 +78,14 @@ CudaDeviceFunction void <?R C(my_velocity_boundaries[ii])?>(){
 for (ii in 1:cube_faces){
 for (j in 1:length(pressure_boundary_types)) {  ?>
 CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boundary_types[j]) ?>(){
-    real_t d = getRho();
-    real_t myPhaseField = 0.0;
 
-    if ((NodeType & NODE_BOUNDARY) == <?%s paste0("NODE_WPressure", pressure_boundary_types[j]) ?>){
-        d = InvasionDrainage == 2? Density_l : Density_h;
-        myPhaseField = InvasionDrainage == 2? 0.0 : 1.0;
-    }
 
-    if ((NodeType & NODE_BOUNDARY) == <?%s paste0("NODE_EPressure", pressure_boundary_types[j]) ?>){
-       d = InvasionDrainage == 2? Density_h : Density_l;
-       myPhaseField = InvasionDrainage == 2? 1.0 : 0.0;
-    }
+    // can't use PhaseF or getRho() because
+    // the values can be streamed from the other side of the
+    // boundary
+    real_t myPhaseField = PhaseField;
+	real_t d =  Density_l + (Density_h-Density_l) * (PhaseField - PhaseField_l)/(PhaseField_h - PhaseField_l);
+
 
 	real_t pstar = Pressure / (d*cs2);
     if (UseExternalPressure > 0){
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 314c1ab05..71f4bf875 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -252,7 +252,7 @@ if (Options$thermo){
 	AddSetting(name="VelocityY", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="VelocityZ", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="Pressure" , default=0.0, comment='inlet/outlet/init density', zonal=T)
-    AddSetting(name='InvasionDrainage', default=0, comment='0 nothing, 1 flooding, -1 drainage')
+    AddSetting(name='InvasionDrainage', default=0, comment='0 nothing, anything bigger is invasion/drainage case")
 	AddSetting(name="GravitationX", default=0.0, comment='applied (rho)*GravitationX')
 	AddSetting(name="GravitationY", default=0.0, comment='applied (rho)*GravitationY')
 	AddSetting(name="GravitationZ", default=0.0, comment='applied (rho)*GravitationZ')

From b34c5b460c3d60a9c74606c5ae70a2b18383a869 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Wed, 11 Sep 2024 14:40:13 +1000
Subject: [PATCH 21/22] Fix typo

---
 models/multiphase/d3q27_pf_velocity/Dynamics.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index 71f4bf875..d5780da1d 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -252,7 +252,7 @@ if (Options$thermo){
 	AddSetting(name="VelocityY", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="VelocityZ", default=0.0, comment='inlet/outlet/init velocity', zonal=T)
 	AddSetting(name="Pressure" , default=0.0, comment='inlet/outlet/init density', zonal=T)
-    AddSetting(name='InvasionDrainage', default=0, comment='0 nothing, anything bigger is invasion/drainage case")
+    AddSetting(name='InvasionDrainage', default=0, comment="0 nothing, anything bigger is invasion/drainage case")
 	AddSetting(name="GravitationX", default=0.0, comment='applied (rho)*GravitationX')
 	AddSetting(name="GravitationY", default=0.0, comment='applied (rho)*GravitationY')
 	AddSetting(name="GravitationZ", default=0.0, comment='applied (rho)*GravitationZ')

From 9d9e58360f0b700e7d07d1621392f7e66a412f32 Mon Sep 17 00:00:00 2001
From: Dmytro Sashko <dmshko@gmail.com>
Date: Tue, 17 Sep 2024 15:25:18 +1000
Subject: [PATCH 22/22] Add more robust version of special boundary points

---
 .../d3q27_pf_velocity/Boundary.c.Rt           | 53 +++++++++----------
 1 file changed, 25 insertions(+), 28 deletions(-)

diff --git a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
index 15467d732..c300fe24b 100644
--- a/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Boundary.c.Rt
@@ -765,35 +765,32 @@ CudaDeviceFunction void calcPhaseGradCloseToBoundary(){
  * yet set value
  */
 CudaDeviceFunction void calcWallPhase_correction() {
-	PhaseF = PhaseF(0,0,0);
-	if (IsSpecialBoundaryPoint == <?%s NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?> ) {
-		// take the phase field calculated already from the node in front.
-		// Might as well calculate using the neighbors, e.g. averaging
-		real_t point_in_the_normal_direction = PhaseF_dyn(nw_x, nw_y, nw_z);
-        // in some cases we might (TODO: investigate) still point into the solid node
-        // that has normal pointing into another, we need to avoid this
-        if (point_in_the_normal_direction == <?%f SPECIAL_POINT_HUGE_MAGIC_NUMBER ?>) {
-            // ignore first node
-            int count = 0;
-            real_t fluid_average = 0.0;
-            real_t x, y, z;
-            real_t denom = 0.0;
-			for (int i=1; i<27 ; i++){
-				x = d3q27_ex[i];
-				y = d3q27_ey[i];
-				z = d3q27_ez[i];
-                if (!IsBoundary_dyn(int(x), int(y), int(z)) && PhaseF_dyn(int(x), int(y), int(z)) > -0.5) {
-                    fluid_average += wg[i]*PhaseF_dyn(int(x), int(y), int(z));
-                    denom += wg[i];
-                    count++;
-                }
-			}
-            // print sum and average and count
-            PhaseF = count == 0 ? 1 : fluid_average / denom;
-        } else {
-            PhaseF = point_in_the_normal_direction;
+    PhaseF = PhaseF(0,0,0);
+    if (IsSpecialBoundaryPoint == <?%s NORMAL_POINTING_INTO_SOLID_ON_NEXT_NODE ?> ) {
+        // We could take the phase field calculated already from the node in front
+        // (solid in the direction of solid normal)
+        // if it is already has the phase field value calculated (from the BC application)
+        // however this seems to be less stable...
+        // More robust is calculate using the neighbors, e.g. averaging.
+        // If this also would not work at some point, more detailed investigation is needed
+        int count = 0;
+        real_t fluid_average = 0.0;
+        real_t x, y, z;
+        real_t denom = 0.0;
+        // ignore first node
+        for (int i=1; i<27 ; i++){
+            x = d3q27_ex[i];
+            y = d3q27_ey[i];
+            z = d3q27_ez[i];
+            if (!IsBoundary_dyn(int(x), int(y), int(z)) && PhaseF_dyn(int(x), int(y), int(z)) > -0.5) {
+                fluid_average += wg[i]*PhaseF_dyn(int(x), int(y), int(z));
+                denom += wg[i];
+                count++;
+            }
         }
-	}
+        // print sum and average and count
+        PhaseF = count == 0 ? 1 : fluid_average / denom;
+    }
 }