# Debugging recommendations

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_

problem = Problem(mesh, problem_parameters)

solver_

solver_

coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

class SS2Ph(PDESubSys

def form(self, Pw, v_Pw, Sw, Sw_1, v_Sw, dt, **kwargs): # v_u is the TestFunction

FPw = inner(grad(Pw), grad(v_Pw))*dx - 1.0/dt*(Sw - Sw_1)*v_Pw*dx

FSw = inner(grad(Pw), grad(v_Sw))*dx - 1.0/dt*(Sw - Sw_1)*v_Sw*dx

return FPw + FSw

# setup the BC

bcs = [DirichletBC(

# add the pde subsystem

coupled.

class InitialConditio

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': InitialConditio

problem.prm['T'] = 1.5

problem.

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:

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

> New question #208509 on CBC.PDESys:

> https:/

>

> 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_

> problem = Problem(mesh, problem_parameters)

>

> solver_

> solver_

> coupled = PDESystem([['Pw', 'Sw']], problem, solver_parameters) # Creates FunctionSpace, Functions etc.

>

> class SS2Ph(PDESubSys

> def form(self, Pw, v_Pw, Sw, Sw_1, v_Sw, dt, **kwargs): # v_u is the TestFunction

> FPw = inner(grad(Pw), grad(v_Pw))*dx - 1.0/dt*(Sw - Sw_1)*v_Pw*dx

> FSw = inner(grad(Pw), grad(v_Sw))*dx - 1.0/dt*(Sw - Sw_1)*v_Sw*dx

> 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(

> # add the pde subsystem

> coupled.

>

> class InitialConditio

> 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': InitialConditio

>

> problem.prm['T'] = 1.5

> problem.

> problem.solve()

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.

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

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:

> https:/

>

> Status: Answered => Solved

>

> 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

>

> --

> You received this question notification because you are an answer

> contact for CBC.PDESys.