You're absolutely right, i missed the y-direction component of the momentum equation in the mass matrix integration.
I guess i have most of the parts i need from here to investigate the remaining on my own! Thank you so much! On Mon, 14 Nov 2016 07:53:45 -0600 Matthew Knepley <[email protected]> wrote: > 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 > > > > >
