[heliana60]

Hello Bruno!!

I hope you are fine.!

I am developing a new fluidity model and using viscoelastic solver to compile it. I dont know if you are familiar with this solver, but I just inserted an equation for fluidity (which is an inverted relaxation time) which the stress depends on.

I can compile the model but when I try to run a simple shear example I get this error almost inmediately:

Courant Number mean: 3.24982e-10 max: 2.12946e-08
Time = 4e-06

DILUPBiCG:  Solving for Ux, Initial residual = 0.0671588, Final residual = 5.32629e-12, No Iterations 5
DILUPBiCG:  Solving for Uy, Initial residual = 0.170076, Final residual = 2.34414e-11, No Iterations 5
DICPCG:  Solving for p, Initial residual = 0.000905894, Final residual = 8.4895e-11, No Iterations 66
time step continuity errors : sum local = 3.7e-26, global = 1.4064e-27, cumulative = 4.73523e-27
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#6  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#7  Foam::SolverPerformance<double> Foam::solve<double>(Foam::tmp<Foam::fvMatrix<double> > const&) at ??:?
#8  Foam::fluiditySSDiffusion::correct() at ??:?
#9  
 at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  
 at ??:?
Floating point exception

So, as far as I was reading there must be a division by zero or something bad in the mathematics.

my equation is a steady state equation, so i am not using time derivaties. it goes like this:

// fluidity model
solve
(
    - a_*(Foam::sqr(flu_)-2*flu_/lambda_+1/sqr(lambda_))+r_ * 0.5 * sqr(mag(twoD)) + fvm::laplacian(D_,flu_, "laplacian(D,flu)")

)

;
//    fluEqn.relax();
//    fluEqn.solve();

////
// Stress transport equation
    fvSymmTensorMatrix tauEqn
    (
        fvm::ddt(tau_)
      + fvm::div(phi(), tau_)
     ==
        etaP_/lambda_*twoD
      + twoSymm(C)
    - fvm::Sp(flu_, tau_)
    );


    tauEqn.relax();
    tauEqn.solve();

I have figured the square is the problem, because when i removed it it worked fine.. but i need it tho

Do you have any idea how I can make this work. If you need more details about the case please let me know ,


Thank you and have a nice day !

Heliana

[wyldckat]

Greetings Heliana,

I've moved your post from the other thread: http://www.cfd-online.com/Forums/ope...ntstack-3.html - because this specifically relates to the development of a new solver/model.

First thing you should do: isolate and conquer. There is a wiki page on the topic of debugging: http://openfoamwiki.net/index.php/HowTo_debugging
And recently I wrote a post with an example on how to do old-school debugging, namely to output as much information as possible, to see what's truly happening: http://www.cfd-online.com/Forums/ope...tml#post512967 post #9

I don't have much time to ask for more details, so the best i can do is to tell you for you to look for the details yourself
More specifically, try to see what exact values do the more relevant fields have before the crash. For example:

Info << "flu_: " << flu_ << endl;

Best regards,
Bruno