Error in form of advection type equation

Question

I am trying to solve the following set of equations:

$$\frac{\partial P_x}{\partial t} + \frac{\partial D}{\partial x} = 0$$

$$\frac{\partial P_y}{\partial t} + \frac{\partial D}{\partial y} = 0$$

Where P_x is the partial derivative of a scalar function P with respect to x. For P with some boundary conditions, where D is known.

The form I use with backward Euler is:

F3 = (1.0/dt)*dot(grad(p_)-grad(p_n), v)*dx + dot(grad(D), v)*dx
a3 = lhs(F3)
L3 = rhs(F3)

A3 = assemble(a3)

When I try running this I get and error in the bcp line:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.981e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
...
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 62 (Invalid argument).

What does the error mean? Am I setting the problem up correctly? The boundary conditions I would like to add are $P_x = 0$ at $x=0$ and $x=2$ and $P_y = 0$ at $y=0$.

The code is here:

from fenics import *

dt = 0.01 # time step size
nElems = 32

# Create mesh and define function spaces
mesh = UnitSquareMesh(nElems, nElems)

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define trial and test functions
v = TestFunction(V)


# Define functions for solutions at previous and current time steps
p_n = Function(Q)
p_  = Function(Q)
D   = Function(Q)

# Define expressions used in variational forms
k   = Constant(dt)


F3 = (1.0/dt)*dot(grad(p_)-grad(p_n), v)*dx + dot(grad(D), v)*dx
a3 = lhs(F3)
L3 = rhs(F3)

A3 = assemble(a3)

b3 = assemble(L3)

solve(A3, p_.vector(), b3)

Answer 1

I got it working by solving the two equations individually treating P_x and P_y as two scalar fields, this also allows me to add the boundary conditions (bcs) as just Dirichlet bcs. For the D partial derivatives I used D.dx(0) and D.dx(1).

I have to use the gradient of P in a ufl expression afterwards.
I found I can use as_vector([P_x, P_y]) and it seems to be working instead of grad(P).
Is this correct usage?

I would still like to know if I can use the original expression I have in the question so that I don't have to find P from grad(P) at the end of the calculation

Comments

I have solved advection diffusion problem like this and this is working hope it will help you.

from dolfin import *

mesh = IntervalMesh(20,0,10)
#plot(mesh, interactive = True)

V = FunctionSpace(mesh, 'CG',2)
def boundary1(x):
    return x[0]==0 
def boundary2(x):
    return x[0]==10

bc1 = DirichletBC(V, Constant(0.0), boundary1)
bc2 = DirichletBC(V, Constant(1.0), boundary2)

bc = [bc1, bc2]

u = TrialFunction(V)
v = TestFunction(V)
uv = as_vector((v,))
du = Function(V)
print shape(v)
print shape(grad(u))
a = ((5.0)*(inner(nabla_grad(v),nabla_grad(u))))*dx +(2.0*(inner(uv, nabla_grad(u))))*dx

solve(lhs(a)==rhs(a), du, bc)
plot(du, interactive = True)