running flow past a dolphin example

Hi,

Sorry about that. I thought it was fixed already. The constant nu is
missing from the PDESubSystem and there is even another typo in the form.
Try

from cbc.pdesys import *

# Set up problem by loading mesh from file
mesh = Mesh('dolfin-2.xml.gz')

# problem_parameters are defined in Problem.py
problem_parameters['time_integration'] = "Transient" # default='Steady'
problem = Problem(mesh, problem_parameters)

# Set up first PDESystem
solver_parameters['space']['u'] = VectorFunctionSpace #
default=FunctionSpace
solver_parameters['degree']['u'] = 2 # default=1
NStokes = PDESystem([['u', 'p']], problem, solver_parameters)

# Use a constant forcing field to drive the flow from right to left
NStokes.f = Constant((-1., 0.))

# No-slip boundary condition for velocity on the dolfin
dolfin = AutoSubDomain(lambda x, on_boundary: on_boundary and not
                       (near(x[0], 0) or near(x[0], 1.) or near(x[1], 0.)
or near(x[1], 1.)))

bc = [DirichletBC(NStokes.V['up'].sub(0), Constant((0.0, 0.0)), dolfin)]

# Set up variational form.
# u_, u_1 are the solution Functions at time steps N and N-1.
# v_u/v_p are the TestFunctions for velocity/pressure in the
MixedFunctionSpace for u and p

NStokes.nu = Constant(0.01)
class NavierStokes(PDESubSystem):
    def form(self, u, v_u, u_, u_1, p, v_p, nu, dt, f, **kwargs):
        U = 0.5*(u + u_1)
        return (1./dt)*inner(u - u_1, v_u)*dx + \
               inner(dot(u_1, nabla_grad(u_1)), v_u)*dx + \
               nu*inner(grad(U), grad(v_u))*dx - \
               inner(p, div(v_u))*dx + inner(div(U), v_p)*dx - \
               inner(f, v_u)*dx

NStokes.pdesubsystems['up'] = NavierStokes(vars(NStokes), ['u', 'p'],
bcs=bc,
                                           reassemble_lhs=False)

# Integrate the solution from t=0 to t=0.5
problem.prm['T'] = 0.5
problem.solve()

# Define a new nonlinear PDESystem for a scalar c
scalar = PDESystem([['c']], problem, solver_parameters)

class Scalar(PDESubSystem):
    def form(self, c, v_c, c_, c_1, U_, dt, nu, **kwargs):
        C = 0.5*(c + c_1)
        return (1./dt)*inner(c - c_1, v_c)*dx + \
                inner(dot(U_, grad(C)), v_c)*dx + \
                nu*(1. + c_**2)*inner(grad(C), grad(v_c))*dx
                # Note nonlinearity in c_ (above)

bcc = [DirichletBC(scalar.V['c'], Constant(1.0), dolfin)]

scalar.U_ = 0.5*(NStokes.u_ + NStokes.u_1) # The Scalar form uses the
velocity
scalar.nu = NStokes.nu
csub1 = Scalar(vars(scalar), ['c'], bcs=bcc, max_inner_iter=5) # Iterate on
c_
scalar.pdesubsystems['c'] = csub1

# Integrate both PDESystems from t=0.5 to t=1.0 using Picard
# iterations on each time step
problem.prm['T'] = 1.0
problem.solve()

# Switch to using the Newton method for the nonlinear variational form
# With these calls we replace c by c_ in the Scalar form and compute the
Jacobian wrt c_
csub1.prm['iteration_type'] = 'Newton'
csub1.define()

# Integrate both PDESystems from T=1.0 to T=1.5 using Newton
# iterations on each time step for the scalar
problem.prm['T'] = 1.5
problem.solve()

Best regards

Mikael

On 5 November 2012 16:51, nicusor <email address hidden>wrote:

> New question #213341 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/213341
>
> I get an error:
>
> Creating new work vector for up
> Adding PDESubSystem: NavierStokes
> Adding ['linear_solver']['up'] = lu to pdesubsystem up
> Adding ['iteration_type'] = Picard to pdesubsystem up
> Traceback (most recent call last):
> File "dolphin.py", line 37, in <module>
> NStokes.pdesubsystems['up'] = NavierStokes(vars(NStokes), ['u', 'p'],
> bcs=bc, reassemble_lhs=False)
> File "/usr/lib/python2.7/cbc/pdesys/PDESubSystems.py", line 331, in
> __init__
> self.define()
> File "/usr/lib/python2.7/cbc/pdesys/PDESubSystems.py", line 337, in
> define
> self.get_form(form_args)
> File "/usr/lib/python2.7/cbc/pdesys/PDESubSystems.py", line 301, in
> get_form
> F = self.form(**form_args)
> TypeError: form() takes exactly 10 arguments (9 given)
>
> Any idea how to solve this.
>
> Thanks
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>