zwicker-group / py-pde

Python package for solving partial differential equations using finite differences.
https://py-pde.readthedocs.io
MIT License
410 stars 52 forks source link

Falls back to basic python backend if user-supplied function cannot be compiled using numba. #498

Closed Waschenbacher closed 9 months ago

Waschenbacher commented 9 months ago

Discussed in https://github.com/zwicker-group/py-pde/discussions/488

Originally posted by **Waschenbacher** November 30, 2023 I am implementing a custom PDE class and would like to add a not-so-easy boundary condition. Is it possible in py-PDE to add a Dirichlet boundary condition with its value being a function of - its fields.data at coordinate (x = L) - f(t) Example: I have a field A on 1-d "grid = CartesianGrid([[0, 10]], [100])]". The boundary condition at left-hand-side should be a function of the field at the right and a function of time t.
Waschenbacher commented 9 months ago

Thank you for the quick feedback. I tried it out but it gave me similar error message when running the example below:

`from pde import ( CartesianGrid, MemoryStorage, plot_kymograph, DiffusionPDE, ScalarField, )

time = 0 integral = 0 time_prev = -1e-6 e_prev = 0

def partial(func, *args, keywords): def newfunc(*fargs, *fkeywords): newkeywords = keywords.copy() newkeywords.update(fkeywords) return func((args + fargs), newkeywords) newfunc.func = func newfunc.args = args newfunc.keywords = keywords return newfunc

def PID(Kp, Ki, setpoint, measurement): global time, integral, time_prev, e_prev

Value of offset - when the error is equal zero

offset = 1.
# PID calculations
e = setpoint - measurement
P = Kp*e
integral = integral + Ki*e*(time - time_prev)
MV = offset + P + integral
# update stored data for next iteration
e_prev = e
time_prev = time

# only output non-negative value
if MV <= 0:
    return 0.
else:
    return MV

user_func = partial(PID, 1e-3, 0., 1.) print(user_func(1.))

bc_closed_loop_left = { "type": "value_expression", "value": "func(value)", "user_funcs": {"func": user_func}, "value_cell": -1, } bc_right = {"derivative": 0} bc = (bc_closed_loop_left, bc_right)

grid = CartesianGrid([[0, 10]], [100]) # generate grid state = ScalarField.random_uniform(grid, 0.2, 0.3)

eq_diffusion_PDE = DiffusionPDE(diffusivity=1e-1, bc=bc) storage = MemoryStorage() sol = eq_diffusion_PDE.solve(state, t_range=20, dt=1e-4, tracker=storage.tracker(0.01)) plot_kymograph(storage)`

Error message is attached as well

Screenshot 2023-12-17 at 11 49 00
david-zwicker commented 9 months ago

Could you please format the code in your last post correctly, so it is easier to read? It's hard to read python code that is not indented correctly :)

Waschenbacher commented 9 months ago

Sorry for that. I tried multiple times with copying my codes using "Code" but it did not work out. I attached the code in the txt file below. debug_py_PDE.txt

david-zwicker commented 9 months ago

I copy the code here for easier reference

from pde import (
    CartesianGrid,
    MemoryStorage,
    plot_kymograph,
    DiffusionPDE,
    ScalarField,
)

def partial(func, *args, **keywords):
    def newfunc(*fargs, **fkeywords):
        newkeywords = keywords.copy()
        newkeywords.update(fkeywords)
        return func(*(args + fargs), **newkeywords)
    newfunc.func = func
    newfunc.args = args
    newfunc.keywords = keywords
    return newfunc

time = 0
integral = 0
time_prev = -1e-6
e_prev = 0

def PID(Kp, Ki, setpoint, measurement):
    global time, integral, time_prev, e_prev
    # Value of offset - when the error is equal zero
    offset = 1.
    # PID calculations
    e = setpoint - measurement
    P = Kp*e
    integral = integral + Ki*e*(time - time_prev)
    MV = offset + P + integral
    # update stored data for next iteration
    e_prev = e
    time_prev = time

    # only output non-negative value
    if MV <= 0:
        return 0.
    else:
        return MV

# # user defined boundary: related to sin func
# def user_func(x: float) -> float:
#     """Time-dependent PDE parameter"""
#     return (1 + np.sin(x)) * 0.5

# user-defined boundary condition
user_func = partial(PID, 1e-3, 0., 1.)
print(user_func(1.))

# boundary condition
bc_closed_loop_left = {
    "type": "value_expression",
    "value": "func(value)",
    "user_funcs": {"func": user_func},
    "value_cell": -1,
}
bc_right = {"derivative": 0}
bc = (bc_closed_loop_left, bc_right)

# grid and state
grid = CartesianGrid([[0, 10]], [100])  # generate grid
state = ScalarField.random_uniform(grid, 0.2, 0.3)

# simulate the pde
eq_diffusion_PDE = DiffusionPDE(diffusivity=1e-1, bc=bc)
storage = MemoryStorage()
sol = eq_diffusion_PDE.solve(state, t_range=20, dt=1e-4, tracker=storage.tracker(0.01))
plot_kymograph(storage)
david-zwicker commented 9 months ago

I can reproduce your problem, but I failed to fix the problem quickly. I will continue to investigate, but it might take a while...