Cavity driven flow : Form (<empty Form>) seems to be zero: cannot compile it.

On Thu, 02 May 2013 11:56:12 -0000
imranal <email address hidden> wrote:
> New question #228011 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/228011
>
> I have implemented a NS solver for a cavity driven flow. The solver
> triggers the following error message (the code comes right after) :
>
> >> python cavity_driven.py
>
> No Jacobian form specified for nonlinear variational problem.
> Differentiating residual form F to obtain Jacobian J = F'.
> Calling FFC just-in-time (JIT) compiler, this may take some time.
> Form (<empty Form>) seems to be zero: cannot compile it.
> Traceback (most recent call last):
> File "cavity_driven.py", line 67, in <module>
> solve(F ==
> Constant((0,0)),up_,bcs,solver_parameters={"newton_solver" :{"maximum_iterations" :
> 30}}) File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py",
> line 250, in solve _solve_varproblem(*args, **kwargs) File
> "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 287,
> in _solve_varproblem
> form_compiler_parameters=form_compiler_parameters) File
> "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 118,
> in __init__ J = Form(J,
> form_compiler_parameters=form_compiler_parameters) File
> "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 54, in
> __init__ common_cell) File
> "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line
> 66, in mpi_jit return local_jit(*args, **kwargs) File
> "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line
> 154, in jit return jit_compile(form, parameters=p,
> common_cell=common_cell) File
> "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 73, in
> jit return jit_form(ufl_object, parameters, common_cell) File
> "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 130, in
> jit_form common_cell=common_cell) File
> "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 150, in
> compile_form analysis = analyze_forms(forms, object_names,
> parameters, common_cell) File
> "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 64, in
> analyze_forms common_cell) for form in forms) File
> "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 64, in
> <genexpr> common_cell) for form in forms) File
> "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 134, in
> _analyze_form "Form (%s) seems to be zero: cannot compile it." %
> str(form)) File "/usr/lib/python2.7/dist-packages/ffc/log.py", line
> 46, in ffc_assert condition or error(*message) File "<string>", line
> 1, in <lambda> File "/usr/lib/python2.7/dist-packages/ufl/log.py",
> line 148, in error raise
> self._exception_type(self._format_raw(*message)) Exception: Form
> (<empty Form>) seems to be zero: cannot compile it.
> ___________________________________________________________________
>
> The code :
>
> """
> Cavity driven flow for a square mesh.
>
> u_top = (1.0,0.0)
> _>_>_>_>_>_>_>_>_>_
> | |
> | |
> | |
> | |
> | |
> |____________________|
> No-slip at walls.
>
>
> In this code we attempt to accomplish the following :
> 1) Solve the coupled problem for velocity and pressure for Re = 100.
>
> """
> from dolfin import *
>
> # Physical parameter
> mu = Constant(100) # dynamic viscosity
> # Dimensionless number
> Re = Constant(100) # Reynold's number
>
> # Velocities at the boundaries
> u_top_lid = Constant((1, 0)) # Cavity driven
> u_bottom_wall = Constant((0, 0)) # No-slip
> u_side_walls = Constant((0, 0)) # No-slip
>
> def bnd_Rside(x,on_boundary):
> return on_boundary and near(x[0],1)
> def bnd_Lside(x,on_boundary):
> return on_boundary and near(x[0],0)
> def bnd_top(x,on_boundary):
> return on_boundary and near(x[1],1)
> def bnd_bottom(x,on_boundary):
> return on_boundary and near(x[1],0)
>
> # Define the mesh and the mixed space
> mesh = UnitSquare(50,50)
> V = VectorFunctionSpace(mesh,"CG",2)
> Q = FunctionSpace(mesh,"CG",1)
> VQ = V*Q
>
> # The boundary conditions
> bc1 = DirichletBC(VQ.sub(0), u_bottom_wall,bnd_bottom)
> bc2 = DirichletBC(VQ.sub(0), u_side_walls, bnd_Rside)
> bc3 = DirichletBC(VQ.sub(0), u_side_walls, bnd_Lside)
> bc4 = DirichletBC(VQ.sub(0), u_top_lid, bnd_top)
> bcs = [bc1,bc2,bc3,bc4]
>
> # Make an initial guess
> u = interpolate(Expression(("x[0]","0")),V)
> p = interpolate(Constant(1.0),Q)

Instead of this do:

w = Function(VQ)
u, p = split(w)

Your version does not work because your F does not depend on function
from VQ space, so when it is differetiated by solution coefficient from
VQ space, derivative is zero.

>
> # The weak formulation
> v,q = TestFunctions(VQ)
> F = mu*inner(grad(v),grad(u))*dx - inner(div(v),p)*dx -
> inner(q,div(u))*dx -\ Re*inner(grad(u)*u, v)*dx
>
> # Solve and store in new functions
> up_ = Function(VQ)
> solve(F ==
> Constant((0,0)),up_,bcs,solver_parameters={"newton_solver" :{"maximum_iterations" :
> 30}})

Rank of F is 0 not 1:

solve(F==0,w,bcs,solver_parameters={"newton_solver" :{"maximum_iterations":30}})

>
> ___________________________________________________
>
> The last triggers the error. What is the error message actually
> telling me ? Is the form F empty ? There is also something else which
> seems strange. If i change F slightly by interchanging the position
> of the nonlinear term :
>
> inner(grad(u)*u, v)*dx

This is correct and means same as
  inner(dot(grad(u), u), v)*dx
or
  inner(dot(u, nabla_grad(u)), v)*dx

>
> as
>
> inner(u*grad(u), v)*dx

This is not correct. UFL does not allow it. Read UFL chapter in FEniCS
book.

>
> I get a new error :
>
> Product can only represent products of scalars.
> Traceback (most recent call last):
> File "cavity_driven.py", line 63, in <module>
> Re*inner(u*grad(u), v)*dx
> File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line
> 120, in _mul return _mult(self, o)
> File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line
> 93, in _mult p = Product(a, b[ii])
> File "/usr/lib/python2.7/dist-packages/ufl/algebra.py", line 177,
> in __new__ error("Product can only represent products of scalars.")
> File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 148, in
> error raise self._exception_type(self._format_raw(*message))
> ufl.log.UFLException: Product can only represent products of scalars.
>
> _________________________________________________________________
>
> I used inner(grad(u)*u, v)*dx initially as this was used in the NS
> demo.
>
>
>

>> # Make an initial guess
>> u = interpolate(Expression(("x[0]","0")),V)
>> p = interpolate(Constant(1.0),Q)
>
>Instead of this do:
>
>w = Function(VQ)
>u, p = split(w)
>
>Your version does not work because your F does not depend on function
>from VQ space, so when it is differetiated by solution coefficient from
>VQ space, derivative is zero.

Then how can I assign spesific guesses for the respective variables ?
Correcting the code as you suggested still yielded the same error message.

>> inner(grad(u)*u, v)*dx
>
>This is correct and means same as
 > inner(dot(grad(u), u), v)*dx
> or
> inner(dot(u, nabla_grad(u)), v)*dx
>
>>
>> as
>>
>> inner(u*grad(u), v)*dx
>
>This is not correct. UFL does not allow it. Read UFL chapter in FEniCS book.

Will do, thanks for reference!

Jan, I am currently running FEniCS 1.0. Specifying the Jacobian did the trick!

On Thu, 2 May 2013 16:03:22 +0200
Jan Blechta <email address hidden> wrote:
> On Thu, 02 May 2013 13:36:00 -0000
> imranal <email address hidden> wrote:
> > Question #228011 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/228011
> >
> > Status: Answered => Open
> >
> > imranal is still having a problem:
> > >> # Make an initial guess
> > >> u = interpolate(Expression(("x[0]","0")),V)
> > >> p = interpolate(Constant(1.0),Q)
> > >
> > >Instead of this do:
> > >
> > >w = Function(VQ)
> > >u, p = split(w)
> > >
> > >Your version does not work because your F does not depend on
> > >function from VQ space, so when it is differetiated by solution
> > >coefficient from VQ space, derivative is zero.
> >
> > Then how can I assign spesific guesses for the respective
> > variables ?
>
> Assignmennt, projection or interpolation to subfunctions is not still
> satisfactorily solved in DOLFIN. There's a blueprint on this. For now
> you can write forms to do a projection like:
>
> w_test, w_trial = TestFunction(VQ), TrialFunction(VQ)
> a_proj = inner(w_trial, w_test)*dx
> L_proj = inner(U, v)*dx + P*q*dx
> solve(a_proj == L_proj)

Of course, here should be:
solve(a_proj == L_proj, w)

>
> where U, P are your initializers - Expression, Constants, Functions
> (of arbitrary function space or mesh, it suffices that it can be
> evaluated on VQ's mesh)
>
> For more help search 'assignment to subfunction' on launchpad
> questions.
>
> > Correcting the code as you suggested still yielded the same error
> > message.
>
> It works form me. What's your FEniCS version? Maybe it is obligatory
> to specify Jacobian to solve function in earlier versions. Try:
>
> J = derivative(F, w)
> solve(F==0, w, bcs, J=J, ...)
>
> >
> > >> inner(grad(u)*u, v)*dx
> > >
> > >This is correct and means same as
> > > inner(dot(grad(u), u), v)*dx
> > > or
> > > inner(dot(u, nabla_grad(u)), v)*dx
> > >
> > >>
> > >> as
> > >>
> > >> inner(u*grad(u), v)*dx
> > >
> > >This is not correct. UFL does not allow it. Read UFL chapter in
> > >FEniCS book.
> >
> > Will do, thanks for reference!
> >
>

By the way, how do you respond to my messages ? I can¨t seem to respond to yours ^^

Thanks Jan Blechta, that solved my question.

On Thu, 02 May 2013 14:21:17 -0000
imranal <email address hidden> wrote:
> Question #228011 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/228011
>
> imranal posted a new comment:
> By the way, how do you respond to my messages ? I can¨t seem to
> respond to yours ^^
>

I'm subscribed to DOLFIN answers mailing list through DOLFIN team
membership. But don't deal woth it because questions are currently
migrating to stackexchange.

Den May 2, 2013 kl. 4:11 PM skrev Jan Blechta:

> Question #228011 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/228011
>
> Status: Open => Answered
>
> Jan Blechta proposed the following answer:
> On Thu, 02 May 2013 13:36:00 -0000
> imranal <email address hidden> wrote:
>> Question #228011 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/228011
>>
>> Status: Answered => Open
>>
>> imranal is still having a problem:
>>>> # Make an initial guess
>>>> u = interpolate(Expression(("x[0]","0")),V)
>>>> p = interpolate(Constant(1.0),Q)
>>>
>>> Instead of this do:
>>>
>>> w = Function(VQ)
>>> u, p = split(w)
>>>
>>> Your version does not work because your F does not depend on function
>>> from VQ space, so when it is differetiated by solution coefficient
>>> from VQ space, derivative is zero.
>>
>> Then how can I assign spesific guesses for the respective variables ?
>
> Assignmennt, projection or interpolation to subfunctions is not still
> satisfactorily solved in DOLFIN. There's a blueprint on this. For now
> you can write forms to do a projection like:
>
> w_test, w_trial = TestFunction(VQ), TrialFunction(VQ)
> a_proj = inner(w_trial, w_test)*dx
> L_proj = inner(U, v)*dx + P*q*dx
> solve(a_proj == L_proj)

Hi Imran,

Happy to see that you're working on your mandatory assignment:-)

BTW, you can also initialize like this

w = interpolate(Expression(("x[0]", "0", "1."), element=VQ.ufl.element()), VQ)

If you're happy with constants you could use Constant instead of Expression above. You can also do

VQ.sub(0).sub(0).dofmap().set(w.vector(), 1)
VQ.sub(0).sub(1).dofmap().set(w.vector(), 2)
VQ.sub(1).dofmap().set(w.vector(), 3)

Best regards

Mikael

>
> where U, P are your initializers - Expression, Constants, Functions
> (of arbitrary function space or mesh, it suffices that it can be
> evaluated on VQ's mesh)
>
> For more help search 'assignment to subfunction' on launchpad questions.
>
>> Correcting the code as you suggested still yielded the same error
>> message.
>
> It works form me. What's your FEniCS version? Maybe it is obligatory to
> specify Jacobian to solve function in earlier versions. Try:
>
> J = derivative(F, w)
> solve(F==0, w, bcs, J=J, ...)
>
>>
>>>> inner(grad(u)*u, v)*dx
>>>
>>> This is correct and means same as
>>> inner(dot(grad(u), u), v)*dx
>>> or
>>> inner(dot(u, nabla_grad(u)), v)*dx
>>>
>>>>
>>>> as
>>>>
>>>> inner(u*grad(u), v)*dx
>>>
>>> This is not correct. UFL does not allow it. Read UFL chapter in
>>> FEniCS book.
>>
>> Will do, thanks for reference!
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Great! So when is the deadline tomorrow ^^ ?