Unexpected behavior from Laplacian operator in Cartesian

[Benj-Ward]

Hello. I am modeling a type of diffusion equation, and I have a spatially dependent concentration proportional to z, which has terms where it is operated on by a Laplacian or gradient operator. I would expect for the Laplacian operator that the term would go to zero, as a second derivative of z => 0; however, this is not the case and very different behavior occurs depending on whether the term is included or not (as shown below). The problem is defined with Dirichlet boundary conditions of 1 on all surfaces except for a subsection of the x=0,y=0 plane. I would appreciate guidance on the mathematical or coding issue going on here. Thanks.

This is the code with the term set to 0 explicitly.

#This is what we had before solved to steady state with the new format. I'm not sure if the software supports cutting off the
#left side of the graph.
import pde 
import numpy as np
import matplotlib.pylab as plt
from pde import CylindricalSymGrid, ScalarField, PDE, CartesianGrid, MemoryStorage, plot_kymograph

grid = CartesianGrid([[0, 32], [0, 32], [0,32]], [32,32,32]) # generate grid
state = ScalarField.random_uniform(grid, 0.2, 0.3) 
#implement an initial guess for initial scalar field. Use final state of above simulation.

trackers = [
    "progress",  # show progress bar during simulation
    "steady_state", #,  # abort when steady state is reached
    #storage.tracker(interval=10),  # store data every simulation time unit
    #pde.PlotTracker(show=True),  # show images during simulation
    #print(some output every 5 real seconds:
    #pde.PrintTracker(interval=pde.RealtimeInterrupts(duration=1))),
]

# Define Boundary conditions, value = dirichlet, derivative = neumann
CONST = np.zeros(32)
Ones = np.ones(32)
heterogeneous = np.array([])
for i in range(0,14):
    heterogeneous = np.append(heterogeneous,[1])
for i in range(0,4):
    heterogeneous = np.append(heterogeneous,[0])
for i in range(0,14):
    heterogeneous = np.append(heterogeneous,[1])

Dirichlet = np.vstack([Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,heterogeneous,
                   heterogeneous,heterogeneous,heterogeneous,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,
                   Ones,Ones,Ones,Ones])
#Neumann = value = np.vstack([CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,
                   #CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST])

    
bc_x_left = {"value":1}
bc_x_right = {"value":1}
bc_y_left = {"value":1}
bc_y_right = {"value":1}
bc_z_left = {"value":Dirichlet}
bc_z_right = {"value":1}

bc_x = [bc_x_left, bc_x_right]
bc_y = [bc_y_left, bc_y_right]
bc_z = [bc_z_left, bc_z_right]

cells = f"1/32 *z"
#term_1 = f"u * laplace({cells})"
term_1 = 0
term_2 = f"dot(gradient({cells}), gradient(u))"
term_3 = f"laplace(u)"

eq = pde.PDE({"u": f"{term_1} + {term_2} + {term_3}"}, bc= [bc_x,bc_y,bc_z])   # Define PDE 


result = eq.solve(state, 10**4, dt=0.01, tracker=trackers)
result.plot()

This is the code with the function (cells = z/32) being operated on by the Laplacian, which I would expect to be identical to the above result. In actuality, very different behavior occurs.

#This is what we had before solved to steady state with the new format. I'm not sure if the software supports cutting off the
#left side of the graph.
import pde 
import numpy as np
import matplotlib.pylab as plt
from pde import CylindricalSymGrid, ScalarField, PDE, CartesianGrid, MemoryStorage, plot_kymograph

grid = CartesianGrid([[0, 32], [0, 32], [0,32]], [32,32,32]) # generate grid
state = ScalarField.random_uniform(grid, 0.2, 0.3) 
#implement an initial guess for initial scalar field. Use final state of above simulation.

trackers = [
    "progress",  # show progress bar during simulation
    "steady_state", #,  # abort when steady state is reached
    #storage.tracker(interval=10),  # store data every simulation time unit
    #pde.PlotTracker(show=True),  # show images during simulation
    #print(some output every 5 real seconds:
    #pde.PrintTracker(interval=pde.RealtimeInterrupts(duration=1))),
]

# Define Boundary conditions, value = dirichlet, derivative = neumann
CONST = np.zeros(32)
Ones = np.ones(32)
heterogeneous = np.array([])
for i in range(0,14):
    heterogeneous = np.append(heterogeneous,[1])
for i in range(0,4):
    heterogeneous = np.append(heterogeneous,[0])
for i in range(0,14):
    heterogeneous = np.append(heterogeneous,[1])

Dirichlet = np.vstack([Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,heterogeneous,
                   heterogeneous,heterogeneous,heterogeneous,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,Ones,
                   Ones,Ones,Ones,Ones])
#Neumann = value = np.vstack([CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,
                   #CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST,CONST])

    
bc_x_left = {"value":1}
bc_x_right = {"value":1}
bc_y_left = {"value":1}
bc_y_right = {"value":1}
bc_z_left = {"value":Dirichlet}
bc_z_right = {"value":1}

bc_x = [bc_x_left, bc_x_right]
bc_y = [bc_y_left, bc_y_right]
bc_z = [bc_z_left, bc_z_right]

cells = f"1/32 *z"
term_1 = f"u * laplace({cells})"
#term_1 = 0
term_2 = f"dot(gradient({cells}), gradient(u))"
term_3 = f"laplace(u)"

eq = pde.PDE({"u": f"{term_1} + {term_2} + {term_3}"}, bc= [bc_x,bc_y,bc_z])   # Define PDE 


result = eq.solve(state, 10**4, dt=0.01, tracker=trackers)
result.plot()

[david-zwicker]

I can reproduce the behavior and I think it is expected (and thus not a bug). You are correct that the Laplace operator applied to a linearly varying field should yield zero. However, at the boundary different values emerge since the boundary condition is applied there, too. This can be seen from the following evaluation of the laplacian alone:

from pde.tools.expressions import evaluate
res = evaluate(f"laplace({cells})", {"u": state}, bc=[bc_x, bc_y, bc_z])
res.plot()

In the plot, you can see that the term vanishes everywhere, except at the boundaries, where the boundary conditions result in non-zero values.

You can define a special operator with different boundary conditions using the following trick:

res2 = evaluate(
    f"laplace_neumann({cells})",
    {"u": state},
    bc=[bc_x, bc_y, bc_z],
    user_funcs={"laplace_neumann": grid.make_operator("laplace", "neumann")},
)
res2.plot()

[Benj-Ward]

Thanks for this response. I followed it quite a bit in trying to answer my questions.