Status of CBC.PDESys?

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.
>

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