merged upstream master

CFD-GO · Apr 29, 2024 · 39adcdf · 39adcdf
2 parents 0402bf2 + 2362a9d
commit 39adcdf
Show file tree

Hide file tree

Showing 9 changed files with 524 additions and 44 deletions.
diff --git a/models/particles/d3q27_PSM/Dynamics.R b/models/particles/d3q27_PSM/Dynamics.R
@@ -71,7 +71,7 @@ AddQuantity(name="U",unit="m/s",vector=T)
 AddQuantity(name="Rho",unit="kg/m3")
 
 AddSetting(name="omegaF", comment='one over F relaxation time and initial relaxation time for kl')
-AddSetting(name="nu", omegaF='1.0/(3*nu+0.5)', default=0.1, comment='kinetic viscosity in LBM unit')
+AddSetting(name="nu", omegaF='1.0/(3*nu+0.5)', default=0.1, comment='kinetic viscosity in LBM unit', unit='m2/s')
 
 if (Options$TRT) {
     AddSetting( name="Lambda", comment="TRT Magic Number")
@@ -120,6 +120,8 @@ AddNodeType(name="BPressure", group="BOUNDARY")
 
 AddNodeType(name="MovingWall_N", group="BOUNDARY")
 AddNodeType(name="MovingWall_S", group="BOUNDARY")
+AddNodeType(name="MovingWall_E", group="BOUNDARY")
+AddNodeType(name="MovingWall_W", group="BOUNDARY")
 
 
 if (Options$particles) {

diff --git a/models/particles/d3q27_PSM/Dynamics.c.Rt b/models/particles/d3q27_PSM/Dynamics.c.Rt
@@ -27,7 +27,6 @@ Code updates:
 
 */
 
-#define SQRT_3 1.732050807568877
 #define PI 3.141592653589793
 
 <?R
@@ -178,6 +177,13 @@ CudaDeviceFunction void Run() {
         	MovingSWall();
                 break;
 
+	case NODE_MovingWall_E:
+        	MovingEWall();
+                break;
+
+	case NODE_MovingWall_W:
+        	MovingWWall();
+                break;
 	}	   
 
 	switch (NodeType & NODE_COLLISION) {
@@ -393,7 +399,7 @@ CudaDeviceFunction void PressureExitB()
 
 CudaDeviceFunction void MovingNWall()
 {
-        f[1]  = f[3];
+        f[4]  = f[3];
         f[16] = f[17];
         f[18] = f[15];
 
@@ -420,6 +426,35 @@ CudaDeviceFunction void MovingSWall(){
         f[7] = f[10] + VelocityX/9.0;
 }
 
+CudaDeviceFunction void MovingEWall()
+{
+        f[2]  = f[1];
+        f[8] = f[9];
+        f[10] = f[7];
+
+        f[20]  = f[25] + VelocityZ/36.0;
+        f[22] = f[23]  + VelocityZ/36.0;
+        f[12] = f[13] + VelocityZ/9.0;
+
+        f[24] = f[21] - VelocityZ/36.0;
+        f[26]= f[19]  - VelocityZ/36.0;
+        f[14] = f[11] - VelocityZ/9.0;
+}
+
+CudaDeviceFunction void MovingWWall()
+{
+        f[1]  = f[2];
+        f[9] = f[8];
+        f[7] = f[10];
+
+        f[25]  = f[20] - VelocityZ/36.0;
+        f[23] = f[22]  - VelocityZ/36.0;
+        f[13] = f[12] - VelocityZ/9.0;
+
+        f[21] = f[24] + VelocityZ/36.0;
+        f[19]= f[26]  + VelocityZ/36.0;
+        f[11] = f[14] + VelocityZ/9.0;
+}
 
 CudaDeviceFunction void CalcEquilibrium(real_t f_tab[27], real_t d, real_t u[3])
 {
@@ -450,11 +485,59 @@ CudaDeviceFunction void CalcF() {
 	C( uP, 0)
 	C( sol, 0)
 ?>
+
+<?R if (Options$MS) { ?>
+
+	for (auto p : SyncParticleIterator(X,Y,Z)) {
+
+        real_t dist = sqrt(p.diff.x*p.diff.x + p.diff.y*p.diff.y + p.diff.z*p.diff.z);
+
+        if ((dist - p.rad) < 0.5) {
+
+            real_t localCoverage = 0.0;
+
+            if ((dist - p.rad) < -1.0){
+                localCoverage = 1.0;
+            } else{
+                localCoverage = (p.rad -  0.084/p.rad + 0.5 - dist);
+            }
+
+            if (localCoverage > 1.0){ localCoverage =  1.0;}
+            if (localCoverage < 0.0){ localCoverage =  0.0;}
+
+            <?R if (Options$KL) { ?>
+
+                real_t omegaF_tmp = 1.0/(3.0*nu_app+0.5); //From previous timestep
+                localCoverage = localCoverage*(1.0/omegaF_tmp - 0.5) / ((1 - localCoverage) + 1.0/omegaF_tmp - 0.5);
+
+            <?R } else { ?>
+
+                localCoverage = localCoverage*(1.0/omegaF - 0.5) / ((1 - localCoverage) + 1.0/omegaF - 0.5); 
+
+            <?R } ?>
+
+<?R
+            localCoverage = PV("localCoverage")
+            C( sol, sol + localCoverage )
+?>
+        }
+
+    }
+
+    real_t sol_factor = 1.0; 
+
+	<?R sol_factor = PV("sol_factor") ?>
+    if (sol > 1.0) {
+        sol_factor = 1.0/sol;
+    } 
+
+<?R } ?>
+
 	for (auto p : SyncParticleIterator(X,Y,Z)) {
 
         real_t dist = sqrt(p.diff.x*p.diff.x + p.diff.y*p.diff.y + p.diff.z*p.diff.z);
 
-        if ((dist - p.rad)<SQRT_3) {
+        if ((dist - p.rad) < 0.5) {
 
             real_t localCoverage = 0.0;
 
@@ -483,21 +566,42 @@ CudaDeviceFunction void CalcF() {
             uparticle = PV("p.cvel.",c("x","y","z"))
             force = PV("force.",c("x","y","z"))
             localCoverage = PV("localCoverage")
-            C( force, d*localCoverage*(u - uparticle) )
-?>
-            p.applyForce(force);
-<?R
-            C( sol, sol + localCoverage )
+
+            if (Options$MS) {
+
+                C( force, d*sol_factor*localCoverage*(u - uparticle) )
+
+            } else {
+
+                C( force, d*localCoverage*(u - uparticle) )
+                C( sol, sol + localCoverage )
+
+            }
+
             C( uP, uP + localCoverage*uparticle )
 ?>
+            p.applyForce(force);
         }
 
     }
 
-	if (sol > 1.0) { sol =  1.0; }
+<?R if (Options$MS) { ?>
+
 	if (sol > 1.0e-8) {
-<?R C( uP, uP*sol^{-1} ) ?>
-	} else {
+        <?R C( uP, uP*sol^{-1} ) ?>
+	    if (sol > 1.0) { sol =  1.0; }
+	} 
+
+<?R } else { ?>
+
+    if (sol > 1.0) { sol =  1.0; }
+    if (sol > 1.0e-8) {
+        <?R C( uP, uP*sol^{-1} ) ?>
+    }
+
+<?R } ?>
+
+    else {
 <?R
 	C(uP, 0)
 	C(sol, 0)

diff --git a/models/particles/d3q27_PSM/conf.mk b/models/particles/d3q27_PSM/conf.mk
@@ -1,3 +1,3 @@
 ADJOINT=0
 TEST=FALSE
-OPT="KL*TRT*(NEBB+SUP+(NEBB+SEP):singlekernel)"
+OPT="MS*KL*TRT*(NEBB+SUP+(NEBB+SEP):singlekernel)"
diff --git a/src/Consts.h.Rt b/src/Consts.h.Rt
@@ -24,21 +24,11 @@
 #define PART_MAR_BOX <?%f max(0,PartMargin-0.5) ?>
 
 <?R
-  big_hex = function(x,bits=16) {
-    ret = ""
-    while (x > 0 | bits > 0) {
-      a = x %% 16
-      ret = sprintf("%01x%s",a,ret)
-      x = (x-a) / 16
-      bits = bits - 4
-    }
-    ret
-  }
   for (n in rows(NodeTypes)) { ?>
-#define <?%-20s n$Index ?> 0x<?%s big_hex(n$value) ?> <?R
+#define <?%-20s n$Index ?> <?%s big_hex(n$value) ?> <?R
   }
   for (n in rows(NodeTypeGroups)) { ?>
-#define <?%-20s n$Index ?> 0x<?%s big_hex(n$mask) ?> <?R
+#define <?%-20s n$Index ?> <?%s big_hex(n$mask) ?> <?R
   }
 ?>
 <?R

diff --git a/src/Geometry.cpp b/src/Geometry.cpp
@@ -688,7 +688,7 @@ inline int Geometry::loadSweep(lbRegion reg, pugi::xml_node node)
 	int order=1;
 	attr = node.attribute("order");
 	if (attr) order = attr.as_int();
-	double dl=1e-3;
+	double dl=1e-4;
 	attr = node.attribute("step");
 	if (attr) dl = attr.as_double();
 	attr = node.attribute("steps");
@@ -719,19 +719,31 @@ inline int Geometry::loadSweep(lbRegion reg, pugi::xml_node node)
 		order = n-1;
 	}
     output("Making a Sweep over region (%ld) with %lf step\n",reg.size(), dl);
-    for (int x = reg.dx; x < reg.dx + reg.nx; x++)
-	for (int y = reg.dy; y < reg.dy + reg.ny; y++)
-	    for (int z = reg.dz; z < reg.dz + reg.nz; z++) {
-		for (double l=0;l<1; l+= dl) {
-			double x0 = bspline(l, nx, order);
-			double y0 = bspline(l, ny, order);
-			double z0 = bspline(l, nz, order);
-			double r = bspline(l, nr, order);
-			if ((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0) < r*r) {
-			    Dot(x, y, z);
+
+	double x0_=-999.0, y0_=-999.0, z0_=-999.0, r0_=-999.0;
+    for (double l=0;l<1; l+= dl) {
+		double x0 = bspline(l, nx, order);
+		double y0 = bspline(l, ny, order);
+		double z0 = bspline(l, nz, order);
+		double r = bspline(l, nr, order);
+
+		if (abs(x0-x0_) < 0.25 && 
+			abs(y0-y0_) < 0.25 && 
+			abs(z0-z0_) < 0.25 && 
+			abs(r-r0_) < 0.25) {
+				continue;		
+		}
+		x0_ = x0; y0_ = y0; z0_ = z0; r0_ = r;
+
+		lbRegion tmpReg=reg.intersect(lbRegion(x0-r-1,y0-r-1,z0-r-1,2*r+2,2*r+2,2*r+2));
+		for (int x = tmpReg.dx; x < tmpReg.dx + tmpReg.nx; x++)
+		for (int y = tmpReg.dy; y < tmpReg.dy + tmpReg.ny; y++)
+		for (int z = tmpReg.dz; z < tmpReg.dz + tmpReg.nz; z++) {
+			if ((.5+x-x0)*(.5+x-x0)+(.5+y-y0)*(.5+y-y0)+(.5+z-z0)*(.5+z-z0) < r*r) {
+				Dot(x, y, z);
 			}
 		}
-	    }
+	}
 
     return 0;	
 }