Status of CBC.PDESys?

Asked by Yaakoub El Khamra on 2012-08-08

So what's the status of the CBC.PDESys package? Is it active? Is it being maintained? I tried to run the flow past a heated dolphin demo here:
http://fenicsproject.org/featured/2011/pdesys.html
and got these errors:
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 "/home/yye00/Dropbox/Code/FEniCS_Code/Python/Sovereign/FlowPastDolfin.py", line 45, in <module>
    reassemble_lhs=False)
  File "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py", line 331, in __init__
    self.define()
  File "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py", line 337, in define
    self.get_form(form_args)
  File "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py", line 301, in get_form
    F = self.form(**form_args)
TypeError: form() takes exactly 10 arguments (9 given)

Any suggestions? Are there alternatives to the CBC.PDESys package that allow easy specification of large systems of coupled non-linear PDE's with ease?

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Yaakoub El Khamra
Solved:
Last query:
Last reply:

Hi Yaakoub,

The CBC.PDESys package is still active. Sorry about the dolphin demo, the
code on the webpage is not entirely up to date and there even seem to be
some typos. I attach a slighly modified and fixed script that runs the flow
past dolphin case, see if you can make that work.

I don't think there are any better option (or at least not simpler) than
using this package for large systems of nonlinear equations. Have a look at
the demos, they should be up to date.

Best regards

Mikael

from cbc.pdesys import *

# Set up problem by loading mesh from file
mesh =
Mesh('/home/mikael/Fenics/dorsal/precise/unstable/share/dolfin/data/meshes/dolfin-2.xml.gz')
mesh = refine(mesh)

# 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.))
NStokes.nu = Constant(0.005) ## This was missing

# 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

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
solver_parameters['familyname'] = 'Scalar'
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
def update(self):
    plot(self.pdesystems['Scalar'].c_)
Problem.update = update

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()

On 8 August 2012 20:45, Yaakoub El Khamra <
<email address hidden>> wrote:

> New question #205330 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/205330
>
>
> So what's the status of the CBC.PDESys package? Is it active? Is it being
> maintained? I tried to run the flow past a heated dolphin demo here:
> http://fenicsproject.org/featured/2011/pdesys.html
> and got these errors:
> 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
> "/home/yye00/Dropbox/Code/FEniCS_Code/Python/Sovereign/FlowPastDolfin.py",
> line 45, in <module>
> reassemble_lhs=False)
> File
> "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py",
> line 331, in __init__
> self.define()
> File
> "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py",
> line 337, in define
> self.get_form(form_args)
> File
> "/home/yye00/Work/FEniCS/lib64/python2.7/site-packages/cbc/pdesys/PDESubSystems.py",
> line 301, in get_form
> F = self.form(**form_args)
> TypeError: form() takes exactly 10 arguments (9 given)
>
> Any suggestions? Are there alternatives to the CBC.PDESys package that
> allow easy specification of large systems of coupled non-linear PDE's with
> ease?
>
>
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Yaakoub El Khamra (yelkhamra) said : #2

Thank you very, very much! The demos work and I am going over them right now.