On Mon, Nov 14, 2016 at 7:35 AM, Matthew Knepley <[email protected]> wrote:
> On Mon, Nov 14, 2016 at 3:24 AM, Julian Andrej <[email protected]> > wrote: > >> I think i'm understanding the concept now. Although i have another set >> of questions. If i am formulating a DAE (like in the linear stokes >> equation for example) it seems i am missing something to formulate the >> correct jacobian. The residual formulation with snes_mf or snes_fd >> works fine in terms of nonlinear convergence. If i use the hand coded >> jacobian the convergence is way slower, which is a reason for an >> incorrect jacobian formulation from my side (at least what i think >> right now). I should be able to test my jacobian with -snes_type test, >> even if i am using TS, right? This gives me a huge difference. I guess >> i am missing something in the formulation and thus in the >> implementation of the callbacks here. >> >> There is a small example attached. >> > > There are a few problems here, but first here are my options for the run > > -vel_petscspace_order 2 -vel_petscspace_poly_tensor 1 > -pres_petscspace_order 1 -pres_petscspace_poly_tensor 1 > -ts_dt 0.0001 -ts_max_steps 5 > -pc_type fieldsplit -pc_fieldsplit_type schur > -pc_fieldsplit_schur_precondition > full > -fieldsplit_velocity_pc_type lu > -fieldsplit_pressure_ksp_rtol 1.0e-9 -fieldsplit_pressure_pc_type svd > -ts_monitor -snes_monitor -snes_view -snes_converged_reason > -ksp_converged_reason -ksp_monitor > > 1) Your g00 term was wrong since it did not handle all velocity components > > void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, > const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], > const PetscScalar u_t[], const PetscScalar u_x[], > const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], > const PetscScalar a_t[], const PetscScalar a_x[], > PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[]) > { > PetscInt d; > for (d = 0; d < dim; ++d) g0[d*dim+d] = u_tShift*1.0; > } > > 2) You do not handle the null space for the pressure solve right now. The > easy way is to use SVD for the PC, however this > is not scalable. Eventually, you want to specify the constant null > space for that solve. I do it in ex62. > > 3) I still do not know what is wrong with the exact solution... > I found it. You are not making the initial condition. Uncomment the line above TSSolve(). Thanks, Matt > Thanks, > > Matt > > >> I appreciate your time and help. >> >> On Thu, Nov 10, 2016 at 4:45 PM, Matthew Knepley <[email protected]> >> wrote: >> > On Thu, Nov 10, 2016 at 9:38 AM, Jed Brown <[email protected]> wrote: >> >> >> >> Matthew Knepley <[email protected]> writes: >> >> >> >> > On Thu, Nov 10, 2016 at 1:31 AM, Julian Andrej <[email protected]> >> >> > wrote: >> >> > >> >> >> I'm getting the correct solution now, thanks! There are still a few >> >> >> open questions. >> >> >> >> >> >> 1. The mass term is necessary _and_ using the u_tShift value if i >> >> >> provide the jacobian. After reading the manual countless times, i >> >> >> still don't get what the u_tShift value tells me in this context. >> Are >> >> >> there any resources which you could point me to? >> >> >> >> >> > >> >> > The manual is always a lagging resource because we are trying things >> >> > out. >> >> >> >> This has been explained in the manual similarly to your email for >> >> several years. If anyone has a suggestion for how to make it better, >> >> we'd like to hear. The typical syndrome is that once you learn it, the >> >> description suddenly makes sense and you wonder why you were ever >> >> confused. It takes feedback from people just learning the material to >> >> make the explanation more clear. >> >> >> >> > I learned about this by reading examples. The idea is the >> >> > following. We have an implicit definition of our timestepping method >> >> > >> >> > F(u, grad u, u_t, x, t) = 0 >> >> >> >> It doesn't make sense to write grad u as a separate argument here. >> > >> > >> > Sorry, that is just how I think of it, not the actual interface. >> > >> >> >> >> Also, is `x` meant to be coordinates or some other statically >> prescribed >> > >> > >> > Coordinates, as distinct from u. Again this is my fault. >> > >> >> >> >> data? It can't be actively changing if you expect a generic TS to be >> >> convergent. (It's not an argument to the function.) So really, you >> >> have: >> >> >> >> F(u, u_t, t) = 0 >> >> >> >> > which is not unlike a Lagrangian description of a mechanical system. >> The >> >> > Jacobian can be considered >> >> > formally to have two parts, >> >> > >> >> > dF/du and dF/du_t >> >> > >> >> > just as in the Lagrangian setting. The u_tShift variable is the >> >> > multiplier >> >> > for dF/du_t so that you actually >> >> > form >> >> > >> >> > dF/du + u_tShift dF/du_t >> >> >> >> We're taking the total derivative of F with respect to u where u_t has >> >> been _defined_ as an affine function of u, i.e., shift*u+affine. >> >> u_tShift comes from the chain rule. When computing the total >> >> derivative, we have >> >> >> >> dF/du + dF/du_t du_t/du. >> >> >> >> >> >> >> 2. Is DMProjectFunction necessary _before_ TSSolve? This acts like >> an >> >> >> initial condition if i'm not mistaken. >> >> >> >> >> > >> >> > Yes, this is the initial condition. >> >> > >> >> > >> >> >> 3. A more in depth question i guess. Where is the "mass action" >> >> >> formed? As i could see by stepping through the debugger, there is >> just >> >> >> a formed action for the residual equation. It seems way over my >> >> >> understanding how this formulation works out. I'd also appreciate >> >> >> resources on that if thats possible. >> >> >> >> >> > >> >> > Jed is the only one who understands this completely. >> >> >> >> That's a stretch -- several other people have developed sophisticated >> TS >> >> implementations. >> > >> > >> > Matt does not understand this completely. >> > >> >> >> >> >> >> > However, I guess by "mass action" you mean the dF/du_t term. In the >> >> > explicit methods, you just have u_t + ..., so >> >> > >> >> > dF/du_t = M >> >> > >> >> > for finite element methods. So you are putting that there in >> >> > FormIJacobian(). In the simplest >> >> > case of backwards Euler then u_tShift would be 1/dt. >> >> >> >> Specifically, for implicit methods, and many semi-implicit methods, we >> >> don't need M separate. >> > >> > >> > Yes. >> > >> > Matt >> > >> > -- >> > What most experimenters take for granted before they begin their >> experiments >> > is infinitely more interesting than any results to which their >> experiments >> > lead. >> > -- Norbert Wiener >> > > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
