# problem with a coupled system

Asked by Yaakoub El Khamra on 2012-10-22

I am very sorry to have to ask this but can you please have a look at my code? The BC's are completely off for this coupled system, nothing like what I would expect it to be.

from cbc.pdesys import *

mesh = UnitSquare(20, 20)

# Change desired items in the problem_parameters dict from cbc.pdesys
problem_parameters['time_integration'] = "Transient" # change this to transient
problem = Problem(mesh, problem_parameters)

#solver_parameters['family'] = 'DG'
solver_parameters['familyname'] = 'Scalar'
solver_parameters['iteration_type'] = 'Picard'
coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

class SS2Ph(PDESubSystem):
def form(self, Pw, v_Pw, Pw_, Sw, Sw_, Sw_1, v_Sw, dt, **kwargs):

return FPw + FSw

# setup the BCs
def left(x):
return x < DOLFIN_EPS
def right(x):
return x > 1.0-DOLFIN_EPS
def top(x):
return x > 1.0-DOLFIN_EPS
def bottom(x):
return x < DOLFIN_EPS

bcs = [DirichletBC(coupled.V['Pw'], (15.0), left),
DirichletBC(coupled.V['Pw'], (5.0), right)]

class InitialConditions(Expression):
def __init__(self):
print "initialized initial conditions"
def eval(self, values, x):
values = 10.0
values = 0.5
def value_shape(self):
return (2,)

def update(self):
plot(self.pdesystems['Scalar'].Sw_, interactive=True)
plot(self.pdesystems['Scalar'].Pw_, interactive=True)

Problem.update = update

problem.q0 = {'PwSw': InitialConditions()}

problem.prm['T'] = 3.0
problem.initialize(coupled)
problem.solve()

## 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:

## This question was reopened

 Mikael Mortensen (mikael-mortensen) said on 2012-10-22: #1

Hi Yakoub,

Not quite sure what you mean, it seems to work for me? What's the problem? Think you should use DirichletBC(coupled.V['PwSw'].sub(0), (15.0), left) though, but it doesn't make a difference here.

Mikael

Den Oct 22, 2012 kl. 6:51 PM skrev Yaakoub El Khamra:

> New question #211982 on CBC.PDESys:
>
> I am very sorry to have to ask this but can you please have a look at my code? The BC's are completely off for this coupled system, nothing like what I would expect it to be.
>
> from cbc.pdesys import *
>
> mesh = UnitSquare(20, 20)
>
> # Change desired items in the problem_parameters dict from cbc.pdesys
> problem_parameters['time_integration'] = "Transient" # change this to transient
> problem = Problem(mesh, problem_parameters)
>
> #solver_parameters['family'] = 'DG'
> solver_parameters['familyname'] = 'Scalar'
> solver_parameters['iteration_type'] = 'Picard'
> coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.
>
>
> class SS2Ph(PDESubSystem):
> def form(self, Pw, v_Pw, Pw_, Sw, Sw_, Sw_1, v_Sw, dt, **kwargs):
>
> return FPw + FSw
>
>
> # setup the BCs
> def left(x):
> return x < DOLFIN_EPS
> def right(x):
> return x > 1.0-DOLFIN_EPS
> def top(x):
> return x > 1.0-DOLFIN_EPS
> def bottom(x):
> return x < DOLFIN_EPS
>
> bcs = [DirichletBC(coupled.V['Pw'], (15.0), left),
> DirichletBC(coupled.V['Pw'], (5.0), right)]
>
>
> # add the pde subsystem
>
> class InitialConditions(Expression):
> def __init__(self):
> print "initialized initial conditions"
> def eval(self, values, x):
> values = 10.0
> values = 0.5
> def value_shape(self):
> return (2,)
>
> def update(self):
> plot(self.pdesystems['Scalar'].Sw_, interactive=True)
> plot(self.pdesystems['Scalar'].Pw_, interactive=True)
>
> Problem.update = update
>
> problem.q0 = {'PwSw': InitialConditions()}
>
> problem.prm['T'] = 3.0
> problem.initialize(coupled)
> problem.solve()
>
>
> --
> contact for CBC.PDESys.

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #2

Well for one thing I would expect the boundary on the left to be 15 (x=0), the right to be 5(x=1). It looks as if the bottom right corner has a value of 5, i.e. as if the BC is rotated by -45 degrees. x=1 and x> 0.5 seems to have a value of 15 (should be 5)

This would be easier if I can attach a screen shot, but I cannot find a way to do it.

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #3

Sorry, accidentally hit problem solved.

 Mikael Mortensen (mikael-mortensen) said on 2012-10-22: #4

Den Oct 22, 2012 kl. 7:21 PM skrev Yaakoub El Khamra:

> Question #211982 on CBC.PDESys changed:
>
>
> Yaakoub El Khamra confirmed that the question is solved:
>
> Well for one thing I would expect the boundary on the left to be 15 (x=0), the right to be 5(x=1). It looks as if the bottom right corner has a value of 5, i.e. as if the BC is rotated by -45 degrees. x=1 and x> 0.5 seems to have a value of 15 (should be 5)

It is 15 on the left and 5 on the right, at least in my plot. Are you looking at the second plot? There is no bc on Sw.

Mikael

>
> This would be easier if I can attach a screen shot, but I cannot find a
> way to do it.
>
> --
> contact for CBC.PDESys.

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #5

The first plot is very weird for me. The second plot is what I am referring to. The second plot should be for Pw, and that's where I apply the BC. I think.

Should I just put the plot results online and paste a link or is there a way to attach a file in launchpad?

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #6

sorry, hit problem solved, again.

 Mikael Mortensen (mikael-mortensen) said on 2012-10-22: #7

Just send it to <email address hidden>

Mikael

Den Oct 22, 2012 kl. 8:45 PM skrev Yaakoub El Khamra:

> Question #211982 on CBC.PDESys changed:
>
>
> Yaakoub El Khamra confirmed that the question is solved:
>
> The first plot is very weird for me. The second plot is what I am
> referring to. The second plot should be for Pw, and that's where I apply
> the BC. I think.
>
> Should I just put the plot results online and paste a link or is there a
> way to attach a file in launchpad?
>
> --
> contact for CBC.PDESys.

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #8

Sent. Thank you very much in advance.

 Yaakoub El Khamra (yelkhamra) said on 2012-10-22: #9

I just sent you an email sir.
using: DirichletBC(coupled.V['PwSw'].sub(0), (15.0), left) and DirichletBC(coupled.V['PwSw'].sub(0), (5.0), right) gives the correct answer. Thank you very much for your help with this.

Regards
Yaakoub