On Thu, Oct 17, 2013 at 12:00 PM, Christophe Ortiz < [email protected]> wrote:
> > > On Wed, Oct 16, 2013 at 4:22 PM, Christophe Ortiz < > [email protected]> wrote: > >> >> >> On Wed, Oct 16, 2013 at 4:10 PM, Christophe Ortiz < >> [email protected]> wrote: >> >>> >>> >>> On Wed, Oct 16, 2013 at 3:36 PM, Jed Brown <[email protected]> wrote: >>> >>>> Christophe Ortiz <[email protected]> writes: >>>> > I found where is the problem. I had >>>> > set TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE); >>>> > I found it produces wrong solutions, either with ARKIMEX3, 1BEE or A2. >>>> > Solution becomes negative. Without this option, all ARKIMEX types >>>> work very >>>> > well and very fast. >>>> >>>> I'm glad they work as IMEX methods (please check that the solutions are >>>> correct), but I would like to find out why they are not working in fully >>>> implicit mode. Is your problem similar to ex25? >>> >>> >>> Yes. The only difference is that I put everything (diffusion and >>> reaction terms) under IFunction and IJacobian. >>> >>> My system is >>> >>> u_t - alpha.u_xx + k.u.v = 0 >>> v_t - beta.v_xx + k.u.v = 0 >>> >>> with Dirichlet boundary conditions >>> u(x=0,t) = u(x=L,t) = cste >>> v(x=0,t) = v(x=L,t) = cste >>> >>> (If I remember >>>> correctly, that example uses a DAE formulation for boundary conditions >>>> and thus has trouble with some methods with explicit stages. There may >>>> be more going on.) >>>> >>> >>> Now you say it, I'm having some troubles. After checking that it worked >>> very well with only diffusion terms with ARKIMEX, I introduced the reaction >>> term. With ARKIMEX 1BEE/A2, it goes very fast up to t=0.27 s, then, steps >>> start being rejected and then it goes very slowly, with small timesteps. >>> See output below. >>> Do you think it is due to the fact that I use IFunction to define >>> boundary conditions, which produces troubles with ARKIMEX methods ? >>> What would be the solution ? >>> >>> Christophe >>> >>> >> I just tried with TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE); It goes much >> better and reaches final time....but solution is clearly not correct. The >> initial gaussian distributions do not change. Like if diffusion was null. >> > > Hi, > > I've tried something else. Since ARKIMEX is able to converge up to t~0.20s > and gives the correct solution for this time, and then comes into trouble, > I thought that it could come from the SNES. > Then, instead of using SNESNEWTONLS, I tried SNESNEWTONTR. I summarize my > findings: > > - With NEWTONTR, ARKIMEX (3, 1BEE or A2) converges. > > -With NEWTONTR, ARKIMEX 3 converges much faster than 1BEE or A2 (much less > timesteps). > > - No conflicts between NEWTONTR and ARKIMEXFullyImplicit as with NEWTONLS. > Give the same results. > > > Then I came back to SNESNEWTONLS > > - With NEWTONLS, ARKIMEX does not converges. > > - When ARKIMEXFullyImplicit is used with NEWTONLS, it converges, but > solution clearly incorrect. It seems it does nothing, initial distributions > do not change. > > > And I played a bit with the tolerances: > > - No influence of rtol (increasing or decreasing it). It stops at same > point. > > - Changing stol from 1e-8 to 1e-7 seems to ease as it goes up to longer > time (6.98 s instead of 0.20) but then stops with an error: > > [0]PETSC ERROR: --------------------- Error Message > ------------------------------------ > [0]PETSC ERROR: Floating point exception! > [0]PETSC ERROR: Vec entry at local location 450 is not-a-number or > infinite at end of function: Parameter number 3! > ... > ... > [0]PETSC ERROR: VecValidValues() line 30 in src/vec/vec/interface/rvector.c > [0]PETSC ERROR: SNESComputeFunction() line 1998 in > src/snes/interface/snes.c > [0]PETSC ERROR: SNESSolve_NEWTONLS() line 162 in src/snes/impls/ls/ls.c > [0]PETSC ERROR: SNESSolve() line 3636 in src/snes/interface/snes.c > [0]PETSC ERROR: TSStep_ARKIMEX() line 765 in src/ts/impls/arkimex/arkimex.c > [0]PETSC ERROR: TSStep() line 2458 in src/ts/interface/ts.c > [0]PETSC ERROR: TSSolve() line 2583 in src/ts/interface/ts.c > [0]PETSC ERROR: main() line 392 in src/ts/examples/tutorials/diffusion.c > > Maybe all this will give you a hint of what's going on... > Saludos, > Christophe > > > I've done additional tests: - NEWTONLS does not work with the default LineSearch, ie bt. - NEWTONLS works very well with LineSearch basic, l2 and cp. Converges fast and results are correct. Seems to be due to the LineSearch bt. Christophe
