Dear Dr. Zwicker,

I am really appreciated that you developed such an amazing package for us to solve PDE. I am really interested in the example, called Brusselator - Using the PDE class. But there are conditions there which I am still confused about. Could you help me to figure out them?

Firstly, does t_range=20 mean t $\in$ [0,20] or mean t$\in$ [0, 20dt]?

Secondly, for bc="natural", does it mean there is no bounday condition?

Thirdly, If I set up bc="neumann", does it means neumann condition is 0?

Thank you so much!

Applying the transpose() methods to a Tensor2Field as

grid = pde.SphericalSymGrid(radius=[0, 5], shape=[100]) # generate grid state = pde.Tensor2Field(grid, 0.1, label='Composition') state.transpose()

returns the error:

ValueError Traceback (most recent call last)
in
1 grid = pde.SphericalSymGrid(radius=[0, 5], shape=[100]) # generate grid
2 state = pde.Tensor2Field(grid, 0.1, label='Composition')
----> 3 state.transpose()
~/Repos/py-pde/pde/fields/tensorial.py in transpose(self, label)
362 """
363 axes = (1, 0) + tuple(range(2, 2 + self.grid.dim))
--> 364 return Tensor2Field(self.grid, self.data.transpose(axes), label=label)
365
366 def symmetrize(
ValueError: axes don't match array

First of all, thanks for sharing this package, it's awesome!

I am trying to produce movies that include stochastic noise. Fluctuations in the field do produce a lot of changes that due to video compression result in very poor quality videos.

Therefore I would like to increase the dpi and bitrate of the video in order to make the images more clear. In the docs I can see that the base Movie class do have options to control these, but apparently, I cannot feed bitrate as a parameter to other functions, such as movie. Maybe this can be accessed from plot_args but I don't know how to do it.

Documentation is not very clear about how to do this. Some things I have tried:

1. Using a MemoryStorage as written in docs, section 3. Result: produces an mp4 file and plays, but I cannot adjust the bitrate.

    storage = MemoryStorage()
    tracker = PlotTracker(interval=1.0, movie=mymovie, plot_args={"vmin":0, "vmax": 0.5, "figsize":(10,5)})
    sol = eq.solve(state, t_range=10.0, dt=1e-2, tracker=["progress", storage.tracker(1)])
    mymovie = movies.movie(storage, "pelicula.mp4", dpi=300, plot_args={"vmin":0, "vmax": 0.5})

2. Try to feed a Movie to a PlotTracker. Result: writes a file, but it cannot be read (I guess because of lack of MemoryStorage)

    mymovie = movies.Movie("pelicula.mp4", dpi=300, bitrate=100)
    tracker = PlotTracker(interval=1.0, movie=mymovie, plot_args={"vmin":0, "vmax": 0.5, "figsize":(10,5)})
    sol = eq.solve(state, t_range=10.0, dt=1e-2, tracker=tracker)

I guess there is an easy way to control movie parameters, but at least for me, it is not clear how to do this. I would appreciate some extra indications on how to increase bitrate and framerate. Thank you very much in advance!

You would better specify the versions of python and anaconda/conda for the codes uploaded. The programs do not seem to be robust. Your tutorial2 seems to have a more serious library issue.

The following error message is typical:

solver = pde.ExplicitSolver(eq)
controller = pde.Controller(solver, t_range=1)
sol1 = controller.run(field, dt=1e-3);
/srv/conda/envs/notebook/lib/python3.7/site-packages/pde/tools/output.py:61: UserWarning: tqdm package is not available. Progress will be indicated by dots.
"tqdm package is not available. Progress will be indicated by dots."

In the example programs of 'Different solvers.ipynb' & 'Solve PDEs.ipynb', this error message is showing up but not in 'Discretized Fields.ipynb'. Several similar messages are showing up in tutorial2.

[openjournals/joss-reviews#2158]

The example gallery is great, you might consider including the notebooks in examples/jupyter in the documentation as well. With the sphinx-nbexamples extension, this should work very similar to the sphinx-gallery setup.

Hi all,

This package seems great! I've been looking for something like this for a while so I'm happy to try it out.

I want to simulate the PDE pictured below.

The sigma field is just some source term that I specify manually. (added equation for sigma so PDE class would recognize variable)

When I integrate this equation

I find that the left boundary condition on the field u is not satisfied (didn't check derivative BC on right side yet). I want u(0) = 0 but when I inspect the solution fields I find that u(0) is not actually zero.

It seems that the BCs have been ignored or I haven't specified them properly. I tried making a tuple for each field
bc = [({"value":0},{"derivative": -0.047}), ({"value":0},{"value":10})] (first tuple has u BCs and second tuple has sigma BCs) but when I do that I get the following error

I've looked through the docs for info on specifying BCs for both fields in a coupled PDE but I wasn't able to find any. Do you know how I can solve this issue?

As part of openjournals/joss-reviews#2158 I found two minor issues with building the sphinx documentation:

• The requirements stated on https://py-pde.readthedocs.io/en/latest/getting_started.html#prerequisites do not include sphinx-gallery.
• The docs/requirements.txt file does not include pillow, but it seems to be necessary to generate thumbnails for the sphinx gallery.

Thanks a lot for this!
I can't seem to find an explanation on how to set up my own boundary conditions as an array though...
I managed to do it using:
bc = {"type" : "value", "value" : <the boundary list>}
though it took some faith in the creators. Maybe it could be added to the documentation.
Thanks!

Thanks for building this awesome package.

I have used the example SIR model as a starting point for my work.
From this example I got to understand how the boundary conditions are applied to the laplace operator: e.g. s.laplace(self.bc).
But that is a boundary condition for the laplace operator applied to the field, how do you specify a boundary condition for the field itself in the framework of the SIR example?

And how do you do that when not all fields have the same boundary conditions, contrary to the SIR example?

From this discussion I understand that bc_ops can used for different boundary conditions for different fields. But that discussion does not use the custom PDE class approach as in the example SIR model, so I am unable to map that to something I could use.

Is there a method for this package to extract the components of the gradient when constructing equations with PDE class, such as eq = PDE({'u':"-gradient_x(u))"})? I did not find a direct way to implement it in the document. Can you add a similar function?

or apply some functions to the point in the domain?

I am having trouble with setting the boundary conditions of the wave equation so that the wave does not reflect off of the edges of the grid. We tried playing around with various configurations of bc_x and bc_y which are all hashed in the code below, like setting the boundaries to be equal to zero, but nothing worked.

import numpy as np
from pde import FieldCollection
from pde import (DiffusionPDE, UnitGrid, ScalarField, MemoryStorage,
                 PlotTracker, PrintTracker, RealtimeIntervals, DataTracker,
                 ExplicitSolver, ScipySolver, Controller)
from pde.visualization import movies
import math
import random
import pde
import pandas
import matplotlib.pyplot as plt; plt.rcdefaults()
import numpy as np
import matplotlib.pyplot as plt


storage = MemoryStorage()

trackers = [
    'progress',                    # show progress bar during simulation
    'steady_state',                # abort when steady state is reached
    storage.tracker(interval=1),   # store data every simulation time unit
    PlotTracker(show=True),        # show images during simulation
]

grid = pde.UnitGrid([64, 64])
field1 = ScalarField.random_uniform(grid, 0, 0)  # generate initial condition
field2 = ScalarField(grid);
field1.add_interpolated([48, 48], 1)
field1.add_interpolated([48,16], 1)
state = FieldCollection([field1,field2]); #first slot is the density field,

#bc_x_left = {'derivative': 0.1}
#bc_x_right = {'value': 'sin(y)'}
#bc_x = [bc_x_left, bc_x_right]
##bc_y = 'periodic'
#bc_y = bc_x
#bc_x = ({'value': 0}, {'derivative': 0})
#bc_y = ({'value': 0}, {'derivative': 0})
#bc_x = ({'value': 0})
#bc_y = ({'value': 0})
#bc_x = ({'derivative' : 0 })
#bc_y = ({'derivative' : 0})
#bc_x = ('dirichlet', {'value' : 0 })
#bc_y = ('dirichlet', {'value' : 0 })
#bc_x = ('dirichlet', {'value' : 0 })
#bc_x = ({'type' : 'dirichlet', 'value' : 0})
bc_y = ({'type': 'mixed', 'value': 0, 'const': 0})
bc_x = ({'type': 'mixed', 'value': 0, 'const': 0})


#eq = pde.WavePDE(speed=1,bc=[[bc_x, bc_y],[bc_x,bc_y]]);
eq = pde.WavePDE(speed=1,bc=[bc_x, bc_y])
eq.check_implementation = False
result = eq.solve(state, t_range=50, dt=0.0001, tracker = trackers, method="implicit")
#solver1 = ScipySolver(eq)
#controller1 = Controller(solver1, t_range=50, tracker = trackers)
#result = controller1.run(state, dt=1e-4);
#result = eq.solve(state, t_range=50, dt=0.0001, tracker = trackers, method="explicit")
result.plot(colorbar = True, show = True)       # plot the resulting field
movies.movie_multiple(storage, './test.mp4')

[openjournals/joss-reviews#2158]

When I run the example jupyter notebooks, the tracker displaying an error bar with tqdm fails if ipywidgets is not installed.

I'm trying to simulate a square drum, i.e. membrane vibrations by evolving the 2D wave equation.
I get the following error if I try to set the boundary to zero

AssertionError:
Not equal to tolerance rtol=1e-07, atol=1e-07
The numba compiled implementation of the right hand side is not compatible with the numpy implementation. This check can be disabled by setting the class attribute `check_implementation` to `False`.
Mismatched elements: 252 / 8192 (3.08%)
Max absolute difference: 0.02450429
Max relative difference: 1.0097242
 x: array([[[ 0.000000e+00,  0.000000e+00,  0.000000e+00, ...,
          0.000000e+00,  0.000000e+00,  0.000000e+00],
        [ 0.000000e+00,  0.000000e+00,  0.000000e+00, ...,...
 y: array([[[ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.      ,
          0.      ],
        [ 0.      ,  0.      ,  0.      , ...,  0.      ,  0.      ,...

If I explicitly set the check_implementation to False I can solve the equation, but what does that mean, can I trust that solution? If I remove the boundary conditions I can solve successfully. Maybe I'm doing some mistake in my implementation, I'm not entirely sure how to interpret the boundary conditions when there are multiple fields. Can someone help? Below is my code.

"""
Simulate waves on membrane fixed to square frame.
"""

import pde
from pde import WavePDE, ScalarField, UnitGrid, FieldCollection

# Square shaped membrane
grid = UnitGrid([64, 64], periodic=[False, False])
u = ScalarField.from_expression(grid, "sin(2*pi*x/64)*sin(2*pi*y/64)", label="Field $u$")
v = ScalarField(grid, label="Field $v$") 
state = FieldCollection([u, v])

# define the pde, fixed frame: boundary values = 0
eq = WavePDE(speed=0.5, bc=[[{'value':0}, {'value':0}], [{'value':0}, {'value':0}]])  

storage = pde.MemoryStorage()

trackers = [
    "progress",  # show progress bar during simulation
    "steady_state",  # abort when steady state is reached
    storage.tracker(interval=1),  # 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.RealtimeIntervals(duration=2)),
]
result = eq.solve(state, 20, dt=0.01, tracker=trackers)

I am using Jupyter Notebook on Windows 10 and when I try to run this code

import pde

grid = pde.UnitGrid([64, 64])                 # generate grid
state = pde.ScalarField.random_uniform(grid)  # generate initial condition

eq = pde.DiffusionPDE(diffusivity=0.1)        # define the pde
result = eq.solve(state, t_range=10)          # solve the pde
result.plot()                                 # plot the resulting field

it gives me the following error message. If you can help me I would really appreciate it.

TypingError Traceback (most recent call last)
in
5
6 eq = pde.DiffusionPDE(diffusivity=0.1) # define the pde
----> 7 result = eq.solve(state, t_range=10) # solve the pde
8 result.plot() # plot the resulting field

~\Anaconda3\lib\site-packages\pde\pdes\base.py in solve(self, state, t_range, dt, tracker, method, ret_info, **kwargs)
312
313 # run the simulation
--> 314 final_state = controller.run(state, dt)
315
316 if ret_info:

~\Anaconda3\lib\site-packages\pde\solvers\controller.py in run(self, state, dt)
142
143 # initialize the stepper
--> 144 stepper = self.solver.make_stepper(state=state, dt=dt)
145
146 # initialize profiling information

~\Anaconda3\lib\site-packages\pde\solvers\scipy.py in make_stepper(self, state, dt)
67 # obtain function for evaluating the right hand side
68 rhs = self._make_pde_rhs(state, backend=self.backend,
---> 69 allow_stochastic=False)
70
71 def rhs_helper(t, state_flat):

~\Anaconda3\lib\site-packages\pde\solvers\base.py in _make_pde_rhs(self, state, backend, allow_stochastic)
123 rhs = self.pde.make_sde_rhs(state, backend=backend)
124 else:
--> 125 rhs = self.pde.make_pde_rhs(state, backend=backend)
126
127 if hasattr(rhs, '_backend'):

~\Anaconda3\lib\site-packages\pde\pdes\base.py in make_pde_rhs(self, state, backend)
86 if backend == 'auto':
87 try:
---> 88 rhs = self._make_pde_rhs_numba(state)
89 except NotImplementedError:
90 backend = 'numpy'

~\Anaconda3\lib\site-packages\pde\pdes\diffusion.py in _make_pde_rhs_numba(self, state)
93
94 @jit(signature)
---> 95 def pde_rhs(state_data: np.ndarray, t: float):
96 """ compiled helper function evaluating right hand side """
97 return diffusivity_value * laplace(state_data)

~\Anaconda3\lib\site-packages\pde\tools\misc.py in (realf)
162 else:
163 # decorator arguments
--> 164 return lambda realf: decorator(realf, *args, **kwargs)
165
166 return new_decorator

~\Anaconda3\lib\site-packages\pde\tools\numba.py in jit(function, signature, parallel, **kwargs)
156 logging.getLogger(name).info('Compile %s with parallel=%s',
157 name, jit_kwargs['parallel'])
--> 158 return nb.jit(signature, **jit_kwargs)(function) # type: ignore
159
160

~\Anaconda3\lib\site-packages\numba\decorators.py in wrapper(func)
184 with typeinfer.register_dispatcher(disp):
185 for sig in sigs:
--> 186 disp.compile(sig)
187 disp.disable_compile()
188 return disp

~\Anaconda3\lib\site-packages\numba\compiler_lock.py in _acquire_compile_lock(*args, **kwargs)
30 def _acquire_compile_lock(*args, **kwargs):
31 with self:
---> 32 return func(*args, **kwargs)
33 return _acquire_compile_lock
34

~\Anaconda3\lib\site-packages\numba\dispatcher.py in compile(self, sig)
691
692 self._cache_misses[sig] += 1
--> 693 cres = self._compiler.compile(args, return_type)
694 self.add_overload(cres)
695 self._cache.save_overload(sig, cres)

~\Anaconda3\lib\site-packages\numba\dispatcher.py in compile(self, args, return_type)
78 return retval
79 else:
---> 80 raise retval
81
82 def _compile_cached(self, args, return_type):

~\Anaconda3\lib\site-packages\numba\dispatcher.py in _compile_cached(self, args, return_type)
88
89 try:
---> 90 retval = self._compile_core(args, return_type)
91 except errors.TypingError as e:
92 self._failed_cache[key] = e

~\Anaconda3\lib\site-packages\numba\dispatcher.py in _compile_core(self, args, return_type)
106 args=args, return_type=return_type,
107 flags=flags, locals=self.locals,
--> 108 pipeline_class=self.pipeline_class)
109 # Check typing error if object mode is used
110 if cres.typing_error is not None and not flags.enable_pyobject:

~\Anaconda3\lib\site-packages\numba\compiler.py in compile_extra(typingctx, targetctx, func, args, return_type, flags, locals, library, pipeline_class)
970 pipeline = pipeline_class(typingctx, targetctx, library,
971 args, return_type, flags, locals)
--> 972 return pipeline.compile_extra(func)
973
974

~\Anaconda3\lib\site-packages\numba\compiler.py in compile_extra(self, func)
388 self.lifted = ()
389 self.lifted_from = None
--> 390 return self._compile_bytecode()
391
392 def compile_ir(self, func_ir, lifted=(), lifted_from=None):

~\Anaconda3\lib\site-packages\numba\compiler.py in _compile_bytecode(self)
901 """
902 assert self.func_ir is None
--> 903 return self._compile_core()
904
905 def _compile_ir(self):

~\Anaconda3\lib\site-packages\numba\compiler.py in _compile_core(self)
888 self.define_pipelines(pm)
889 pm.finalize()
--> 890 res = pm.run(self.status)
891 if res is not None:
892 # Early pipeline completion

~\Anaconda3\lib\site-packages\numba\compiler_lock.py in _acquire_compile_lock(*args, **kwargs)
30 def _acquire_compile_lock(*args, **kwargs):
31 with self:
---> 32 return func(*args, **kwargs)
33 return _acquire_compile_lock
34

~\Anaconda3\lib\site-packages\numba\compiler.py in run(self, status)
264 # No more fallback pipelines?
265 if is_final_pipeline:
--> 266 raise patched_exception
267 # Go to next fallback pipeline
268 else:

~\Anaconda3\lib\site-packages\numba\compiler.py in run(self, status)
255 try:
256 event("-- %s" % stage_name)
--> 257 stage()
258 except _EarlyPipelineCompletion as e:
259 return e.result

~\Anaconda3\lib\site-packages\numba\compiler.py in stage_nopython_frontend(self)
513 self.args,
514 self.return_type,
--> 515 self.locals)
516 self.typemap = typemap
517 self.return_type = return_type

~\Anaconda3\lib\site-packages\numba\compiler.py in type_inference_stage(typingctx, interp, args, return_type, locals)
1122
1123 infer.build_constraint()
-> 1124 infer.propagate()
1125 typemap, restype, calltypes = infer.unify()
1126

~\Anaconda3\lib\site-packages\numba\typeinfer.py in propagate(self, raise_errors)
925 if errors:
926 if raise_errors:
--> 927 raise errors[0]
928 else:
929 return errors

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<function _make_laplace_numba_2d..laplace at 0x00000208DBD26708>) with argument(s) of type(s): (array(float64, 2d, C), out=array(float64, 2d, C))

• parameterized
  In definition 0:
  TypingError: Failed in nopython mode pipeline (step: nopython frontend)
  Invalid use of Function(<function BoundaryPair.make_region_evaluator..region_evaluator at 0x00000208DBD0E438>) with argument(s) of type(s): (array(float64, 2d, C), tuple(int64 x 2))
• parameterized
  In definition 0:
  TypingError: Failed in nopython mode pipeline (step: nopython frontend)
  Invalid use of Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>) with argument(s) of type(s): (array(float64, 1d, A), int64, (..., int64))
• parameterized
  In definition 0:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  In definition 1:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  This error is usually caused by passing an argument of a type that is unsupported by the named function.
  [1] During: resolving callee type: Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>)
  [2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py (379)

File "..\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py", line 379:
def region_evaluator(arr, idx: Tuple[int, ...])

arr_1d, i_point, bc_idx = get_arr_1d(arr, idx)
return (ap_low(arr_1d, i_point, bc_idx),
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

In definition 1:
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>) with argument(s) of type(s): (array(float64, 1d, A), int64, (..., int64))

• parameterized
  In definition 0:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  In definition 1:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  This error is usually caused by passing an argument of a type that is unsupported by the named function.
  [1] During: resolving callee type: Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>)
  [2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py (379)

File "..\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py", line 379:
def region_evaluator(arr, idx: Tuple[int, ...])

arr_1d, i_point, bc_idx = get_arr_1d(arr, idx)
return (ap_low(arr_1d, i_point, bc_idx),
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: resolving callee type: Function(<function BoundaryPair.make_region_evaluator..region_evaluator at 0x00000208DBD0E438>)
[2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\operators\cartesian.py (351)

File "..\Anaconda3\lib\site-packages\pde\grids\operators\cartesian.py", line 351:
def laplace(arr, out=None):

j = 0
arr_x_l, arr_c, arr_x_h = region_x(arr, (i, j))
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

In definition 1:
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<function BoundaryPair.make_region_evaluator..region_evaluator at 0x00000208DBD0E438>) with argument(s) of type(s): (array(float64, 2d, C), tuple(int64 x 2))

• parameterized
  In definition 0:
  TypingError: Failed in nopython mode pipeline (step: nopython frontend)
  Invalid use of Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>) with argument(s) of type(s): (array(float64, 1d, A), int64, (..., int64))
• parameterized
  In definition 0:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  In definition 1:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  This error is usually caused by passing an argument of a type that is unsupported by the named function.
  [1] During: resolving callee type: Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>)
  [2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py (379)

File "..\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py", line 379:
def region_evaluator(arr, idx: Tuple[int, ...])

arr_1d, i_point, bc_idx = get_arr_1d(arr, idx)
return (ap_low(arr_1d, i_point, bc_idx),
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

In definition 1:
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>) with argument(s) of type(s): (array(float64, 1d, A), int64, (..., int64))

• parameterized
  In definition 0:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  In definition 1:
  KeyError: "<class 'numba.targets.cpu.CPUTargetOptions'> does not support option: 'inline'"
  raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\targets\options.py:20
  This error is usually caused by passing an argument of a type that is unsupported by the named function.
  [1] During: resolving callee type: Function(<function BCBase1stOrder.make_adjacent_evaluator..adjacent_point at 0x00000208DBD0E828>)
  [2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py (379)

File "..\Anaconda3\lib\site-packages\pde\grids\boundaries\axis.py", line 379:
def region_evaluator(arr, idx: Tuple[int, ...])

arr_1d, i_point, bc_idx = get_arr_1d(arr, idx)
return (ap_low(arr_1d, i_point, bc_idx),
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: resolving callee type: Function(<function BoundaryPair.make_region_evaluator..region_evaluator at 0x00000208DBD0E438>)
[2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\grids\operators\cartesian.py (351)

File "..\Anaconda3\lib\site-packages\pde\grids\operators\cartesian.py", line 351:
def laplace(arr, out=None):

j = 0
arr_x_l, arr_c, arr_x_h = region_x(arr, (i, j))
^

raised from C:\Users\matth\Anaconda3\lib\site-packages\numba\typeinfer.py:927

This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: resolving callee type: Function(<function _make_laplace_numba_2d..laplace at 0x00000208DBD26708>)
[2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\tools\numba.py (259)

File "..\Anaconda3\lib\site-packages\pde\tools\numba.py", line 259:
def f_with_allocated_out(arr, out):

""" helper function allocating output array """
return f_jit(arr, out=np.empty_like(arr))
^

[1] During: resolving callee type: type(CPUDispatcher(<function _make_laplace_numba_2d..laplace at 0x00000208DBD268B8>))
[2] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\pdes\diffusion.py (97)

[3] During: resolving callee type: type(CPUDispatcher(<function _make_laplace_numba_2d..laplace at 0x00000208DBD268B8>))
[4] During: typing of call at C:\Users\matth\Anaconda3\lib\site-packages\pde\pdes\diffusion.py (97)

File "..\Anaconda3\lib\site-packages\pde\pdes\diffusion.py", line 97:
def pde_rhs(state_data: np.ndarray, t: float):

""" compiled helper function evaluating right hand side """
return diffusivity_value * laplace(state_data)
^

Hello,
I am just getting started with py-pde, I am attempting to go through many of the examples. I did a clean install on a fresh environment in conda, but I am getting a similar error on some, but not all the examples. Here is the error, for example on the "Simple Diffusion Equation" example:

`~/anaconda3/envs/pdeenv/lib/python3.7/site-packages/pde/fields/base.py in _apply_operator(self, operator, bc, out, label, **kwargs)
1620 # apply the operator with boundary conditions
1621 op_with_bcs = self.grid.make_operator(operator_info, bc=bc, **kwargs)
-> 1622 out.data[:] = op_with_bcs(self.data)
1623
1624 return out

TypeError: expected dtype object, got 'numpy.dtype[float64]'`

Dear Dr. Zwicker,

I tried some examples by py-pde today and found something different about plotting. For the previous version, for any given t_range and interval, say t_range=20 and interval=1, the program will draw the pictures of solutions at every interge time from 0 to 20. In the current version, if I still use the same setting, the program will only draw the pictures of solutions at some certain integer times but not every step.
How could I do for recovering the previous way, like plotting for every time step?

Thank you so much!

Looking at the SIR model example what would be the easiest way to limit the values of one of the fields?

By means of evolution_rate the three fields s, i and r evolve.
What if one of those three fields cannot exceed some value because of some constraint?

Say that i cannot exceed some value alpha. So if i became larger than alpha it would have to be set equal to alpha.

How would you implement that?

Hi,

it seems that the FileStorage class can only handle real-valued fields, a simple modification of the Schrödinger equation example to use FileStorage instead of MemoryStorage should reproduce the problem. I think it stems from the fact that in _create_hdf_dataset() of FileStorage in pde/storage/file.py the dtype is always the default np.double.
If anybody else has the same problem, a workaround might be to subclass StorageTracker and reimplement the handle() function to store only the absolute value of the field instead of complex.
Would be cool if this can be fixed, the package looks very nice with a very good high-level API, I like it a lot!

Greetings from Austria,
Gregor

Dear Dr Zwicker

Thanks for creating this public pde solver. I am new to using it but it looks great!

I have an issue trying to simulate an exponential decay in 2 dimensions. The code is
"
import numpy as np
from pde import PDE, FieldCollection, PlotTracker, ScalarField, UnitGrid, CartesianGrid
import functools

tauE=5
eq = PDE(
{
"E": f"-E/{tauE}",
},
bc="periodic"
)

#Define Cartesian Grid
grid = CartesianGrid([[0, 1], [0, 1]], [32, 32], periodic=[True,True]) # generate grid
E = ScalarField(grid,1, label="Field $E$")
state = FieldCollection([E])

simulate the pde

tracker = PlotTracker(interval=0.5, plot_args={"vmin": -3, "vmax": 3})
sol = eq.solve(state, t_range=10, dt=1e-5, tracker=tracker)
"
I expect the value at every point in the grid to decay to 0. However instead the value decays to 0 and then blows up to large negative values. Can you see where the error is in my code?

after i install the package, i try to run the example code 2.1. Plotting a vector field. i copied the code and put it into my jupyter notebook, i got this message
module 'sympy.core.core' has no attribute 'basic'

many other examples can not run, can anyone help me figure out why and solve it?

thank you very much.

[openjournals/joss-reviews#2158]

For consistency, I'd suggest mentioning the dependencies for running tests already on the Getting Started page of the documentation (in the same way as for the documentation dependencies).

Given that "nobody reads the documentation", you might also want to consider a check for the presence of pytest etc. directly in the shell scripts with a helpful error message for the user. For example, if the user does not have pytest-xdist installed, running tests_parallel.sh will currently fail with "unknown option -n".

Thank you for publishing this package. However, after reading through the documentation I don't think I can use it to model my system. Hopefully you can tell me I'm mistaken!

I need to model a system where a sequential reaction (e.g. compound 1 ➞ compound 2 ➞ compound 3) is occurring in the aqueous phase while the chemical compounds involved in that sequential reaction are also diffusing into a solid phase. The reaction only occurs in the aqueous phase and the solid phase is only present on one wall of the reaction vessel, so the problem involves the 1D diffusion equation.

When there is no solid phase present I can model the system using scipy.integrate.solve_ivp to solve a set of ODEs describing the reaction kinetics:

where c is a vector of aqueous chemical concentrations.

Can this package:

1. Perform multiple solutions of the 1D diffusion equation (one for each chemical compound) in parallel?
2. Incorporate Neumann boundary conditions that are a determined by a system of ODEs?

Hi,
I am trying to solve a simple 2D advection-diffusion equation using py-pde

The equation is:

I understand how to implement first part of the equation with:

eq = pde.PDE({"C": f"{D} * laplace(C)"})

But cannot find any description or examples that deal with vector fields to implement the divergence part of the equation.

I am asking for pointers on how to implement the advection part in py-pde, if possible. I assume I need to write my own PDEBase class but I can't find a comparable example in docs.

Best regards.

When I try to use the wave equation the amplitude of the wave is amplified by a lot over the duration of the simulation. Here is my code:

import numpy as np
from pde import FieldCollection
from pde import (DiffusionPDE, UnitGrid, ScalarField, MemoryStorage,
                 PlotTracker, PrintTracker, RealtimeIntervals, DataTracker)
from pde.visualization import movies
import math
import random
import pde
import pandas
import ffmpeg
import matplotlib.pyplot as plt; plt.rcdefaults()
import numpy as np
import matplotlib.pyplot as plt


storage = MemoryStorage()

trackers = [
    'progress',                    # show progress bar during simulation
    'steady_state',                # abort when steady state is reached
    storage.tracker(interval=1),   # store data every simulation time unit
    PlotTracker(show=True),        # show images during simulation
]

grid = pde.UnitGrid([64, 64])
field1 = ScalarField.random_uniform(grid, 0, 0)  # generate initial condition
field2 = ScalarField(grid);
field1.add_interpolated([48, 48], 1)
field1.add_interpolated([48,16], 1)
state = FieldCollection([field1,field2]); #first slot is the density field, 

eq = pde.WavePDE(speed=1);
result = eq.solve(state, t_range=200, dt=1, tracker = trackers)
result.plot(colorbar = True, show = True)       # plot the resulting field
movies.movie_multiple(storage, filename = r'C:\Users\matth\URECA Simulation\testfile.mp4')

and here is the file of the movie. In the movie the colorbar is at 10^90 and the wave can only be seen at the end of the movie, which means that the amplitude is increasing throughout the simulation. I looked at the documentation and I am not sure why this is happening.

from math import sqrt,pi
import pde
from pde import WavePDE, CartesianGrid, ScalarField, MemoryStorage, FieldCollection, PDE, plot_kymograph, PlotTracker, plot_magnitudes, FileStorage
import matplotlib.pyplot as plt


#Set the Geometry
#Generate grid
grid = CartesianGrid( [[-1, 1],[-1,1]] , 256 )


#Initial condition
p=2
q=1
state = ScalarField.from_expression(grid, f"sin({p}*pi*x)*sin({q}*pi*y)" )

#Boundary condition
#specified within eq call
nu =0.01/pi
eq = PDE(
{
        "u": f"{nu} * (getx(laplace(u))+gety(laplace(u))) - u * (getx(gradient(u))+gety(gradient(u)))"
}, user_funcs={"get_x": lambda arr: arr[0]},{"get_y": lambda arr: arr[1]}],
    
    bc={"value": 0}
 )

# store intermediate information of the simulation
storage = FileStorage("2D_burger_data.npz")

# solve the pde
result = eq.solve(state, t_range = 0.99, tracker=storage.tracker(0.01))

pde.movie(storage, filename="2D_burger_movie.mov")
# plot the resulting field
result.plot()

# visualize the result in a space-time plot
plot_kymograph(storage)

# slice data
data=state
smoothed = data.smooth()  # Gaussian smoothing to get rid of the noise
sliced = smoothed.slice({"y": 0})

# create two plots of the field and the modifications
fig, axes = plt.subplots(nrows=1, ncols=2)
data.plot(ax=axes[0], title="Original field")
sliced.plot(ax=axes[1], title="Slice of smoothed field at $y=1$")
plt.subplots_adjust(hspace=1.5)

Thank you for all your efforts! I'm trying to solve the 2D burger's equation. I ran into this error:

 user_functions.update(user_funcs)
AttributeError: 'list' object has no attribute 'update'

What is the right way to define and solve the equation?

Here's my reference:

In my case v is just a constant :0.01/pi.

Thanks in advance!

Matplotlib reports a broken pipe when the plot tracker should write a movie to a folder that does not exist. We should catch this and emit a more useful error.

Hi,

I'd like to try py-pde for some school projects, but I'm not able to make it plotting any results. This is the only thing I get, no error message, nothing. I pip installed py-pde, version is 0.12.1, all other packages are the latest ones.

One of examples in Jupyter Notebook:

In Jupyterlab I get extra:
<pde.tools.plotting.PlotReference at 0x7feef0b071c0>

Seems like something is wrongly configured at my side, but I can't find what it is :-( Could you please help me with this?

Edit: It works fine in ipython.

Hello, I have a simulation I run with mpiexec -n 16 python py_pde_test.py

every other time I run it I get the error File data\storage-1108_191954\ could not be opened for reading from some of the mpi threads.

The code:

import numpy as np
import numba as nb
import pde
from pde.tools import mpi
from pde.solvers import Controller, ExplicitMPISolver
import time


class ModelPDE(pde.PDEBase):
    def __init__(self, bc="auto_periodic_neumann", nu=10 / 3, eta=3.5, rho=0.95, gamma=50 / 3, delta_b=1 / 30,
                 delta_w=10 / 3, delta_h=1000 / 3, a=33.33, q=0.05, f=0.1, p=0.5):
        self._nu = nu
        self._eta = eta
        self._rho = rho
        self._gamma = gamma
        self._delta_b = delta_b
        self._delta_w = delta_w
        self._delta_h = delta_h
        self._a = a
        self._q = q
        self._f = f
        self._p = p

        self._bc = bc

    def _make_pde_rhs_numba(self, state):
        a = self._a
        q = self._q
        f = self._f
        eta = self._eta
        nu = self._nu
        rho = self._rho
        gamma = self._gamma
        delta_b = self._delta_b
        delta_w = self._delta_w
        delta_h = self._delta_h
        p = self._p

        laplace = state.grid.make_operator("laplace", bc=self._bc)
        div = state.grid.make_operator("divergence", bc=self._bc)
        grad = state.grid.make_operator("gradient", bc=self._bc)

        zeta = np.fromfunction(lambda x, _: x/32, state.grid.shape) # TODO: this is bad

        @nb.jit
        def pde_rhs(state_data, t):
            # state
            b = state_data[0]
            w = state_data[1]
            h = state_data[2]

            rate = np.empty_like(state_data)
            # Calculate constants
            I = a * (b + q * f) / (b + q)  # Infiltration Rate
            L2 = np.float_power(1 + eta * b, 2)
            Gb = nu * w * L2
            Gw = gamma * b * L2

            # Calculate time derivatives
            rate[0] = Gb * b * (1 - b) - b + delta_b * laplace(b)
            rate[1] = I * h - nu * (1 - rho * b) * w - Gw * w + delta_w * laplace(w)

            J = -2*delta_h*h*grad(h+zeta)
            rate[2] = p - I * h - div(J)

            return rate

        return pde_rhs

    def evolution_rate(self, state: pde.FieldBase, t: float = 0) -> pde.FieldBase:
        b = state[0]
        w = state[1]
        h = state[2]

        zeta = np.fromfunction(lambda x, _: x/32, state.grid.shape) # TODO: this is bad

        # Calculate constants
        I = self._a * (b + self._q * self._f) / (b + self._q)  # Infiltration Rate
        L2 = np.float_power(1 + self._eta * b, 2)
        Gb = self._nu * w * L2
        Gw = self._gamma * b * L2

        # Calculate time derivatives
        b_t = Gb * b * (1 - b) - b + self._delta_b * b.laplace(bc=self._bc)
        w_t = I * h - self._nu * (1 - self._rho * b) * w - Gw * w + self._delta_w * w.laplace(bc=self._bc)

        J = -2*self._delta_h*h *(h + zeta).gradient(bc=self._bc)
        h_t = self._p - I * h - J.divergence(bc=self._bc)

        return pde.FieldCollection([b_t, w_t, h_t])


def terrain(coords):
    return coords[:, :, 0] / 32

def main():
    shape = (128, 128)
    # grid = pde.UnitGrid([128, 128], periodic=[True, True])
    grid = pde.CartesianGrid([(0, 32), (0, 32)], shape, periodic=[True, False])
    b = pde.ScalarField(grid, 1)
    w = pde.ScalarField(grid, 0)
    h = pde.ScalarField(grid, 0)
    state = pde.FieldCollection([b, w, h])

    bc_x = "periodic"
    bc_y = [{"value": 0} ,{"curvature": 0}]
    bc = [bc_x, bc_y]

    years = 2

    t = time.localtime()
    timestamp = time.strftime("%m%d_%H%M%S", t)
    BACKUP_NAME = "data\storage-" + timestamp
    storage = pde.FileStorage(BACKUP_NAME, write_mode="append")

    eq = ModelPDE(p=0.4, bc=bc)

    solver = ExplicitMPISolver(eq, backend="numba")
    controller = Controller(solver, t_range=years, tracker=["progress", storage.tracker(1e-1)])
    sol = controller.run(state, dt=1e-4)
    if mpi.is_main:
        VIDEO_NAME = f"results\movie-{timestamp}.mp4"
        pde.movie(storage, filename=VIDEO_NAME, plot_args={"vmin": 0, "vmax": 1})


if __name__ == '__main__':
    main()

u_t = a^2 (u_xx + u_yy + u_zz)

https://en.wikipedia.org/wiki/Heat_equation

Hello!

I'm a beginner using py-pde. I'm wondering how I can define custom Dirichlet conditions for a grid. Is it possible to pass a numpy array with my config? See the image below for clarity.

Thanks in advance!

Hi, @david-zwicker and collaborators!

Firstly, I would like to say that py-pde is a really nice and well written Python package. The API is a charm. Congratulations for this work.

Well, getting back to business. Would you be interested in extending installation options with conda? I put py-pde there, using the conda-forge channel. Instructions to install py-pde using conda is provided in this feedstock.

If you have interest, I would happily open a PR updating installation instructions.

Cheers!
Diego

Hi!

I am trying to solve the heat equation on a grid with polar symmetry subject to the boundary condition ∂u/∂n = -0.5 · (u − 1). (This boundary condition is equivalent to Newton’s law of cooling with a heat-transfer coefficient of 0.5 and an ambient temperature of 1.)

The following four examples illustrate my problem:

First I try to do this on a one-dimensional Cartesian grid:

grid = CartesianGrid([[0, 1]], [32], periodic=False)
state = ScalarField.from_expression(grid, '0')
eq = PDE({'T': 'laplace(T)'}, bc=[{'derivative': 0}, {'derivative_expression': '-.5*(T-1)'}])
storage = MemoryStorage()
result = eq.solve(state, t_range=3, tracker=['progress', storage.tracker(.1)])

I get an error message saying that T is not in the expression signature ['value', 'dx', 'x', 't']. I assume value means the value of the solution at the boundary. Substituting value for T, I get a reasonably-looking solution:

grid = CartesianGrid([[0, 1]], [32], periodic=False)
state = ScalarField.from_expression(grid, '0')
eq = PDE({'T': 'laplace(T)'}, bc=[{'derivative': 0}, {'derivative_expression': '-.5*(value-1)'}])
storage = MemoryStorage()
result = eq.solve(state, t_range=3, tracker=['progress', storage.tracker(.1)])
for data in storage.data:
    plt.plot(data)

Now I repeat for a polar grid:

grid = PolarSymGrid(radius=1, shape=32)
state = ScalarField.from_expression(grid, '0')
eq = PDE({'T': 'laplace(T)'}, bc=[{'derivative': 0}, {'derivative_expression': '-.5*(T-1)'}])
storage = MemoryStorage()
result = eq.solve(state, t_range=3, tracker=['progress', storage.tracker(.1)])

I get the same error message: T is not in the expression signature ['value', 'dx', 'r', 't']. However, when I substitute value for T:

grid = PolarSymGrid(radius=1, shape=32)
state = ScalarField.from_expression(grid, '0')
eq = PDE({'T': 'laplace(T)'}, bc=[{'derivative': 0}, {'derivative_expression': '-.5*(value-1)'}])
storage = MemoryStorage()
result = eq.solve(state, t_range=3, tracker=['progress', storage.tracker(.1)])

I get an error message saying that the interpreter failed in nopython mode pipeline. Although this issue is similar to issue 12, my numba version is up to date at version 0.56.2 as seen by calling pde.environment():

{'package version': '0.21.0',
 'python version': '3.8.10 (default, Jun 22 2022, 20:18:18) \n[GCC 9.4.0]',
 'platform': 'linux',
 'config': {'numba.debug': False,
  'numba.fastmath': True,
  'numba.parallel': True,
  'numba.parallel_threshold': 65536},
 'mandatory packages': {'matplotlib': '3.3.3',
  'numba': '0.56.2',
  'numpy': '1.19.4',
  'scipy': '1.6.3',
  'sympy': '1.7.1'},
 'matplotlib environment': {'backend': 'module://ipykernel.pylab.backend_inline',
  'plotting context': 'JupyterPlottingContext'},
 'optional packages': {'h5py': '3.7.0',
  'napari': 'not available',
  'pandas': '1.2.1',
  'pyfftw': 'not available',
  'tqdm': '4.64.1'},
 'numba environment': {'version': '0.56.2',
  'parallel': True,
  'fastmath': True,
  'debug': False,
  'using_svml': False,
  'threading_layer': 'omp',
  'omp_num_threads': None,
  'mkl_num_threads': None,
  'num_threads': 4,
  'num_threads_default': 4,
  'cuda_available': False,
  'roc_available': False}}

The full error message is attached as a file.

Is this an error in py-pde? Is there a workaround for polar grids?

(Note that I have used a user-defined PDE in these examples. However, using the predefined DiffusionPDE gives the same results for the second example both in the Cartesian and polar case.)

Thanks for any help to resolve this!

Hi,
I didn't see it mentioned anywhere and it doesn't appear in any of the examples.
Can I define a PDE of the form:

where D(x) could be a polynomial of x for example

Thanks
Itay

There are currently two places where expressions in forms of strings can be used:

1. We can treat expressions using pde.tools.expressions.ScalarExpression, which represent functions of multiple variables, but have no notion of space (they do not depend on any grid).
2. We can express the right hand side of a PDE as an expression using the PDE class, which also allows using differential operators.

It would be great to have some way to evaluate expressions involving fields. The syntax could be something like

pde.evaluate("expression_string", {"field": ScalarField})

where the first argument would be the expression to evaluate (which might contain differential operators) and the second argument would be a dictionary of fields that are involved in the expression. We might not need to compile the expression, assuming that it is not called often. One use case would be to test expressions before actually using them in PDE.

Discussed in #191

I am trying to use $\delta(x)$ as initial condition but i am getting the following error:

ValueError Traceback (most recent call last)
ValueError: Error from parse_expr with transformed code: "Symbol ('sympy' ).DiracDelta (Symbol ('x' ))"

The above exception was the direct cause of the following exception: AttributeError: 'Symbol' object has no attribute 'DiracDelta'

Below is my code. What am I doing wrong? Even if i use 'DiracDelta(x)' instead of 'sympy.DiracDelta' since expressions are passed using the sympy package, I still have the same error. I need help please.

grid = pde.CartesianGrid([[-1, 1]], 64) # generate grid
field = pde.ScalarField.from_expression(grid,'sympy.DiracDelta(x)') # generate initial condition

eq = pde.PDE({"p": "d_dx(p * (a+(b*x))) + (0.05 * laplace(p))"},
consts={'a':0.02,
'b':0.13},
bc={"value": 0})
res = eq.solve(field, 15, dt=0.01,)

I am trying to create a movie in Jupyter Notebook on Windows 10 using the following code

import numpy as np
from pde import FieldCollection
from pde import (DiffusionPDE, UnitGrid, ScalarField, MemoryStorage,
                 PlotTracker, PrintTracker, RealtimeIntervals)
from pde.visualization import movies
import pde

storage = MemoryStorage()

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

grid = pde.UnitGrid([64, 64])
field2 = ScalarField.random_uniform(grid, 1, 1)  # generate initial condition
field1 = ScalarField.random_harmonic(grid)  # generate initial condition
state = FieldCollection.scalar_random_uniform(2, grid)
state = FieldCollection([field1,field2]);


eq = pde.WavePDE(speed=0.1);
result = eq.solve(state, t_range=3, dt=1, tracker = trackers)
result.plot(show=True)                                 # plot the resulting field
movies.movie_multiple(storage, filename = r'C:\Users\matth\URECA Simulation\testfile.mp4')

it gives me this error:
FileNotFoundError: [WinError 2] The system cannot find the file specified.
I have tried saving it to files without a path, files with a path, and files in the jupyter notebook environment but I have been unable to get it to work.

[openjournals/joss-reviews#2158]

The documentation states "The code is tested under Linux, Windows, and macOS", but the repository currently only runs the test suite on Linux via Travis. Travis also supports testing under OS X, and there are other (free for FOSS) services like https://www.appveyor.com/ or https://azure.microsoft.com/en-us/services/devops/pipelines/ that test across all three OS.

Can't resolve a system of partial differential equations with numba acelleration.

Going to solve PDE system
dCa/dt = - Fv * dCa/dV - k * Ca
dCb/dt = - Fv * dCb/dV + k * Ca

where the term k * Ca comes from a pre-defined user function.
The PDE system is implemented using a custom PDE class.
The resolution of the problem using the method evolution_rate(self, state, t=0) works (Code below).
Turn reactor.resolve(Numba=False), at the end, to not use numba.

Thank you very much! :)

from pde import (PDEBase, FieldCollection, CartesianGrid, ScalarField, 
                 ImplicitSolver, Controller
                )
from pde.tools.numba import jit
from numba import njit
import matplotlib.pyplot as plt
from pde import config

config['numba.debug'] = True
# =============================================================================
#  Clases
# -----------------------------------------------------------------------------

class ReactiveSystem:
    def __init__(self, ks, inflows, reactions):
        self.ks = ks
        self.inflows = inflows
        self.reactions = reactions
        self.net_flow = sum(inflows)

    @jit
    def react(self, state):
        return [
            r(self.ks, state) for r in self.reactions
        ]
    
    def __len__(self):        
        return len(self.reactions)


class ReactorSolver(PDEBase):
    def __init__(
        self, reactive_system, V, bcs,
        n_points=10, nt=10, t_end=10, t=0, max_it=1000, tol=1e-5
    ):
        # Number of reactives (hence, balances)
        self.n_reactives = len(reactive_system)

        # System dimension
        grid = CartesianGrid([[0, V]], n_points)
        states = [
            ScalarField(grid, data=0)
            for i in range(self.n_reactives)
        ]
        self.state = FieldCollection(states)

        # Border conditions
        self.bcs = bcs

        # Reactive system
        self.reactive_system = reactive_system

        # Time ranges
        self.nt = nt
        self.t_end = t_end
        self.dt = t_end/nt

        # Solver tolerance
        self.max_iter = max_it
        self.tol = tol       
    
    def evolution_rate(self, state, t=0):
        # Agregué el índice cero en el gradiente, hay que indicar que es
        #  la derivada en X más allá de que sea la única variable

        # Calculate reaction rates
        n_reactions = len(self.reactive_system)
        ks = self.reactive_system.ks
        reactions = [
            react(ks, state) for react in self.reactive_system.reactions
            ]

        # Net flow (of course it shouldn't be like this)
        net_flow = self.reactive_system.net_flow

        # Whole balances
        system = [
            - net_flow * state[i].gradient(self.bcs[i])   # Flow balance
            + reactions[i]                                # Reactions
            for i in range(self.n_reactives)
        ]
        return FieldCollection(system)

    def _make_pde_rhs_numba(self, state):
        #attributes locally available

        #Net flow (this is bad implemented - supposed to be vol flow)
        net_flow = self.reactive_system.net_flow

        #Number of reactives
        n_reactives = len(self.reactive_system)

        #kinetic constans
        ks = self.reactive_system.ks

        #Reaction function of each compound

        reactions = [
            react(ks, state).data for react in self.reactive_system.reactions
            ]
        
        #Operands definition
        gradient_c = [
            state.grid.make_operator("gradient", bc=self.bcs[i]) for i 
            in range(n_reactives)
            ]

        @jit
        def pde_rhs(state_data, t=0):
            """ compiled helper function evaluating right hand side """
            state_grad = []

            for i in range(n_reactives):
                gradient = gradient_c[i]
                state_grad.append(gradient(state_data[i]))

            #Reaction function of each compound
          
            """system = [
                - net_flow * state_grad[i]   # Flow balance
                + reactions[i]               # Reactions
                for i in range(n_reactives)
                ]"""

            system = []
            for i in range(n_reactives):
                system.append(- net_flow * state_grad[i])
            return system

        return FieldCollection(pde_rhs)

    def resolve(self, numba = False):
        if numba == False:
            solver = ImplicitSolver(self, self.max_iter, self.tol, 
                                    backend='numpy')
            controller = Controller(solver, t_range=self.t_end, 
                                    tracker=None)
            return controller.run(self.state)
        else:
            solver = ImplicitSolver(self, self.max_iter, self.tol, 
                                    backend='numba')
            controller = Controller(solver, t_range=self.t_end, 
                                    tracker=None)
            return controller.run(self.state)

# =============================================================================
# =============================================================================
#  Planteo de problema
# -----------------------------------------------------------------------------

# Variables generales
kr = 1
t0 = 0
tf = 10
V = 1
Fv = 1
Ca0 = 1

# Definicion de sistema reactivo
rs = ReactiveSystem(
    ks=[kr],
    inflows=[Fv, 0],
    reactions=[
        lambda ks, concentrations: -ks[0] * concentrations[0],
        lambda ks, concentrations:  ks[0] * concentrations[0],        
    ]
)


# Border conditions
##  - x0 = CaO
##  - xf = dCa/dx = 0
bc_a = [
    {"value": Ca0},
    {"derivative": 0}
    ] 
bc_b = [
    {"value": 0},
    {"derivative": 0}
]

# Reactor
reactor = ReactorSolver(
    reactive_system=rs,
    V=V,
    bcs=[bc_a, bc_b],
    n_points = 100,
    t_end=10,
    nt = 1000,
)
# ==============================================================================

results = reactor.resolve(numba=True)


plt.plot(results[0].data)
plt.plot(results[1].data)
plt.show()

The latest pytest version does have a bug where conftest.py is not discovered properly: pytest-dev/pytest#9767
This breaks our tests, so we temporarily fixed the pytest version to 7.0.1 in tests/requirements.txt and .github/workflows/tests_minversion.yml.
This change should be reverted once pytest is fixed.

Are there any parameters in method field.plot() to output clear pictures for visualizing physical fields, such as dpi=200 ?

Currently, we can only pass a scalar diffusivity. It'd be very useful if we could optionally pass a scalar field for diffusivity (each element corresponding to a single point/control-volume in the grid).

Discussed in #120

I am trying to solve 4th order ODEs containing a biharmonic operator. The documentation describes solving the Kuramoto-Sivashinsky equation and the Cahn-Hilliard equation.

However, there is no option to specify enough boundary conditions to solve a 4th order ODE. For a linear ODE we need 4 and this is generally true for nonlinear ODEs also. With Neumann BCs I expect to be able to specify values and derivatives at the same point, with e.g. something like for a 2d system:
bcs = [{'value': -1, 'derivative': 0}, {'value': 1, 'derivative': 0}]
However, this always throws an error saying this feature is not supported. So I am curious how 4th order ODEs are listed as examples. Ultimately I'm interested in something similar to the Cahn-Hilliard equation so look to that as an example.

Looking at the implementation of the Cahn-Hilliard equation, it looks like this would only work with periodic boundary conditions because this implicitly sets values and all derivatives. That being said, I am never able to get convergent results for e.g. the binodal so I suspect either I am using py-pde wrong or there is a problem with the underlying solver. Below is a minimal working example that illustrates this where I try to obtain the binodal. I start with two interfaces around x=+-0.5 to respect periodic boundaries at x=+-1. The solver quickly fails with a "Field was not finite" error:

from pde import CahnHilliardPDE, ScalarField, CartesianGrid
grid = CartesianGrid([[-1, 1], [-1, 1]], [1000,50])
state = ScalarField.from_expression(grid, 'tanh(cos(pi*x)/0.1)')
eq = CahnHilliardPDE()
result = eq.solve(state, 1, dt=1e-6, method='implicit')
result.plot()

My expected behaviour is that it should converge to tighten up the interfaces with the correct interface width (rather than the dummy value of 0.1 I gave as an IC).

Is this a bug, or am I using it wrong? If the latter could you please advise me on how I should use the Cahn Hilliard (and similar) solvers? Thank you.

Thanks for the very beneficial code. I'm actually using it as one of my base stones in my research!

I 'm trying to implement the wave equation with a time-dependent source placed at the boundary.

speed = 0.01
frequency = 10
source = f"sin(pi*x) * exp(-2*pi*{frequency}*t)"
eq = PDE(
    {
        "u": f"v",
        "v": f"{speed}**2 * laplace(u) " 
    }, bc=[{"value": f"sin(pi*x) * exp(-2*pi*{frequency}*t)"}, {"derivative":0}]
)

but I get this error:

 f"Arguments {args} were not defined in expression signature {signature}"
RuntimeError: Arguments {'t'} were not defined in expression signature ['x']

It will be a novel addition to your code if it becomes supported!

Thanks again.

Hi,
I was wondering if it's possible to define custom geometries in py-pde, for insance a 2D grid with a solid block take out of the domain (see below):

o o o o o o o o o o
o o o o o o o o o o
o o o x x x x o o o
o o o x x x x o o o
o o o x x x x o o o

Thanks!

Thank you for the magnificent work in py-pde!

I'm simulating the 1D wave equation (with a time function source) in py-pde, and I want to compare the results to another results I got from another python code output.
I want to know how can I use FileStorage to save the fields:
x, t , u

• How can I store the data in npz/ h5py file?
• How can I extract the data in my other code?

In my other code, I'm saving my data in an npz file:

data_file = 'mydata'
np.savez(data_file, t_input=t, x_input=x, y_input=y, u_output=exact )

and then loading them to check/plot the values:

data = np.load(data_file+'.npz')
t_i, x_i, y_i, exact_i = data["t_input"], data["x_input"], data["y_input"], data["u_output"]

How can I do something similar in py-pde?

Thanks in advance!

This my code:

from math import sqrt,pi
import pde
from pde import WavePDE, CartesianGrid, ScalarField, MemoryStorage, FieldCollection, PDE, plot_kymograph, PlotTracker

#Set the Geometry
#Generate grid
grid = CartesianGrid( [[-1, 1]] , 256 )


#Initial condition
u = ScalarField.from_expression(grid, "sin(pi*x)", label="Field $u$" )
v = ScalarField(grid, 0, label="Field $v$")
initial_state = FieldCollection([u, v])


#Boundary condition
#specified within eq definition
#eq = pde.WavePDE(speed = 0.01, bc=[{"value": 0}, {"derivative": 0}])
speed = 0.01
frequency = 10
source = f"sin(2*pi*{frequency}*t)"
eq = PDE(
    {
        "u": f"v",
        "v": f"{speed}**2 * laplace(u) + {source}"
    }, bc=[{"value": 0}, {"derivative":0}]
)

# store intermediate information of the simulation
storage = MemoryStorage()

# solve the pde
result = eq.solve(initial_state, t_range = 1, tracker=[storage.tracker(0.01)])
    
#create a simulation movie
pde.movie(storage, filename="1D_wave_movie.mov")

# plot the resulting field
result.plot()

# visualize the results as a space-time plot
plot_kymograph(storage, scalar="norm_squared")