[omican]

Hello all!

I am trying to solve the following system in 2D:

dc1/dt = b*c2*c2 + L1 laplace c1 + n L4 laplace phi
dc2/dt = 2a - 2b*c2*c2 + L3 laplace c1 + L4 laplace phi
laplace phi = - L2 * c2

for which I made a straightforward solver like this:

for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readSIMPLEControls.H"
scalar maxResidual = 0;

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix EqC1
(
fvm::ddt(C1)
==
b*C2*C2
+
fvm::laplacian(L1, C1)
+
fvc::laplacian(n*L4, elstat_pot)
);
EqC1.relax();

scalar eqnResidualC1 = EqC1.solve().initialResidual();
maxResidual = max(eqnResidualC1, maxResidual);

fvScalarMatrix EqC2
(
fvm::ddt(C2)
==
2*a-2*b*C2*C2
+
fvc::laplacian(L3, C1)
+
fvc::laplacian(L4, elstat_pot)
);
EqC2.relax();

scalar eqnResidualC2 = EqC2.solve().initialResidual();
maxResidual = max(eqnResidualC2, maxResidual);

fvScalarMatrix EqPot
(
fvm::laplacian(elstat_pot)
==
-L2*C2
);
EqPot.relax();

scalar eqnResidualPot = EqPot.solve().initialResidual();
maxResidual = max(eqnResidualPot, maxResidual);

}

# include "write.H"

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

It works fine with some cases, but crashes for others. For example when I use the following values of parameters:

a=200
b=6.324555e-6
n=1.5
L1=6.0e-10
L2=1.089733454e+9
L3=1.19762e-10
L4=3.10929e-5

the solver works pretty fine. But when I raise L2 to 1.08973e+15 (which is physically more reasonable), only a few time levels are computed and then the solver crashes saying

#0 Foam::error::printStack(Foam::Ostream&) in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<d ouble>&, Foam::Field<double> const&) const in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#5 Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMa trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#6 Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/sporit/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libfiniteVolume.so"
#8 main in "/home/sporit/OpenFOAM/sporit-1.6/applications/bin/linuxGccDPOpt/cellFoamForum"
#9 __libc_start_main in "/lib/libc.so.6"
#10 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/home/sporit/OpenFOAM/sporit-1.6/applications/bin/linuxGccDPOpt/cellFoamForum"
Floating point exception(SIGFPE)

I've tried various numerical schemes as well as various linear-solvers, but the only settings that really helps is to reduce deltaT. However, for some values (like the ones shown above) deltaT has to be set to something like 5e-11 which makes the process extremely slow. Any advice or recommendation how to speed up the computation in this case will be highly appreciated.

Thanks a lot in advance,
Ondra