How to input this PDE without errors?

[githubbuhtig10]

I wish to solve the following set of 2D coupled PDEs:

du/dt = -d/dx(alpha*v) + d2u/dx2 + c^2(1/y * du/dy + d2u/dy2)

dv/dt = au - b^2 * d/dx(alphau) + d2v/dx2 + c^2(1/y * dv/dy + d2v/dy2)

alpha is a function of x. a, b, and c are constants. 0<=x, y<=1 are the domains of x and y. Boundary conditions are homogeneous Dirichlet for x=1, y=0,1 and homogeneous Neumann for x=0 for both u and v. Initial conditions are u=1 and v=0.

I tried to modify the Brusselator example from the documentation to achieve this but kept facing errors. I have included the code I am trying to implement this with and kept commented out sections which cause errors.

I sometimes get an error saying dimensions don't match when using CartesianGrid rather than UnitGrid, specifically when trying to solve systems of equations.

Another error occurs when I try to change the boundary conditions to anything other than "natural," including manually typing what natural boundary conditions represent. The error states that laplace is not recognised.

The last error occurs when get_x/get_y is used and the message is : Untyped global name 'get_x': Cannot determine Numba type of <class 'function'> .

Any help to get this working would be appreciated.

The code I am trying is below:

from pde import PDE, FieldCollection, PlotTracker, ScalarField, UnitGrid, CartesianGrid

#Desired Boundary Conditions
#bc_x_left = {"derivative": [0, 0]}
#bc_x_right = {"value": [0, 0]}
#bc_y_bottom = {"value": [0, 0]}
#bc_y_top = {"value": [0, 0]}
#bc_x = [bc_x_left, bc_x_right]
#bc_y = [bc_y_bottom, bc_y_top]
#bc = [bc_x, bc_y]


# Wanted PDE:
#a, b, c = -3, 7, 1
#eq = PDE(
#    {
#        "u": f"-get_x(gradient(x * v)) + get_x(gradient(get_x(gradient(u)))) + {c} * {c} * ((1/y) * get_y(gradient(u)) + get_y(gradient(get_y(gradient(u)))))",
#        "v": f"{a} * u - {b} * {b} * get_x(gradient(x * u)) + get_x(gradient(get_x(gradient(v)))) + {c} * {c} * ((1/y) * get_y(gradient(v)) + get_y(gradient(get_y(gradient(v)))))",
#    },
#    bc="natural",
#    user_funcs={"get_x": lambda arr: arr[0], "get_y": lambda arr: arr[1]},
#)

#Test PDE
a, b, c = 1, 3, 1
eq = PDE(
    {
       "u": f"0.1*laplace(u)",
       "v": f"0.1*laplace(v)",
    },
    bc="natural",
)


#Test UnitGrid
grid = UnitGrid([64, 64])
u = a/b + 0.1 * ScalarField.random_normal(grid, label="Field $u$")
v = b / a + 0.1 * ScalarField.random_normal(grid, label="Field $v$")
state = FieldCollection([u, v])

#Wanted Cartesian Grid
#grid = CartesianGrid([[0, 1], [0, 1]], [32, 32])  # generate grid
#u = ScalarField(grid, 1, label="Field $u$")
#v = ScalarField(grid, 0, label="Field $v$")
#state = FieldCollection([u, v])

# Simulating PDE
sol = eq.solve(state, t_range=3)

#Plotting solution
sol.plot()

#Attempts to extract data and plot
#print(sol[1].slice({"x": 0.5}).data[:]) #prints data at x =30 across y at last time
#sol[1].plot()
#sol[1].slice({"x": 0.5}).plot(title='v against y at x = 30') #plots above

[david-zwicker]

I had a look at your problem and think I got it to work by using a different method for implementing the spatial derivative:

import functors
from pde import PDE, FieldCollection, PlotTracker, ScalarField, UnitGrid, CartesianGrid

#Desired Boundary Conditions
#bc_x_left = {"derivative": [0, 0]}
#bc_x_right = {"value": [0, 0]}
#bc_y_bottom = {"value": [0, 0]}
#bc_y_top = {"value": [0, 0]}
#bc_x = [bc_x_left, bc_x_right]
#bc_y = [bc_y_bottom, bc_y_top]
#bc = [bc_x, bc_y]

from pde.grids.operators.cartesian import _make_derivative
from pde.grids.cartesian import CartesianGridBase

make_gradient_x = functools.partial(_make_derivative, axis=0)
CartesianGridBase.register_operator('gradient_x', make_gradient_x)

make_gradient_y = functools.partial(_make_derivative, axis=1)
CartesianGridBase.register_operator('gradient_y', make_gradient_y)

# Wanted PDE:
a, b, c = -3, 7, 1
eq = PDE(
   {
       "u": f"-gradient_x(x * v) + gradient_x(gradient_x(u)) + {c} * {c} * (gradient_y(u) / y + gradient_y(gradient_y(u)))",
       "v": f"{a} * u - {b} * {b} * gradient_x(x * u) + gradient_x(gradient_x(v)) + {c} * {c} * (gradient_y(v) / y + gradient_y(gradient_y(v)))",
   },
   bc="natural"
)

#Wanted Cartesian Grid
grid = CartesianGrid([[0, 1], [0, 1]], [32, 32])  # generate grid
u = ScalarField.random_uniform(grid, 0.9, 1.1, label="Field $u$")
v = ScalarField(grid, 0, label="Field $v$")
state = FieldCollection([u, v])

# Simulating PDE
sol = eq.solve(state, t_range=3)

#Plotting solution
sol.plot()

#Attempts to extract data and plot
#print(sol[1].slice({"x": 0.5}).data[:]) #prints data at x =30 across y at last time
#sol[1].plot()
#sol[1].slice({"x": 0.5}).plot(title='v against y at x = 30') #plots above

Note that I actually register 1d-gradient operators for both Cartesian grids. To specify boundary conditions more finely grained, you would need to use the bc_ops argument of the PDE class.

[githubbuhtig10]

Thanks for your answer, I very much appreciate it. I noticed that I had to add 'import functools' for it to work when I copied it into the IDE but afterwards it worked flawlessly.

My only issue now is inputting the boundary conditions. My problem does not make sense with "natural" periodic boundary conditions unfortunately. You mentioned bc_ops in your answer and, in looking at the documentation, I am still uncertain of what I would input in order to overwrite "natural" in the bcs argument in order to attain my desired boundary conditions of homogeneous Dirichlet on three sides and homogeneous Neumann on the last. It seems bc_ops is about specifying boundary conditions for operators but I am unsure how I would relate this to the boundary conditions I want to input. In trying, I have found that anything other than "natural" in bcs causes operators in rhs to no longer be recognised. Would using bc_ops get around this? If so, do you know what I would input to achieve the boundary conditions I am looking for?

[david-zwicker]

I know that the interface for setting boundary conditions (BCs) is quite complex, but you need to remember that in principle different BCs can be specified for each differential operator and each condition can be different for each axis. To make life easier, the PDE class uses the following convention for specifying BCs via bc_ops, which is a dictionary, where they keys define the operators to wich the BC is applied and the value uses the normal convention of lists/dicts to specify different BCs on different axis. To allow specifying BCs for multiple operators, wildcards are supported. For instance, if you want to only specify the BCs for your u variable, you use the following setting: bc_ops={"u:*": [...]}, where [...] specifies the BCs for all axes. In case this does not work, please open another discussion with a minimal working example to separate it from this thread here.