Is it possible to solve a system of PDEs with Danckwerts boundary conditions?

[jAniceto]

Hi. I am trying to solve the following problem with two coupled PDEs.

$C$ and $q$ are the dependent variables and I am solving in 2 dimensions, time ($t$) and length ($z$).

Can this be solved with py-pde? The package looks more "physics oriented", which is not my field, so I am a bit overwhelmed with the examples/nomenclature. For instance, I am not sure it is possible to implement the Danckwerts boundary conditions.

Thanks in advance.

[david-zwicker]

Yes, such a simulation should be rather easy to do. The boundary condition looks like a Robin boundary condition to me; see documentation.

[jAniceto]

Hi @david-zwicker. Thanks for the quick reply.

I've been going through the docs and examples and they are of great quality.
I think I managed to make an implementation of the problem above.

# Spacial discretization: 1-D grid (axial position, z)
grid = CartesianGrid([[0, L]], 128, periodic=False)

# Initial conditions for C and q
C = ScalarField(grid, 0, label='C')   # t = 0, C = 0 for all z
q = ScalarField(grid, 0, label='q')   # t = 0, q = 0 for all z
state = FieldCollection([C, q])

# Boundary conditions
bc_z_left = {"type": "mixed", "value": f"-{v} / {D}", "const": f"-{v} * {C0} / {D}"}   # left, C - D/v * dCdz = C0
bc_z_right = {"derivative": 0}   # right, dCdz = 0
bc_z = [bc_z_left, bc_z_right]

# Define the coupled PDEs
eqs = PDE(
    {
        "C": f"{D} * laplace(C) - {v} * d_dx(C) - {F} * {K} * ({H} * C - q)",
        "q": f"{F} * {K} * ({H} * C - q)",
    }, 
    bc=bc_z
)

# Solve PDEs
storage = MemoryStorage()  # store intermediate information of the simulation
sol = eqs.solve(state, t_range=tfinal, dt=1e-3, tracker=["progress", "plot", storage.tracker(0.1)])

This runs without throwing any errors, however the results don't seem to make sense. A few questions:

1. Looks like integration stops at 0.006. I'm specifying above t_range=tfinal where tfinal is 7.

2. I'm not sure I am representing the second derivative in my first equation above correctly. I'm inputting laplace(C). What I want to input is $\frac{d^2C}{dz^2}$.

Apologies if these are silly questions.

[david-zwicker]

1. I don't really know why the simulation stops early. Could you provide a full working example? I could then try to investigate.
2. Using laplace is certainly fine, since you're doing a simulation in 1D, where laplace(C) = d^2C/dx^2. Alternatively, you could use d2_dx2.

[jAniceto]

Thanks.

I was misunderstanding the end time. It does integrate until tfinal. It looks mostly as expected.

The result I am more interested in is the $C$ over time at the final position $z$. I was extracting those values incorrectly.
Is there a quick way to extract this? I did manage to do it but in what seems a very convoluted way :p. This is probably more of a general Python problem though.

Thanks a lot for the very quick support. The package seems good for my use, I will continue to explore it.

[jAniceto]

Sorry to bother again. I got the case above working and I have been trying other stuff. Although, I'm having trouble making this one work.

I want C0 in my 1st boundary condition to depend on time so that C0 = 3 if t < 0.0125 else 0. I've tried two different approaches but none seem to work.

First, using heaviside, I changed the code above to:

bc_z_left = {"type": "mixed", "value": f"-{v} / {D}", "const": f"-{v} / {D} * {C0} * heaviside(0.0125-t, 0)"}

but I get the following error:

RuntimeError: ("Arguments {'t'} were not defined in expression signature ['x']", "Problems in boundary condition `[{'type': 'mixed', 'value': '-12.732395447351628 / 0.00557', 'const': '-12.732395447351628 / 0.00557 * 3.6 * heaviside(0.0125-t, 0)'}, {'derivative': 0}]` for operator `d2_dx2` in PDE for `C`")

Second, I tried defining the boundary condition with a function:

def bc_value(adjacent_value, dx, x, t):
    """define the boundary condition for a pulse"""
    if t <= tpulse:
        C0 = C0
    else:
        C0 = 0
    return -v / D * C0

bc_z_left = {"type": "mixed", "value": f"-{v} / {D}", "const": bc_value}

from where I get the following error:

TypeError: ("ufunc 'isnan' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''", "Problems in boundary condition `[{'type': 'mixed', 'value': '-12.732395447351628 / 0.00557', 'const': <function bc_value at 0x000001F544DA7880>}, {'derivative': 0}]` for operator `d_dx` in PDE for `C`")

Looking at section 3.3.1 and example 2.10 from the docs, it looks like this should be possible. What am I missing?

[david-zwicker]

Time-dependent Robin boundary conditions are not yet implemented. Please open an issue in the issue tracker and I'll see whether I can implement this anytime soon.

[jAniceto]

Ah, thanks. Just created issue #465.

[david-zwicker]

I just merged #464, which introduces mixed boundary conditions with expressions. You should be able to use the new boundary condition mixed_expression, e.g., bc_z_left = {"type": "mixed_expression", "value": f"-{v} / {D}", "const": bc_value}, where the expressions can now depend on time t.