Corrected Darcy flow simulation

demos/Darcy/00_darcy_flow.py

@@ -6,6 +6,7 @@
 
 """
 Updecs on the Darcy-Flow equation with RBFs
+Challenge: The RBF is always continuous. Darcy's solution presents disocntinuities
 """
 
 import time
@@ -47,9 +48,11 @@ def perm_rhs_operator(x, centers=None, rbf=None, fields=None):
 cloud_perm = SquareCloud(Nx=Nx, Ny=Ny, facet_types=facet_types_perm, noise_key=key, support_size=SUPPORT_SIZE)
 
 def perm_func(x, y):
-    tmp = jnp.where(x>=y, 0., 0.90)
-    tmp = jnp.where(x<=0.5, tmp, 0)
-    return jnp.where(y**2>=0.1, tmp, 0.90)
+    tmp = jnp.where(y**2>=x, 0.1, 1.20)
+    tmp = jnp.where(0.25+y**1>=x, tmp, 0.1)
+    return 10*jnp.abs(tmp)
+    # tmp = jnp.where(x<=0.75, tmp, 0.2)
+    # return jnp.where(y**2<=x, tmp, 0.90)
 xy = cloud_perm.sorted_nodes
 permeability = jax.vmap(perm_func)(xy[:,0], xy[:,1])
 
@@ -62,8 +65,10 @@ def perm_func(x, y):
                     max_degree=MAX_DEGREE,)
 
 
-fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(3.75*2,3))
+fig, ax1 = plt.subplots(1, 1, figsize=(3.6, 3))
 cloud_perm.visualize_field(field=perm_field.vals, title="Permeability field", ax=ax1);
+plt.savefig(DATAFOLDER+"darcy_perm.png", dpi=300, bbox_inches="tight");
+
 
 
 # %%
@@ -76,21 +81,26 @@ def perm_func(x, y):
 # cloud.visualize_cloud(s=0.1, figsize=(4,3));
 
 def my_diff_operator(x, center=None, rbf=None, monomial=None, fields=None):
-    lap = nodal_laplacian(x, center, rbf, monomial)
-    return lap
+    perm_val = fields[0]
+    # perm_val = 1.
+    lap = nodal_div_grad(x, center, rbf, monomial, (perm_val, perm_val))
+    return -lap
+
 
 def my_rhs_operator(x, centers=None, rbf=None, fields=None):
-    return value(x, fields[:,0], centers, rbf)
+    # return value(x, fields[:,0], centers, rbf)
+    return 1.
 
 d_zero = lambda x: 0.
 boundary_conditions = {"South":d_zero, "West":d_zero, "North":d_zero, "East":d_zero}
 
 start = time.time()
 
 ufield = pde_solver_jit(diff_operator=my_diff_operator,
-                    rhs_operator = my_rhs_operator,
+                    rhs_operator=my_rhs_operator,
+                    diff_args=[perm_field.vals],
                     rhs_args=[perm_field.vals],
-                    cloud = cloud,
+                    cloud=cloud,
                     boundary_conditions = boundary_conditions, 
                     rbf=RBF,
                     max_degree=MAX_DEGREE,)
@@ -101,10 +111,11 @@ def my_rhs_operator(x, centers=None, rbf=None, fields=None):
 seconds = walltime % 60
 print(f"Walltime: {minutes} minutes {seconds:.2f} seconds")
 
-cloud.visualize_field(ufield.vals, cmap="jet", title=f"Solution field", ax=ax2, levels=200)
-plt.show()
+fig, ax2 = plt.subplots(1, 1, figsize=(3.65, 3))
+cloud.visualize_field(ufield.vals, cmap="jet", title=f"Solution field", ax=ax2, levels=200);
+# plt.show()
 
-# %%
+plt.draw()
+plt.savefig(DATAFOLDER+"darcy_sol.png", dpi=300, bbox_inches="tight");
 
-# plt.draw()
-plt.savefig(DATAFOLDER+"darcy_flow.png", dpi=300, bbox_inches="tight");
+# %%

docs/assets/NextRelease.md

@@ -1,2 +1,4 @@
 For the next release v1.1.0
-- [X] Added colorbar to animate fields
+- [X] Added colorbar to animate fields
+- [X] Fixed the args inputs to construct the local matrix for nodal_div_grad: (if array, else, etc.)
+- [X] Implemented the Darcy flow problem

updes/operators.py

@@ -92,11 +92,18 @@ def nodal_div_grad(x, center=None, rbf=None, monomial=None, args=None):
         center (Float[Array, "dim"]): The centroid of the RBF if used. (Currently mandadatory, despite the signature.)
         rbf (callable): The rbf to use. (Currently mandadatory, despite the signature.)
         monomial (callable): The monomial to use. (Currently mandadatory, despite the signature.)
+        args (list): The values to pre-multiply the gradients before taking the divergence
 
     Returns:
         float: The value of the laplacian rbf or polynomial at x
     """
-    matrix = jnp.stack((args, args), axis=-1)
+
+    ## Trick to duplicate the args into a matricx
+    if isinstance(args, jnp.ndarray):       ## a 2D array
+        matrix = jnp.stack((args, args), axis=-1)   
+    else:                                   ## a tuple
+        matrix = jnp.array((args, args))
+
     if center != None:
         return jnp.nan_to_num(jnp.trace(matrix * core_laplacian_rbf(rbf)(x, center)), posinf=0., neginf=0.)
     elif monomial != None: