# Initial conditions

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_

problem = Problem(mesh, problem_parameters)

# setup the solver parameters

solver_

solver_

# setup the pdesystem

poisson = PDESystem([['u']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

class Poisson(

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(

# add the pdesubsystem

poisson.

# Initial conditions setup

class InitialConditio

def __init__(self):

print "initialized initial conditions"

def eval(self, values, x):

values = Expression(

def value_shape(self):

return (1,)

problem.q0 = {'u':InitialCon

def update(self):

plot(

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:

Yaakoub El Khamra (yelkhamra) said : | #1 |

So turns out I am missing:

problem.

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.

Yaakoub El Khamra (yelkhamra) said : | #2 |

Turns out if I do this:

Init = Expression(

problem.q0 = {'u':Init}

then everything works just fine.

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

>

> Yaakoub El Khamra gave more information on the question:

> Turns out if I do this:

> Init = Expression(

> problem.q0 = {'u':Init}

>

> then everything works just fine.

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.

Hi,

Your code actually made me realize I had made an error in my recent fix of solve_Newton_

Best regards

Mikael

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

> Question #208399 on CBC.PDESys changed:

> https:/

>

> Yaakoub El Khamra gave more information on the question:

> Turns out if I do this:

> Init = Expression(

> problem.q0 = {'u':Init}

>

> then everything works just fine.

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.

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_

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

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.

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

Regarding:

# Initial conditions setup

class InitialConditio

def __init__(self):

print "initialized initial conditions"

def eval(self, values, x):

values = Expression(

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 InitialConditio

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

Yaakoub El Khamra (yelkhamra) said : | #7 |

Mikael thank you very, very much!

Regards

Yaakoub

Yaakoub El Khamra (yelkhamra) said : | #8 |

Thanks Mikael Mortensen, that solved my question.