CBC.PDESys

Initial conditions

Asked by Yaakoub El Khamra

Minor problem with initial conditions: I am supposed to use the q0 dictionary but I am not sure if I am doing it correctly. The initial conditions __init__ is called, I know that much, but the value of u_1 is never set. I am sure I am doing this wrong as it does work in the CahnHilliard example. Any idea what I am doing wrong?

from cbc.pdesys import *

mesh = UnitSquare(10, 10)
# setup the problem parameters
problem_parameters['time_integration'] = "Transient" # default='Steady'
problem = Problem(mesh, problem_parameters)
# setup the solver parameters
solver_parameters['familyname'] = 'Scalar'
solver_parameters['iteration_type'] = 'Newton'
# setup the pdesystem
poisson = PDESystem([['u']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

class Poisson(PDESubSystem):
    def form(self, u, v_u, u_, u_1, dt, **kwargs): # v_u is the TestFunction
        return inner(grad(u), grad(v_u))*dx + 1.0/dt*u*v_u*dx -1.0/dt*u_1*v_u*dx

bcs = DirichletBC(poisson.V['u'], (0.0), DomainBoundary())
# add the pdesubsystem
poisson.pdesubsystems['u'] = Poisson(vars(poisson), ['u'], bcs=[bcs])

# Initial conditions setup
class InitialConditions(Expression):
    def __init__(self):
        print "initialized initial conditions"
    def eval(self, values, x):
        values = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
    def value_shape(self):
        return (1,)

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

def update(self):
    plot(self.pdesystems['Scalar'].u_)
Problem.update = update

problem.prm['T'] = 0.3
problem.solve()

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:
Revision history for this message
Yaakoub El Khamra (yelkhamra) said : 		#1

So turns out I am missing:
problem.initialize(poisson)

before the problem.solve() but now I just get:
Error in initialize! Provide tuples of strings, Constants or Expressions.

Not quite sure how to fix that, or why it happened to be honest.

Revision history for this message
Yaakoub El Khamra (yelkhamra) said : 		#2

Turns out if I do this:
Init = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
problem.q0 = {'u':Init}

then everything works just fine.

Revision history for this message
Mikael Mortensen (mikael-mortensen) said : 		#3

You can even do

problem.q0 = {'u': "10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02"}

Mikael

Den Sep 13, 2012 kl. 1:35 AM skrev Yaakoub El Khamra:

> Question #208399 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/208399
>
> Yaakoub El Khamra gave more information on the question:
> Turns out if I do this:
> Init = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
> problem.q0 = {'u':Init}
>
> then everything works just fine.
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Revision history for this message
Mikael Mortensen (mikael-mortensen) said : 		#4

Hi,

Your code actually made me realize I had made an error in my recent fix of solve_Newton_system. I fixed it back and now the code should work for your problem using either Picard or Newton. For your linear problem the two approaches should actually be identical.

Best regards

Mikael

Den Sep 13, 2012 kl. 1:35 AM skrev Yaakoub El Khamra:

> Question #208399 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/208399
>
> Yaakoub El Khamra gave more information on the question:
> Turns out if I do this:
> Init = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
> problem.q0 = {'u':Init}
>
> then everything works just fine.
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Revision history for this message
Yaakoub El Khamra (yelkhamra) said : 		#5

Thank you very much for your help sir!
I applied the update and everything works fine. I still do not understand why the InitialConditions() fails in my example but works just fine in the CahnHilliard example. Can you please shed some light on that?

Also the output now looks as follows:

Creating new work vector for u
Adding PDESubSystem: Poisson
Adding ['linear_solver']['u'] = lu to pdesubsystem u
Adding ['iteration_type'] = Newton to pdesubsystem u
Adding ['omega']['u'] = 1.0 to pdesubsystem u
    Iter 1 error | u | 2.2797e+03 2.4678e+03 |
Time = 0.01, End time = 0.1
    Iter 2 error | u | 1.1192e+03 6.8668e+02 |
Time = 0.02, End time = 0.1
    Iter 3 error | u | 6.2009e+02 4.7639e+02 |
Time = 0.03, End time = 0.1
    Iter 4 error | u | 4.5160e+02 3.7356e+02 |
Time = 0.04, End time = 0.1
    Iter 5 error | u | 3.5907e+02 3.0462e+02 |
Time = 0.05, End time = 0.1
    Iter 6 error | u | 2.9413e+02 2.5169e+02 |
Time = 0.06, End time = 0.1
    Iter 7 error | u | 2.4340e+02 2.0887e+02 |
Time = 0.07, End time = 0.1
    Iter 8 error | u | 2.0210e+02 1.7360e+02 |
Time = 0.08, End time = 0.1
    Iter 9 error | u | 1.6800e+02 1.4435e+02 |
Time = 0.09, End time = 0.1
    Iter 10 error | u | 1.3971e+02 1.2005e+02 |
Time = 0.1, End time = 0.1

Pardon my ignorance but if those are errors, they are pretty high are they not? Or am I missing something?

Regards
Yaakoub

Revision history for this message
Best Mikael Mortensen (mikael-mortensen) said : 		#6

Hi Yaakoub

Your equation is linear so the error after one iteration is actually zero. The error that is plotted above is the norm of the rhs and the norm of the vector of change. To get a fair estimate of the error you need to recompute b after the iteration step. Problem is this is expensive and I'm not sure I can justify putting in such a check for all cases in general. The reason it did not work for me yesterday was that I had tried to incorporate such a check. It worked well for steady problems, but not for transient like yours. I will look into it and see if I can find a way around this, meanwhile you can use just one iteration on the time step or put this

poisson.prm['max_iter'] = 2

before solving. You will see that the error after two time steps is zero.

Regarding:

# Initial conditions setup
class InitialConditions(Expression):
    def __init__(self):
        print "initialized initial conditions"
    def eval(self, values, x):
        values = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
    def value_shape(self):
        return (1,)

The eval method is called for each and every element in the mesh and you cannot put an Expression on the rhs of values. You should have done something like this:

# Initial conditions setup
class InitialConditions(Expression):
    def __init__(self):
        print "initialized initial conditions"
    def eval(self, values, x):
        values[0] = 10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)

but I guess that to make this work you need to use the numpy exp, not the ufl exp that you get when you do from dolfin import *.

Best regards

Mikael

Revision history for this message
Yaakoub El Khamra (yelkhamra) said : 		#7

Mikael thank you very, very much!

Regards
Yaakoub

Revision history for this message
Yaakoub El Khamra (yelkhamra) said : 		#8

Thanks Mikael Mortensen, that solved my question.