Hi, Thanks Matthew and Jed for your input. I indeed envision an implicit solver in the sense Jed mentioned - Jiri Blazek's book is a nice intro to this concept.
Matthew, I do not know exactly what to change right now because although I understand globally what the DMPlexComputeXXXX_Internal methods do, I cannot say for sure line by line what is happening. In a structured code, I have a an implicit FVM solver with PETSc but I do not use any of the FV structure, not even a DM - I just use C arrays that I transform to PETSc Vec and Mat and build my IJacobian and my preconditioner and gives all that to a TS and it runs. I cannot figure out how to do it with the FV and the DM and all the underlying "shortcuts" that I want to use. Here is the top method for the structured code : int total_size = context.npoints * solver->nvars ierr = TSSetRHSFunction(ts,PETSC_NULL,PetscRHSFunctionImpl,&context); CHKERRQ(ierr); SNES snes; KSP ksp; PC pc; SNESType snestype; ierr = TSGetSNES(ts,&snes); CHKERRQ(ierr); ierr = SNESGetType(snes,&snestype); CHKERRQ(ierr); flag_mat_a = 1; ierr = MatCreateShell(MPI_COMM_WORLD,total_size,total_size,PETSC_DETERMINE, PETSC_DETERMINE,&context,&A); CHKERRQ(ierr); context.jfnk_eps = 1e-7; ierr = PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL); CHKERRQ(ierr); ierr = MatShellSetOperation(A,MATOP_MULT,(void (*)(void)) PetscJacobianFunction_JFNK); CHKERRQ(ierr); ierr = MatSetUp(A); CHKERRQ(ierr); context.flag_use_precon = 0; ierr = PetscOptionsGetBool(PETSC_NULL,PETSC_NULL,"-with_pc",(PetscBool*)(& context.flag_use_precon),PETSC_NULL); CHKERRQ(ierr); /* Set up preconditioner matrix */ flag_mat_b = 1; ierr = MatCreateAIJ(MPI_COMM_WORLD,total_size,total_size,PETSC_DETERMINE, PETSC_DETERMINE, (solver->ndims*2+1)*solver->nvars,NULL, 2*solver->ndims*solver->nvars,NULL,&B); CHKERRQ(ierr); ierr = MatSetBlockSize(B,solver->nvars); /* Set the RHSJacobian function for TS */ ierr = TSSetIJacobian(ts,A,B,PetscIJacobian,&context); CHKERRQ(ierr); Thibault Bridel-Bertomeu — Eng, MSc, PhD Research Engineer CEA/CESTA 33114 LE BARP Tel.: (+33)557046924 Mob.: (+33)611025322 Mail: thibault.bridelberto...@gmail.com Le jeu. 20 août 2020 à 18:43, Jed Brown <j...@jedbrown.org> a écrit : > Matthew Knepley <knep...@gmail.com> writes: > > > I could never get the FVM stuff to make sense to me for implicit methods. > > Here is my problem understanding. If you have an FVM method, it decides > > to move "stuff" from one cell to its neighboring cells depending on the > > solution to the Riemann problem on each face, which computed the flux. > This > > is > > fine unless the timestep is so big that material can flow through into > the > > cells beyond the neighbor. Then I should have considered the effect of > the > > Riemann problem for those interfaces. That would be in the Jacobian, but > I > > don't know how to compute that Jacobian. I guess you could do everything > > matrix-free, but without a preconditioner it seems hard. > > So long as we're using method of lines, the flux is just instantaneous > flux, not integrated over some time step. It has the same meaning for > implicit and explicit. > > An explicit method would be unstable if you took such a large time step > (CFL) and an implicit method will not simultaneously be SSP and higher than > first order, but it's still a consistent discretization of the problem. > > It's common (done in FUN3D and others) to precondition with a first-order > method, where gradient reconstruction/limiting is skipped. That's what I'd > recommend because limiting creates nasty nonlinearities and the resulting > discretizations lack h-ellipticity which makes them very hard to solve. >
