Add basic physical scaling and Knudsen quality criterion

The paper "Automatic grid refinement criterion for lattice Boltzmann method" [2015]
by Lagrava et al. describes a criterion for measuring the local simulation quality
using a comparison of the theoretical Knudsen number and the quotient of the cells's
non-equilibrium and equilibrium function.

While this criterion was developed to enable automatic selection of areas to be refined,
it also offers a interesting and unique perspective on the fluid structure.

As the criterion requires calculation of the modeled Reynolds-, Mach- and Knudsen-numbers
I took the time to set up the basics for scaling the simulation to actually model a physical
system. Or rather calculating which physical model is represented by the chosen resolution
and relaxation time.

[2015]: https://arxiv.org/abs/1507.06767

2 files changed, 57 insertions, 10 deletions

diff --git a/src/shader/code/collide.glsl b/src/shader/code/collide.glsl
index 3bfd332..9ff6ab3 100644
--- a/src/shader/code/collide.glsl
+++ b/src/shader/code/collide.glsl
@@ -11,6 +11,12 @@ uniform uint nX;
 uniform uint nY;
 uniform uint iT;
 
+/// Fluid characteristics
+
+const float physCharLength   = 1.0;
+const float physCharVelocity = 1.0;
+const float physViscosity    = 0.1;
+
 /// LBM constants
 
 const uint  q         = 9;
@@ -19,9 +25,25 @@ const float weight[q] = float[](
 	1./9. , 4./9., 1./9. ,
 	1./36 , 1./9., 1./36.
 );
+const float invCs2 = 1./3.;
+
+const float relaxationTime      = 0.6;
+const float relaxationFrequency = 1 / relaxationTime;
+
+/// Unit conversion
+
+const float convLength    = physCharLength / nX;
+const float convTime      = (relaxationTime - 0.5) / invCs2 * convLength*convLength / physViscosity;
+const float convVelocity  = convLength / convTime;
+const float convViscosity = convLength * convLength / convTime;
+
+const float latticeCharVelocity = physCharVelocity / convVelocity;
+
+/// Emergent fluid numbers
 
-const float tau   = 0.65;
-const float omega = 1/tau;
+const float Re = physCharVelocity * physCharLength / physViscosity;
+const float Ma = latticeCharVelocity * sqrt(invCs2);
+const float Kn = Ma / Re;
 
 /// Vector utilities
 
@@ -65,12 +87,18 @@ void set(uint x, uint y, int i, int j, float v) {
 	streamCells[indexOfLatticeCell(x,y) + indexOfDirection(i,j)] = v;
 }
 
-void setFluid(uint x, uint y, vec2 v) {
+void setFluidVelocity(uint x, uint y, vec2 v) {
 	const uint idx = indexOfFluidVertex(x, y);
 	fluidCells[idx + 0] = v.x;
 	fluidCells[idx + 1] = v.y;
 }
 
+void setFluidQuality(uint x, uint y, float knudsen, int quality) {
+	const uint idx = indexOfFluidVertex(x, y);
+	fluidCells[idx + 0] = knudsen;
+	fluidCells[idx + 1] = quality;
+}
+
 int getMaterial(uint x, uint y) {
 	const uint idx = indexOfFluidVertex(x, y);
 	return int(fluidCells[idx + 2]);
@@ -149,19 +177,29 @@ void main() {
 			d = 1.0;
 		}
 
-		setFluid(x,y,v);
+		//setFluidVelocity(x,y,v);
+
+		float knudsen = 0.0;
 
 		for ( int i = -1; i <= 1; ++i ) {
 			for ( int j = -1; j <= 1; ++j ) {
-				set(x+i,y+j,i,j, get(x,y,i,j) + omega * (equilibrium(d,v,i,j) - get(x,y,i,j)));
+				const float feq  = equilibrium(d,v,i,j);
+				const float fneq = get(x,y,i,j) - feq;
+				knudsen += abs(fneq / feq);
+
+				set(x+i,y+j,i,j, get(x,y,i,j) + relaxationFrequency * (feq - get(x,y,i,j)));
 			}
 		}
+
+		knudsen /= q;
+
+		setFluidQuality(x,y,knudsen,int(round(log2(knudsen / Kn))));
 	}
 
 	if ( isBounceBackCell(material) ) {
 		for ( int i = -1; i <= 1; ++i ) {
 			for ( int j = -1; j <= 1; ++j ) {
-				set(x+(-1)*i,y+(-1)*j,(-1)*i,(-1)*j, get(x,y,i,j) + omega * (equilibrium(d,v,i,j) - get(x,y,i,j)));
+				set(x+(-1)*i,y+(-1)*j,(-1)*i,(-1)*j, get(x,y,i,j) + relaxationFrequency * (equilibrium(d,v,i,j) - get(x,y,i,j)));
 			}
 		}
 	}
diff --git a/src/shader/code/vertex.glsl b/src/shader/code/vertex.glsl
index 676b68e..7851292 100644
--- a/src/shader/code/vertex.glsl
+++ b/src/shader/code/vertex.glsl
@@ -10,7 +10,8 @@ out VS_OUT {
 uniform uint nX;
 uniform uint nY;
 
-const float displayAmplifier = 3.0;
+const float velocityDisplayAmplifier = 3.0;
+const int   qualityDisplayRestrictor = 6;
 
 float unit(float x) {
 	return 1.0/(1.0+exp(-x));
@@ -40,6 +41,14 @@ bool isWallFrontier(int material) {
 	return material == 2 || material == 3;
 }
 
+float restrictedQuality(float quality) {
+	if ( quality < 0.0 ) {
+		return 0.0;
+	} else {
+		return min(1.0, quality / qualityDisplayRestrictor);
+	}
+}
+
 void main() {
 	const vec2 idx = fluidVertexAtIndex(gl_VertexID);
 
@@ -58,9 +67,9 @@ void main() {
 		vs_out.color = vec3(0.0, 0.0, 0.0);
 	} else {
 		vs_out.color = mix(
-			vec3(-0.5, 0.0, 1.0),
-			vec3( 1.0, 0.0, 0.0),
-			displayAmplifier * norm(VertexPosition.xy)
+			vec3(0.0, 1.0, 0.0),
+			vec3(1.0, 0.0, 0.0),
+			restrictedQuality(VertexPosition.y)
 		);
 	}
 }