convection equation

[vnikoofard]

Hi,

trying to solve the convection equation with the following code

import numpy as np
import pde
from pde import PDE, CartesianGrid, MemoryStorage, ScalarField

nx = 151  # number of spatial discrete points
L = 2.0  # length of the 1D domain
dx = L / (nx - 1)  # spatial grid size
nt = 10  # number of time steps
dt = 0.001  # time-step size
x = np.linspace(0.0, L, num=nx)
u0 = np.ones(nx)
mask = np.where(np.logical_and(x >= 0.5, x <= 1.0))
u0[mask] = 2.0

grid = CartesianGrid([[0., 2.]], nx)
field = ScalarField(grid, data=u0)

eq = PDE({"u": "-gradient(u)"})

# solve the equation and store the trajectory
storage = MemoryStorage()
res = eq.solve(field, t_range=5, tracker=storage.tracker(0.1), 
             dt=dt, 
             method='scipy',
             
             )

# plot the trajectory as a space-time plot
res.plot()

I receive the following error

AssertionError                            Traceback (most recent call last)
Cell In[14], line 5
      3 # solve the equation and store the trajectory
      4 storage = MemoryStorage()
----> 5 res = eq.solve(field, t_range=5, tracker=storage.tracker(0.1), 
      6              dt=dt, 
      7              method='scipy',
      8              
      9              )
     11 # plot the trajectory as a space-time plot
     12 res.plot()

File [~/miniconda3/lib/python3.10/site-packages/pde/pdes/base.py:597](https://file+.vscode-resource.vscode-cdn.net/home/vahid/GDrive/UERJ/My_Lectures/numerical_methods_using_python/~/miniconda3/lib/python3.10/site-packages/pde/pdes/base.py:597), in PDEBase.solve(self, state, t_range, dt, tracker, method, ret_info, **kwargs)
    594 controller = Controller(solver, t_range=t_range, tracker=tracker)
    596 # run the simulation
--> 597 final_state = controller.run(state, dt)
    599 # copy diagnostic information to the PDE instance
    600 if hasattr(self, "diagnostics"):

File [~/miniconda3/lib/python3.10/site-packages/pde/solvers/controller.py:322](https://file+.vscode-resource.vscode-cdn.net/home/vahid/GDrive/UERJ/My_Lectures/numerical_methods_using_python/~/miniconda3/lib/python3.10/site-packages/pde/solvers/controller.py:322), in Controller.run(self, initial_state, dt)
    319 # decide whether to call the main routine or whether this is an MPI client
    320 if mpi.is_main:
    321     # this node is the primary one
--> 322     self._run_single(state, dt)
...
         -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,
         -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,...
 y: array([ -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,
        -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,
        -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,  -0.  ,...

I tried different solvers but the erro persists.

[david-zwicker]

The issue is that you define the time-derivative of a scalar field using a vectorial quantity (the gradient of the same field). This is simply incompatible, but unfortunately the error message is difficult to improve because of the automatic parsing of the right hand side. Since you're considering a 1D problem, you should be able to fix the problem using eq = PDE({"u": "-d_dx(u)"}).

[vnikoofard]

The issue is that you define the time-derivative of a scalar field using a vectorial quantity (the gradient of the same field). This is simply incompatible, but unfortunately the error message is difficult to improve because of the automatic parsing of the right hand side. Since you're considering a 1D problem, you should be able to fix the problem using eq = PDE({"u": "-d_dx(u)"}).

Thanks for your help. You are right.

I don't know if I have to open another topic or here I can about the result of the code. Applying your suggestion the code works but the result is strange. Even with the implicit method there are instability in the final result.

[david-zwicker]

This is a common issue with advective problems. As far as I know there is no completely general solution, which ensures stability for general advection equations. You might gain some leverage by adjusting the spatial and temporal discretization, though.

[david-zwicker]

I moved the issue to a discussion, where it fits better.