wip porous

example/porous_media/fromPoro2Perm.py

@@ -2,60 +2,116 @@
 import matplotlib.pyplot as plt
 import cv2
 
-image = cv2.imread('porous_media.png') 
+
+# %% 
+path= '/home/grzegorz/GITHUB/LBM/TCLB/example/porous_media/'
+image = cv2.imread(path+'porous_media_contrast_1024x256.png') 
 # image = cv2.imread('japan_1024x640.png') 
 
 gray_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # convert to grayscale
-np_gray_image = np.array(gray_image)
+# gray_image = image
+# np_gray_image = np.array(gray_image)
+# Convert to float
+float_image = gray_image.astype(np.float32)
+
+# Rescale to 0-1 range
+normalized_image = float_image / 255.0
+
 shape = gray_image.shape
+print(f"input shape: {shape}")
 
 
+# %% 
 # Parameters for the Kozeny-Carman equation
-K = 200  # Kozeny constant, arbitrary choice for demonstration
-S = 1  # Specific surface area, arbitrary choice for demonstration
+K = 1.  # Kozeny constant, arbitrary choice for demonstration
+S = 1.  # Specific surface area, arbitrary choice for demonstration
 
 # Epsilon values to prevent division by zero in the denominator
-epsilon = 1E-3
+epsilon = 1E-2
 
 # Porosity range from 0 to 1 (not including 1 to avoid division by zero in the original equation)
-phi = np_gray_image #np.linspace(0.0, 0.9, 200)
+phi = normalized_image #np.linspace(0.0, 0.9, 200)
 
 # Prepare the plot
-plt.figure(figsize=(10, 7))
+
 
 # Calculate and plot permeability for each epsilon value
 
+def make_plot(x, title, cmap):
+    plt.figure(figsize=(10, 7))
+    plt.imshow(x, cmap=cmap)
+    plt.colorbar()
+    plt.title(f'{title}')
+    plt.show()
+    plt.close()
+
+def save_img(x, h, w, cmap, title='generated_permability.png'):
+
+#     plt.figure(figsize=(3.841, 7.195), dpi=100)
+# ( your code ...)
+# plt.savefig('myfig.png', dpi=1000)
+    # plt.figure(figsize=(5.12, 1.024), dpi=100)
+    # plt.axis('off')
+    # plt.imshow(x, cmap='gist_gray') #  vmax=1000
+
+    # plt.savefig(f'{title}', bbox_inches='tight', pad_inches=0, dpi=100)
+    # plt.close()
+
+    # w = 1000
+    # h = 512
+    # im_np = numpy.random.rand(h, w)
+
+
+    my_dpi = 100
+    fig = plt.figure(figsize=(w/my_dpi, h/my_dpi), dpi=my_dpi)
+
+    ax = plt.Axes(fig, [0., 0., 1., 1.])
+    ax.imshow(x, cmap=cmap)
+    fig.add_axes(ax)
+    # plt.imshow(x, cmap='gist_gray')
+
+    plt.axis('off')
+    plt.savefig(f'{title}', bbox_inches='tight', pad_inches=0, dpi=my_dpi)
+
+
+
+def rescale(background, new_min, new_max):
+    # new_min = 0.01 +  1* np.min(background) # Define new minimum value
+    # new_max = 0.99 * np.max(background)  # Define new maximum value
+    min  = np.min(background)
+    max = np.max(background)
+    new_backgroud = new_min + ((background - min) * (new_max - new_min)) / (max - min)
+    return new_backgroud
+
+
+# Define the smooth limiter function based on the given C code
+def smooth_limiter(x, max=10, steepness = 0.1):
+    result = x*(3*x + 1)/(x+1)**2
+    # result = max / (1 + np.exp(-steepness * (x - max)))
+    return result
+
 
     # Adjusted Kozeny-Carman equation to include epsilon in the denominator
 def Kozeny_Carman(phi,K,S):
     return (phi**3) / (K * ((1 - phi)**2 + epsilon)) * (1/S**2)
 
-k = Kozeny_Carman(phi,K,S)
-
-
-# Define the smooth limiter function based on the given C code
-def smooth_limiter(x):
-    limit = 1000.0
-    steepness = 0.1  # Adjust steepness for a smoother or sharper transition
-    result = limit / (1 + np.exp(-steepness * (x - limit)))
-    return result
+rephi = rescale(phi, 0.01, 0.99)
+k = Kozeny_Carman(rephi,K,S)
 
+make_plot(rephi, 'Input Porosity', 'gist_gray')
+make_plot(k, 'Permability', 'gist_gray')
 
-# plt.imshow(gray_image, cmap='gist_gray')
-# plt.colorbar()
-# plt.title(f'Porosity')
-# plt.show()
+# ks = rescale(k, 0.01, 1000) 
+# make_plot(ks, 'RescaledPermability')
+# h, w, c = k.shape
+h, w = k.shape
+save_img(k, h, w, 'gist_gray')
 
-# plt.imshow(smooth_limiter(k), cmap='gist_gray')
-# plt.colorbar()
-# plt.title(f'Clipped Permability')
-# plt.show()
+# %%
 
-plt.axis('off')
-plt.imshow(k, cmap='gist_gray')
-# plt.colorbar()
-# plt.title(f'Permability')
-# plt.show()
+# x = np.linspace(0,10,1000)
+# y = smooth_limiter(x, max=10, steepness= 1)
 
-plt.savefig('generated_permability.png', bbox_inches='tight', pad_inches=0, dpi=100)
-plt.close()
+# plt.plot(x,y)
+# plt.plot(x,x)
+# plt.grid()

models/porous_media/d2q9_porous_heat/Dynamics.R

@@ -90,6 +90,7 @@ AddQuantity(name="H", 		unit="J" )
 AddQuantity(name="T", 		unit="K")
 AddQuantity(name="Porosity", unit="1")
 AddQuantity(name="NSPermability", unit="1")
+AddQuantity(name="HPermability", unit="1")
 
 AddQuantity(name="Hflux", unit="1", ,vector=T)
 AddQuantity(name="HDflux", unit="1", ,vector=T)

models/porous_media/d2q9_porous_heat/Dynamics.c.Rt

@@ -575,6 +575,8 @@ CudaDeviceFunction void relax_and_collide_ADE_CM_HIGHER(real_t rho, real_t omega
 	h120 = temp110 + temp120 - temp210 - temp220;
 	h220 = temp110 + temp120 + temp210 + temp220;
 
+	// h100 += fluxHDarcy.x/2.;
+	// h010 += fluxHDarcy.y/2.;
 
 	//central moments from raw moments
 	temp000 = h000;
@@ -789,8 +791,8 @@ CudaDeviceFunction void CollisionCM_HIGHER(){
 	vector_t fluxHDarcy = getHDarcyForce(Hflux, permability);
 	// fluxHDarcy.x = 0;
 	// fluxHDarcy.y = 0;
-	real_t omega_k = GetOmegaTRTPlus(conductivity, permability);
-	// real_t omega_k = 1.0/(3*conductivity+0.5);
+	// real_t omega_k = GetOmegaTRTPlus(conductivity, permability);
+	real_t omega_k = 1.0/(3*conductivity/rho +0.5);
 
 
 

models/porous_media/d2q9_porous_heat/ForcingTerms.c.Rt

@@ -18,7 +18,7 @@ CudaDeviceFunction real_t getPermability(real_t Kazeny_coeff)
 	real_t ssa = 1. ; // the specific surface area (surface area per unit volume of solid).
 	real_t permability;
 	#ifdef OPTIONS_PORO121PERM
-		permability = poro0;
+		permability = poro0/Kazeny_coeff;
 	#else
 		real_t eps = 1E-3;
 		permability = poro0*poro0*poro0; 

models/porous_media/d2q9_porous_heat/VTKOutputs.c.Rt

@@ -28,18 +28,42 @@ CudaDeviceFunction real_t getNSPermability()
 	// The Kozeny-Carman equation provides a way to estimate the permeability of a porous medium based on its porosity and specific surface area.
 	real_t ssa = 1. ; // the specific surface area (surface area per unit volume of solid).
 	real_t permability;
-	#ifdef OPTIONS_PORO121PERM
-		permability = poro0;
-	#else
-		permability = poro0*poro0*poro0; 
-		permability /= (Kazeny_NScoeff*(1-poro0)*(1-poro0)+1E-3);
-		permability /= (ssa*ssa);
-	#endif
+
+	permability = getPermability(Kazeny_NScoeff);
+
+	// #ifdef OPTIONS_PORO121PERM
+	// 	permability = poro0/Kazeny_NScoeff;
+	// #else
+
+	// 	permability = poro0*poro0*poro0; 
+	// 	permability /= (Kazeny_NScoeff*(1-poro0)*(1-poro0)+1E-3);
+	// 	permability /= (ssa*ssa);
+	// #endif
 
 	return permability;
 }
 
 
+CudaDeviceFunction real_t getHPermability()
+{	
+	// Use this function is only for vtk output.
+
+	// Darcy Forcing
+	// The Kozeny-Carman equation provides a way to estimate the permeability of a porous medium based on its porosity and specific surface area.
+	real_t ssa = 1. ; // the specific surface area (surface area per unit volume of solid).
+	real_t permability;
+	permability = getPermability(Kazeny_Hcoeff);
+	// #ifdef OPTIONS_PORO121PERM
+	// 	permability = poro0/Kazeny_Hcoeff;
+	// #else
+	// 	permability = poro0*poro0*poro0; 
+	// 	permability /= (Kazeny_Hcoeff*(1-poro0)*(1-poro0)+1E-3);
+	// 	permability /= (ssa*ssa);
+	// #endif
+
+	return permability;
+}
+
 
 CudaDeviceFunction real_t getRho(){
 	// Use this function is only for vtk output.