Thank you, that conditional worked. I have realized that I should in fact have trial functions, (equations are linear in (n+1) terms), I have been mistakenly following example scripts for nonlinear problems. I have changed v, p, phi, mu = split(u) rho = Function(P1) to v, p, phi, mu = TrialFunctions(W) rho = TrialFunction(P1) and I now I do have a bilinear form so I changed a = derivative(a1 + a2 + a3 + a4, u) back to a = a1 + a2 + a3 + a4 and b1 = (1/dt)*eta*rho*dx b = derivative(b1, rho) to b = (1/dt)*eta*rho*dx Everything else I kept unchanged. Now when I try to run the code I get the following: Traceback (most recent call last): File "navier_stokes_CH2.py", line 115, in <module> solve(a == L, u, bcs=[bc0, bc1, bc2, bc3], solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'}) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 145, in _solve_varproblem options_prefix=options_prefix) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 252, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 126, in __init__ ctx = solving_utils._SNESContext(problem) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 106, in __init__ for problem in problems) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 106, in <genexpr> for problem in problems) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", line 66, in assemble inverse=inverse, nest=nest) File "<decorator-gen-296>", line 2, in _assemble File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/assemble.py", line 100, in _assemble inverse=inverse) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", line 243, in compile_form number_map).kernels File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line 203, in __new__ obj = make_obj() File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/caching.py", line 193, in make_obj obj.__init__(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/tsfc_interface.py", line 165, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/tsfc/driver.py", line 42, in compile_form do_estimate_degrees=True) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/compute_form_data.py", line 362, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 143, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 137, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/corealg/map_dag.py", line 84, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/ufl/algorithms/check_arities.py", line 55, in product raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x.number(), x)) ufl.algorithms.check_arities.ArityMismatch: Multiplying expressions with overlapping form argument number 1, argument is v_1. Exception AttributeError: "'LinearVariationalSolver' object has no attribute '_parameters'" in <bound method LinearVariationalSolver.__del__ of <firedrake.variational_solver.LinearVariationalSolver object at 0x457c750>> ignored ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Miklós Homolya <m.homolya14@imperial.ac.uk> Sent: 06 April 2016 14:50 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Integral is missing an integration domain. 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<mailto:firedrake-bounces@imperial.ac.uk> <firedrake-bounces@imperial.ac.uk><mailto:firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk><mailto:lawrence.mitchell@imperial.ac.uk> Sent: 04 April 2016 15:37 To: firedrake@imperial.ac.uk<mailto: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><mailto: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<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake