Dear Fryderyk,

You are trying to do a symbolic operation on numeric value, and, of course, it doesn't work as you expect it.

f = Function(P1).interpolate(FreeEnergyPotential())
This construct f as the numerical values in the given function space.
dfdphi = diff(f, phi0)
This tries to symbolically differentiate f according to phi0, but f does not depend on phi0, so you get zero, hence the error message.

You should define f symbolically:

f = conditional(phi0 > 1, (phi0 - 1)**2, conditional(phi0 < -1, (phi0 + 1)**2, 0.25*(phi0**2 - 1)**2))

This way the symbolic differentiation will work correctly.

On 06/04/16 13:47, Fryderyk Wilczynski wrote:
Dear all,

I want to define a piecewise expression, the derivative of which will be used on the RHS in my weak for and needs to be updated every timestep.
So where I previously had: 

phi = variable(phi)
f    = 0.25*(phi**2 - 1)**2
dfdphi = diff(f, phi)


I want to define f as follows: 
f =      (phi +1)^2                 if     phi < -1 
      0.25* ( phi^2 - 1 )^2     if   -1 < phi < 1 
      (phi - 1)^2                     if     phi > 1 


I have tried the following

phi0 = variable(phi0)
class FreeEnergyPotential(Expression):
def eval(self, values, phi0):
if phi0 > 1:
values[:] = (phi0-1)**2
elif phi0 < -1:
values[:] = (phi0+1)**2
else:
values[:] = 0.25*(phi0**2 -1)**2
f= Function(P1).interpolate(FreeEnergyPotential())
dfdphi = diff(f, phi0)

But that didn't work and I get the familiar exception:
This integral is missing an integration domain.
Traceback (most recent call last):
  File "navier_stokes_CH.py", line 96, in <module>
    L2 = - (beta/Cn)*phi0*psi*dx + (1/Cn)*dfdphi*psi*dx 
  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/measure.py", line 398, in __rmul__
    error("This integral is missing an integration domain.")
  File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.

How can I implement it correctly?
Full code attached.

Thanks,
Fryderyk


________________________________________
From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
Sent: 04 April 2016 15:37
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] Integral is missing an integration domain.


On 4 Apr 2016, at 14:59, Fryderyk Wilczynski <scfw@leeds.ac.uk> wrote:

Thank you for spotting my mistakes!

I have implemented the suggested corrections and added the no-slip boundary conditions.

...


New code:

rho_hat = Function(P1)

Here you make a function rho_hat which is initialised to zero.


a2 = psi*mu*dx - (beta/Cn)*phi*psi*dx - Cn*inner(grad(phi),grad(psi))*dx - (1/rho_hat)*drho_hat*psi*p*dx

This term has a 1/rho_hat -> 1/0.0 -> NAN


      # compute density using algebraic expression rho_hat = 0.5*( (1+phi) + (rho2/rho1)*(1-phi))
      rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))

...

This just assigns the name rho_hat the expression on the RHS, it does not magically go back and substitute it in to the expression for a2 above.

My suggestion was to define the expression before using it when defining the variational forms:

rho_hat = 0.5*((1 + phi0) + (rho2/rho1)*(1 - phi0))


a1 = ...
a2 = ...
a3 = ...

Then no need to do anything in the while loop.

      # Solve the system of equations
      solve(a == L, u, bcs=[bc0, bc1, bc2, bc3]) #
      # update rho
      rho -= dt*div(rho*v0)

This line doesn't do what you think it does.  You cannot compute div(rho*v0) point wise.

This is your update for equation (10), right?

rho^{n+1} = rho^n - dt div (rho^n v^n)

You will need to solve this equation weakly in the normal way.  Else wise (although I am not sure if this will destroy the conservation properties of the scheme) you could interpolate the expression:

tmp_rho = Function(P1)

while ...:

    ...
    tmp_rho.interpolate(rho - dt*div(rho*v0))
    rho.assign(tmp_rho)


Thanks,

Lawrence



_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake