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

        FPw = inner(grad(Pw),grad(v_Pw))*dx
        FSw = inner(grad(Pw),grad(v_Sw))*dx + Sw/dt*v_Sw*dx - Sw_1/dt*v_Sw*dx
        return FPw + FSw

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

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

# add the pde subsystem
coupled.add_pdesubsystem(SS2Ph, ['Pw', 'Sw'], bcs=bcs)

class InitialConditions(Expression):
    def __init__(self):
        print "initialized initial conditions"
    def eval(self, values, x):
        values[0] = 10.0
        values[1] = 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:
Last reply:

This question was reopened

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:
> https://answers.launchpad.net/cbcpdesys/+question/211982
>
> 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):
>
> FPw = inner(grad(Pw),grad(v_Pw))*dx
> FSw = inner(grad(Pw),grad(v_Sw))*dx + Sw/dt*v_Sw*dx - Sw_1/dt*v_Sw*dx
> return FPw + FSw
>
>
> # setup the BCs
> def left(x):
> return x[0] < DOLFIN_EPS
> def right(x):
> return x[0] > 1.0-DOLFIN_EPS
> def top(x):
> return x[1] > 1.0-DOLFIN_EPS
> def bottom(x):
> return x[1] < DOLFIN_EPS
>
> bcs = [DirichletBC(coupled.V['Pw'], (15.0), left),
> DirichletBC(coupled.V['Pw'], (5.0), right)]
>
>
> # add the pde subsystem
> coupled.add_pdesubsystem(SS2Ph, ['Pw', 'Sw'], bcs=bcs)
>
> class InitialConditions(Expression):
> def __init__(self):
> print "initialized initial conditions"
> def eval(self, values, x):
> values[0] = 10.0
> values[1] = 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()
>
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Yaakoub El Khamra (yelkhamra) said : #2

Well for one thing I would expect the boundary on the left to be 15 (x[0]=0), the right to be 5(x[0]=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[0]=1 and x[1]> 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 : #3

Sorry, accidentally hit problem solved.

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

> Question #211982 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/211982
>
> Status: Answered => Solved
>
> 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]=0), the right to be 5(x[0]=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[0]=1 and x[1]> 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.
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Yaakoub El Khamra (yelkhamra) said : #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 : #6

sorry, hit problem solved, again.

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:
> https://answers.launchpad.net/cbcpdesys/+question/211982
>
> Status: Answered => Solved
>
> 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?
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Yaakoub El Khamra (yelkhamra) said : #8

Sent. Thank you very much in advance.

Yaakoub El Khamra (yelkhamra) said : #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