On Tue, 14 Mar 2006, Nils Erik Svangård wrote:
Hi, Yes, there is a difference between the old 3-stage RK-solver and the SNESsolver using RK. In the old a timestep was calculatade for each
^^^^^^^^^^^^^^^
Are you using PETSc's explicit RK? If so that does not use SNES.
cell based on a stability criteria and if chosen it could run in time-accurate mode eg. the timestep is the same for each cell, in SNES
^^^^^^^^^ Do you mean for backward Euler for the explicit RK? Or both?
I use a timestep that is the same for all cells and changes only if it returns non-physical answers or the solution diverges.
I think you need to put the "proper" computation of the time-step based on the stability criteria that you have in the "old" code. Unfortunately this is currently hardwired in src/ts/impls/explicit/rk/rk.c (search for "Computing next stepsize"). I think you need to use the minimum of the timestep computed here and what is needed for the "stability criteria".
Another interesting observation is that when using back euler at the moment when it bugs out petsc has assigned neigboring cells pressure values that differs about 1 ATM, which is a very non physical value,
Likely some stability problem. Barry
but I have no clue as to why petsc want such a value in thesolution. /nisse On 3/3/06, Barry Smith <bsmith@mcs.anl.gov> wrote:There must be a difference between the two RKs. The initial timestep it uses, or algorithm for changing time step or slight difference in algorithms? Barry On Fri, 3 Mar 2006, Nils Erik Svangård wrote:A interesting thing is when I change TS_BEULER to TS_RUNGE_KUTTA the solver exits fast due to the non physical values. My original solver is a 3-step Runge-Kutta which works flawlessly. The error when running RK: ./ImplicitG3D on a linux-gnu named sethnx004.vac.com by yy26539 Fri Mar 3 10:31:27 2006 Libraries linked from /home/yy26539/work/NISSE/petsc-2.3.0/lib/linux-gnu Configure run at Thu Oct 13 08:23:56 2005 Configure options --with-cc=gcc --with-fc="f77 -N109" --download-mpich=1 --download-mpich-pm=gforker --download-f-blas-lapack=1 --with-shared=0 ----------------------------------------------------------------------- [0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [0]PETSC ERROR: likely location of problem given in stack below [0]PETSC ERROR: --------------- Stack Frames --------------- [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available, [0]PETSC ERROR: INSTEAD the line number of the start of the function [0]PETSC ERROR: is given. [0]PETSC ERROR: [0] TS user right-hand-side function line 0 unknownunknown [0]PETSC ERROR: [0] TSComputeRHSFunction line 282 src/ts/interface/ts.c [0]PETSC ERROR: [0] TSRkqs line 256 src/ts/impls/explicit/rk/rk.c [0]PETSC ERROR: [0] TSStep_Rk line 322 src/ts/impls/explicit/rk/rk.c [0]PETSC ERROR: [0] TSStep line 1520 src/ts/interface/ts.c [0]PETSC ERROR: -------------------------------------------- [0]PETSC ERROR: User provided function() line 0 in unknown directory unknown file [0]PETSC ERROR: Signal received! [0]PETSC ERROR: ! aborting job: application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0 On 3/3/06, Nils Erik Svangård <nilserik@gmail.com> wrote:Hello again,Is it just not compiling because PETSC_ERR_ARG_DOMAIN is not defined? Or does it still not work even when you have put the number in? What is the error? Or did you get it all working with SNES?It is not compiling because PETSC_ERR_ARG_DOMAIN is not defined: error on line 585 of SNESG3D2.F: data type is undefined for variable PETSC_ERR_ARG_DOMAIN I tried putting a number in, then it compiles but complains while linking that the Push() and Pop() functions is not defined. I made a workaround by adding a variable in a common block and setting it before returning from a non physical value in the FormFunction. Now I can compile and when there is a non-physical value it returns and resumes with a smaller timestep. SNES still give non physical values which are hard to add checks for in the program (e.g sqrt(RO(L)-RO(L+1))). Did I understand correctly that SNES does a linesearch and probably goes for a solution which isnt the one I looking for? Does SNES use the linesearch algorithm when I use -snes_mf (which I use)?It may work automatically with TS; you would just set TS with a smaller timestep and call TSSolve again?OK, I'll add a check and put TSSolve in a while loop. My SNESSolve loop looks like this right now: call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr) call VecDuplicate(x,r,ierr) call VecDuplicate(x,xb,ierr) call SNESCreate(PETSC_COMM_SELF,snes,ierr) call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr) call SNESSetFromOptions(snes,ierr) call FormInitialGuess(x,ierr) call VecCopy(x,xb,ierr) write(6,*)zero,dt,tmax,i1000 step=0 timestep=dt bang=0 C While loop 10 if (step .le. 100) then C call PetscExceptionPush(87,ierr) C call PetscExceptionPush(PETSC_ERR_ARG_DOMAIN,ierr) call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr) ierr=bang C call PetscExceptionPop(87,ierr2) IF ( ierr .EQ. 87) THEN timestep=timestep/2 call VecCopy(xb,x,ierr) write(6,*) "No physical solution, trying", . " with smaller timestep: ",timestep bang=0 ELSE call SNESGetConvergedReason(snes,reason,ierr) IF ( reason .LT. 0 ) THEN timestep=timestep/2 call VecCopy(xb,x,ierr) write(6,*)"Solution diverged, trying with", . " smaller timestep: ",timestep ELSE timestep=dt step=step+1 call VecCopy(x,xb,ierr) write(6,*) "Solve success, going for step: ",step write(6,*) "Reason: ",reason," ierr: ",ierr reason=0 ierr=0 ENDIF ENDIF goto 10 endif C End While loop And the check in the FormFunction: if (ET(L) .LT. 0) then ierr = 87 bang = 87 C SETERRQ(87,"Non physical value!",ierr) return elseif ( RO(L) .LT. 0 ) then ierr = 87 bang = 87 return endif /nisseBarry I'm adding a petscerror.h for Fortran that includes these values Sorry I forgot we hadn't given them for Fortran. On Thu, 2 Mar 2006, Nils Erik Svangård wrote:Hi again, I realised that my implementation was a bit flawed so I have rewritten it again, the push and pop() functions dont seem to work, and they are keeping me from compiling the program, this is the main timestep loop: call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr) call VecDuplicate(x,r,ierr) call VecDuplicate(x,xb,ierr) call SNESCreate(PETSC_COMM_SELF,snes,ierr) call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr) call SNESSetFromOptions(snes,ierr) call FormInitialGuess(x,ierr) call VecCopy(x,xb,ierr) write(6,*)zero,dt,tmax,i1000 step=0 timestep=dt C While loop 10 if (step .le. 100) then call PetscExceptionPush(87,ierr) call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr) call PetscExceptionPop(87,ierr2) call SNESGetConvergedReason(snes,reason,ierr) IF ( reason .LT. 0 ) THEN timestep=timestep/2 VecCopy(xb,x,ierr) write(6,*)"Solution diverged, trying with", . " smaller timestep: ",timestep ELSE IF ( ierr .EQ. 87) THEN timestep=timestep/2 VecCopy(xb,x,ierr) write(6,*) "No physical solution, trying", . " with smaller timestep: ",timestep ELSE timestep=dt step=step+1 VecCopy(x,xb,ierr) write(6,*) "Solve success, going for step: ",step ENDIF ENDIF goto 10 endif C End While loop And in my Formfunction I have added the following check: if (ET(L) .LT. 0) then ierr = 87 return endif PETSC_ERR_ARG_DOMAIN didnt work: error on line 6500 of SNESG3D2.F: data type is undefined for variable PETSC_ERR_ARG_DOMAIN Even thought I had petsc.h included. Can I do this check when I use TS also? I have the sameproblem with the values being negative when using TS. /nisse On 2/27/06, Barry Smith <bsmith@mcs.anl.gov> wrote:SNES works by computing p = -approxinv(J)*F(uold) and then does a line search on unew = uold + lambda*p to get the new u. First it uses a test value of 1 for lambda so it tries to compute F(uold + p). It is possible that uold + p has some "non-physical" values in it. There are two ways you can try handling it: 1) Before doing the linesearch SNES calls a "precheck" function, that can change the step if it decides there is a problem with the step (like it is too long). You can provide your own precheck function with SNESLineSearchSetPreCheck() it could, for example, shrink the Newton direction to make it remain physical. or 2) when your form function detects an illegal value it calls SETERRQ(PETSC_ERR_ARG_DOMAIN,"Nonphysical function input"); Note: this is before it takes the square root of the negative number, so check the number and call the error before calling SETERRQ(). (If using Fortran then simply set ierr = PETSC_ERR_ARG_DOMAIN and then immediately return). Then replace your call to SNESSolve with ierr = PetscExceptionTry1(SNESSolve(snes,b,x),PETSC_ERR_ARG_DOMAIN); if (PetscExceptionCaught(ierr,PETSC_ERR_ARG_DOMAIN)) { /* this means your function found a non-physical value so cut your time-step and continue through the loop again. */ Put code to do this here. } If using fortran then do call PetscExceptionPush(PETSC_ERR_ARG_DOMAIN,ierr) call SNESSolve(snes,b,x,ierr) call PetscExceptionPop(PETSC_ERR_ARG_DOMAIN,anotherierr) if (ierr == PETSC_ERR_ARG_DOMAIN) then non-physical value so cut the time-step and try again else ! everything is normal so take the next time-step Barry The Fortran interface may be missing PetscExceptionPush() and Pop() if so let us know and we'll provide the patch. On Mon, 27 Feb 2006, Nils Erik Svangård wrote:Hi all, I have problems solving some CFD problems using SNES and my custom back-euler. I have 7 equations that I want to solve. First my FormFunction copies the values from the Vec that SNES uses to the variables that my code use: RO(L)=xx(1,L) RU(L)=xx(2,L) RV(L)=xx(3,L) RW(L)=xx(4,L) ET(L)=xx(5,L) RQ(L)= xx(6,L) REPS(L)=xx(7,L) Here I also print the values of ET(2) for debugging: write(6,*)"ET(2) = xx(5,2) : ",ET(2)," = ",xx(5,2) Then I get the fluxes by calling custom functions call AUXVR call VGRAD call FLUX call KESRC Then I perform back-euler save the new value to so that I can use it in the next iterationi (TSF(L) is Time step function, which is set by hand and is the same for all L): ff(1,L) = RO(L)-OLD(1,L)-TSF(L)*DRO(L) ff(2,L) = RU(L)-OLD(2,L)-TSF(L)*DRU(L) ff(3,L) = RV(L)-OLD(3,L)-TSF(L)*DRV(L) ff(4,L) = RW(L)-OLD(4,L)-TSF(L)*DRW(L) ff(5,L) = ET(L)-OLD(5,L)-TSF(L)*DET(L) ff(6,L) = RQ(L)-OLD(6,L)-TSF(L)*DRQ(L) ff(7,L) = REPS(L)-OLD(7,L)-TSF(L)*DREPS(L) And save the new value of RO-REPS to use in the next iteration of back-euler: OLD(1,L)=RO(L) OLD(2,L)=RU(L) OLD(3,L)=RV(L) OLD(4,L)=RW(L) OLD(5,L)=ET(L) OLD(6,L)=RQ(L) OLD(7,L)=REPS(L) Here I print the values of L=2 as above for debugging: write(6,*)ff(5,2)," = ",ET(2),"-",OLD(5,2),"-",TSF(2),"*",DET(2) The program exit abnormaly after 33 runs of the FormFunction. The cause of this is that AUXVR tries to perform sqrt(ET(2)) when ET(2) is negative, ET is the total energy and should never be negative. It seems that all of a sudden the PETSc SNES solver decides to supply the FormFunction with a negative ET value. I need help understanding why and how to fix it (if it is fixable). Here is the output of: write(6,*)"ET(2) = xx(5,2) : ",ET(2)," = ",xx(5,2) Which is in the beginning of my FormFunction. And: write(6,*)ff(5,2)," = ",ET(2),"-",OLD(5,2),"-",TSF(2),"*",DET(2) Which is in the end of my Formfunction. 1. ET(2) = xx(5,2) : 253250. = 253250.000000000 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 2. ET(2) = xx(5,2) : 253250. = 253250.000010259 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 3. ET(2) = xx(5,2) : 253250. = 253249.999999988 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4132419437260 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941324E+08 4. ET(2) = xx(5,2) : 253250. = 253250.000048432 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 5. ET(2) = xx(5,2) : 253250. = 253249.999999872 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4177283437784 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941773E+08 6. ET(2) = xx(5,2) : 253250. = 253250.000071113 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 7. ET(2) = xx(5,2) : 253250. = 253249.999980115 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4453347441010 = 253250. - 253250.000000000 - 1.000000E-07 * -2.944533E+08 8. ET(2) = xx(5,2) : 253250. = 253250.000037298 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4082723436679 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940827E+08 9. ET(2) = xx(5,2) : 253250. = 253249.999942706 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4082723436679 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940827E+08 10. ET(2) = xx(5,2) : 253250. = 253249.999948703 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 11. ET(2) = xx(5,2) : 253250. = 253250.000007808 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.7283811474088 = 253250. - 253250.000000000 - 1.000000E-07 * -2.972838E+08 12. ET(2) = xx(5,2) : 253250. = 253250.000018146 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 13. ET(2) = xx(5,2) : 253250. = 253250.000001865 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 28.3742819315846 = 253250. - 253250.000000000 - 1.000000E-07 * -2.837428E+08 14. ET(2) = xx(5,2) : 253250. = 253249.999997296 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 15. ET(2) = xx(5,2) : 253250. = 253250.000000049 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 30.6548963582361 = 253250. - 253250.000000000 - 1.000000E-07 * -3.065490E+08 16. ET(2) = xx(5,2) : 253250. = 253250.000022870 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 17. ET(2) = xx(5,2) : 253250. = 253249.999999716 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.5326051451209 = 253250. - 253250.000000000 - 1.000000E-07 * -2.953260E+08 18. ET(2) = xx(5,2) : 253250. = 253250.000074328 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4165187437643 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941652E+08 19. ET(2) = xx(5,2) : 253250. = 253250.000049809 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4085059436707 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940851E+08 20. ET(2) = xx(5,2) : 253250. = 253250.000019253 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 21. ET(2) = xx(5,2) : 253250. = 253250.000000189 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.8851075492403 = 253250. - 253250.000000000 - 1.000000E-07 * -2.988511E+08 22. ET(2) = xx(5,2) : 253250. = 253249.999997897 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4085059436707 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940851E+08 23. ET(2) = xx(5,2) : 253250. = 253249.999975142 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.3961539435263 = 253250. - 253250.000000000 - 1.000000E-07 * -2.939615E+08 24. ET(2) = xx(5,2) : 253250. = 253250.000018324 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4085059436707 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940851E+08 25. ET(2) = xx(5,2) : 253250. = 253249.999987073 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4165187437643 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941652E+08 26. ET(2) = xx(5,2) : 253250. = 253249.999980982 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4165187437643 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941652E+08 27. ET(2) = xx(5,2) : 253250. = 253249.999994645 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 28. ET(2) = xx(5,2) : 253250. = 253249.999962364 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4165187437643 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941652E+08 29. ET(2) = xx(5,2) : 253250. = 253249.999999921 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 30. ET(2) = xx(5,2) : 253250. = 253250.000000058 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 28.4019971319085 = 253250. - 253250.000000000 - 1.000000E-07 * -2.840200E+08 31. ET(2) = xx(5,2) : 253250. = 253250.000040348 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4085059436707 = 253250. - 253250.000000000 - 1.000000E-07 * -2.940851E+08 32. ET(2) = xx(5,2) : 253250. = 253250.000005245 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 33. ET(2) = xx(5,2) : 253250. = 253250.000005245 ff(5,2) = ET(2)-OLD(5,2)-TSF(2)*DET(2) : 29.4126563437192 = 253250. - 253250.000000000 - 1.000000E-07 * -2.941266E+08 34. ET(2) = xx(5,2) : -848141. = -848141.388090847 -------------------------------------------------------------------------- Petsc Release Version 2.3.0, Patch 32, April, 26, 2005 See docs/changes/index.html for recent updates. See docs/faq.html for hints about trouble shooting. See docs/index.html for manual pages. ----------------------------------------------------------------------- ./SNESG3D2 on a linux-gnu named sethnx004.vac.com by yy26539 Mon Feb 27 14:21:13 2006 Libraries linked from /home/yy26539/work/NISSE/petsc-2.3.0/lib/linux-gnu Configure run at Thu Oct 13 08:23:56 2005 Configure options --with-cc=gcc --with-fc="f77 -N109" --download-mpich=1 --download-mpich-pm=gforker --download-f-blas-lapack=1 --with-shared=0 ----------------------------------------------------------------------- [0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [0]PETSC ERROR: likely location of problem given in stack below [0]PETSC ERROR: --------------- Stack Frames --------------- [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available, [0]PETSC ERROR: INSTEAD the line number of the start of the function [0]PETSC ERROR: is given. [0]PETSC ERROR: [0] SNES user function line 0 unknownunknown [0]PETSC ERROR: [0] SNESComputeFunction line 788 src/snes/interface/snes.c [0]PETSC ERROR: [0] SNESLineSearchCubic line 514 src/snes/impls/ls/ls.c [0]PETSC ERROR: -------------------------------------------- [0]PETSC ERROR: User provided function() line 0 in unknown directory unknown file [0]PETSC ERROR: Signal received! [0]PETSC ERROR: ! Suggestions and tips are very welcome! /nisse -- Nils-Erik Svang�rd MSN: schweingaard@hotmail.com Skype: schweingaard-- Nils-Erik Svangård MSN: schweingaard@hotmail.com Skype: schweingaard-- Nils-Erik Svangård E-Mail: nilserik@gmail.com MSN: schweingaard@hotmail.com Skype: schweingaard Mobil: +46-(0)70-3612178-- Nils-Erik Svangård E-Mail: nilserik@gmail.com MSN: schweingaard@hotmail.com Skype: schweingaard Mobil: +46-(0)70-3612178-- Nils-Erik Svangård E-Mail: nilserik@gmail.com MSN: schweingaard@hotmail.com Skype: schweingaard Mobil: +46-(0)70-3612178