Hi Daniel, The vectors you have in this equation P(u+h*v) =..., which of those have constrains distributed and which zeroed? If you assembly matrices with ConstraintMatrix.distribute_local_to_global() the diagonal elements corresponding to constrained DoFs will be dummy (and positive) just to make SLAE positive-definite. When you use such solution further in matrix-vector products you would most likely need to set constrained DoFs of the solution vector to zero because your matrix is assembled with constrained already distributed. Say in my case I need to compute the scalar (u, Ku) and that's where i would zero constrained DoFs. The same value could be obtained from (u', K' u') where by ' i denote object which don't take into account constrains, i.e. solution vector with distribute constrains and matrix assembled with empty constrains. On another hand, if you use this vector to evaluate the unknown field and update quadrature data, then you would need to distribute constrains to make the field continuous.
Along the same lines: does you test work with non-homogeneous Dirichlet BC? The discussion above applies there as well. Regards, Denis. On Thursday, February 23, 2017 at 1:49:08 AM UTC+1, Daniel Shapero wrote: > > Hi all -- I'm writing a library that involves solving a nonlinear elliptic > PDE, which I'll write as f(u) = 0. There is an exact solution for this PDE > for a certain simplified geometry. To test everything, I check that my > numerical solution is tolerably close to the analytic solution, with both > bilinear and biquadratic elements, and on a mesh that has been adaptively > refined in part of the domain. This all works fine. > > The PDE can be derived from an action principle, i.e. there's a convex > functional P such that f = dP. I've decided to reimplement some of my code > to explicitly use the action principle, so I wrote a routine to calculate > P. The action is analytically computable for the exact solution I'm already > using to test the code, so this is already one unit test. However, I > decided to take it a step further, and check that the nonlinear > differential operator f is the derivative of P, and that its linearization > df is the Hessian operator of P. To check that, I take some vector field u, > another vector field v, and check that > > P(u + h * v) = P(u) + h*f(u) . v + h^2 * (df(u)v) . v / 2 + O(h^3) > > This test passes in the simple geometry with both bilinear and biquadratic > elements, but *not* for an adaptively refined mesh. The errors in both the > local linear approximation to the action functional and the local quadratic > approximation do not go to 0 as the increment h is reduced. In addition, > the value of the action is roughly the same for each h for both the > adaptively refined and uniform meshes. The local quadratic approximation > also converges as h goes to 0 on a realistic, complex geometry with real > data, although we don't have an analytic solution for the problem there. > > I'm very confused by this result; since the solution of the PDE on an > adaptively-refined mesh agrees with the analytic solution, it's kind of > baffling that the local quadratic approximation to the action functional > isn't working right in that case. The only other clue I have is that the > error of the numerical solution against the analytic solution is actually > somewhat worse on the adaptively-refined mesh than for the uniform mesh, > but before I assumed this just had something to do with the linear solver. > > Any ideas of what I might have missed? The code is available here > <https://github.com/danshapero/icepack/tree/action-functional> if that > helps. > > Daniel > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
