# Debugging recommendations

Asked by Yaakoub El Khamra on 2012-09-13

Mikael I am very sorry to keep bothering you with incessant questions but I need some help debugging a simplified coupled system. I am getting an [0]PETSC ERROR: Floating point exception! error from PETSc when I do the solve and quite frankly I think I have my linear forms wrong. Can you please have a look and let me know if you spot anything? Also do you have any recommendations as to how I would go about debugging errors of this nature? Thank you very much in advance.

from cbc.pdesys import *

mesh = UnitSquare(10, 10)
# 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['familyname'] = 'Scalar'
solver_parameters['iteration_type'] = 'Newton'
coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

class SS2Ph(PDESubSystem):
def form(self, Pw, v_Pw, Sw, Sw_1, v_Sw, dt, **kwargs): # v_u is the TestFunction
return FPw + FSw

# setup the BC
bcs = [DirichletBC(coupled.V['PwSw'], (1.0, 1.0), "on_boundary")]

class InitialConditions(Expression):
def __init__(self):
print "initialized initial conditions"
def eval(self, values, x):
values[0] = 1.0
values[1] = 1.0
def value_shape(self):
return (2,)

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

problem.prm['T'] = 1.5
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:
 Mikael Mortensen (mikael-mortensen) said on 2012-09-13: #1

Den Sep 13, 2012 kl. 9:31 PM skrev Yaakoub El Khamra:

> New question #208509 on CBC.PDESys:
>
> Mikael I am very sorry to keep bothering you with incessant questions but I need some help debugging a simplified coupled system. I am getting an [0]PETSC ERROR: Floating point exception! error from PETSc when I do the solve and quite frankly I think I have my linear forms wrong. Can you please have a look and let me know if you spot anything? Also do you have any recommendations as to how I would go about debugging errors of this nature? Thank you very much in advance.

No worries, just happy someone can find this useful:-)

>
> from cbc.pdesys import *
>
> mesh = UnitSquare(10, 10)
> # 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['familyname'] = 'Scalar'
> solver_parameters['iteration_type'] = 'Newton'
> coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.
>
> class SS2Ph(PDESubSystem):
> def form(self, Pw, v_Pw, Sw, Sw_1, v_Sw, dt, **kwargs): # v_u is the TestFunction
> return FPw + FSw
>

I don't understand this? What does your system of equations look like? Shouldn't you at least have 1.0/dt*(Pw - Pw_1)*v_Pw in the first line here? Other than the weird equations I can't really see anything 'wrong' with the code.

Best regards

Mikael

> # setup the BC
> bcs = [DirichletBC(coupled.V['PwSw'], (1.0, 1.0), "on_boundary")]
> # add the pde subsystem
>
> class InitialConditions(Expression):
> def __init__(self):
> print "initialized initial conditions"
> def eval(self, values, x):
> values[0] = 1.0
> values[1] = 1.0
> def value_shape(self):
> return (2,)
>
> problem.q0 = {'PwSw': InitialConditions()}
>
> problem.prm['T'] = 1.5
> problem.initialize(coupled)
> problem.solve()
>
> --
> contact for CBC.PDESys.

 Yaakoub El Khamra (yelkhamra) said on 2012-09-13: #2

Thank you, FEniCS and the CBC PDESys are more than helpful. I am re-writing my dissertation to use CBC PDESys exclusively. Essentially I am solving a simplified version of the Darcy equations for two phase flow. I was also about to write you regarding the error, turns out PEBKAC. The equations I am trying to solve are simplified as:
Laplacian(Pw) = d(Sw)/dt
Laplacian(Pw) = d(So)/dt
Where So = 1- Sw. And that's where my mistake was, I forgot to reverse the sign for d(Sw)/dt for the second equation. That takes care of the Floating point exception error.

Once I am done I will distribute a flow through porous media toolkit, along the lines of RANS. Thank you very much for putting up with me Mikael!

Regards
Yaakoub

 Mikael Mortensen (mikael-mortensen) said on 2012-09-13: #3

PEBKAC, lol, had to look it up in an urban dictionary so I must be a lesser nerd than I thought I was.

Mikael

Den Sep 13, 2012 kl. 10:11 PM skrev Yaakoub El Khamra:

> Question #208509 on CBC.PDESys changed:
>
>
> Yaakoub El Khamra confirmed that the question is solved:
>
> Thank you, FEniCS and the CBC PDESys are more than helpful. I am re-writing my dissertation to use CBC PDESys exclusively. Essentially I am solving a simplified version of the Darcy equations for two phase flow. I was also about to write you regarding the error, turns out PEBKAC. The equations I am trying to solve are simplified as:
> Laplacian(Pw) = d(Sw)/dt
> Laplacian(Pw) = d(So)/dt
> Where So = 1- Sw. And that's where my mistake was, I forgot to reverse the sign for d(Sw)/dt for the second equation. That takes care of the Floating point exception error.
>
> Once I am done I will distribute a flow through porous media toolkit,
> along the lines of RANS. Thank you very much for putting up with me
> Mikael!
>
> Regards
> Yaakoub
>
> --