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
        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(coupled.V['PwSw'], (1.0, 1.0), "on_boundary")]
# 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] = 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:
Last reply:

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

> New question #208509 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/208509
>
> 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
> 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(coupled.V['PwSw'], (1.0, 1.0), "on_boundary")]
> # 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] = 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()
>
> --
> 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://answers.launchpad.net/cbcpdesys/+question/208509
>
> 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.