On Thu, 3 Mar 2011, Boyce Griffith wrote: > On 3/2/11 11:55 AM, Roy Stogner wrote: > >> For the inhomogeneous cases: we basically turn K*x = F, x=C*y into >> C^T*K*C*y=C^T*F. That's still a symmetric matrix. It's a smaller >> symmetric matrix in the case where length(y) < length(x), but we don't >> actually condense the matrix, so we end up having indeterminate >> hanging node entries in y. In the symmetric case we make the matrix a >> diagonal for those entries, and then we fix them up later after the >> linear solve with enforce_constraints_exactly. In the asymmetric case >> we reuse those entries' lines with the constraint equations so that >> the solver will enforce them (if only inexactly). >> >> So now we've got x=C*y+H >> >> Which means if we do things analogously in the symmetric case we want >> >> C^T*K*(C*y+H)=C^T*F >> C^T*K*C*y = C^T*(F-K*H) >> >> Hmm... this isn't what your patch does, and when I think about it we >> probably *don't* want to do things analogously (it's nice that in the >> homogeneous case x is also just y with some entries omitted, but the >> "x=C*y+H" model doesn't preserve that property). Is there any >> similarly clean way to formulate what your patch is doing? > > I'll try to think about this next week. I will also need to review what the > patch actually does. :-)
After running into a problem where penalty boundary constraints were quite poor, I'm renewing my interest in fixing this. I've got an updated version of your patch locally (libMesh has diverged too much from what you diffed against to apply the original patch without updates, thanks mostly to my NodeConstraints stuff) and I've been thinking about the issues. Let's expand last year's notation by assuming we've got a nonlinear problem f(x) = 0 to solve, and we know that x is a function of the *real* independent variables y; i.e. x = m(y) gives us x in some higher dimensional space. We'll have derivatives K := df/dx (our Jacobian) and C := dx/dy (our constraint matrix) We'll use linearizations f(x) := K*x-F and m(y)=C*y+H to double-check against intuition, but it's critical to make sure we've got something that works for nonlinear physics and preferably even nonlinear Dirichlet conditions. Imagine an isothermal boundary in a conserved-variables flow code, to pick a non-hypothetical example of the latter. Ben, I'd appreciate your corrections here, since you've looked into that example much more deeply than I have; in particular there must be a nice library-encapsulatable way to handle the implicit nonlinear Dirichlet constraints case without turning it into an explicit local Newton solve m() or adding Lagrange multiplier variables? But ignoring that for the moment: Instead of solving a weighted problem (f(x), dx) = 0 \forall dx, we want to solve some (g(y), dy) = 0 \forall dy. dx = dx/dy * dy = C*dy, so we can turn our original weighted-f equations into (f(m(y)),C*dy) = 0, so (g(y),dy) := (C^T*f(m(y)),dy) = 0. Plug in linear homogeneous constraint equations and that gives us residual g(y) = C^T*f(C*y) and jacobian dg/dy = C^T*K*C which is just what we were solving before. Reassuring. For nonlinear problems, we have to do the x = m(y) (by hand with enforce_constraints_exactly()?) before evaluating each residual and/or Jacobian to get the nonlinearities correct, do another enforce_constraints_exactly() at the end of solves for consistency in the inexact solve and nonlinear constraints cases, and then constrain_element_matrix does C^T*K*C just as before, and constrain_element_vector does C^T*F just as before. Hooray! For linear problems, we don't have to worry about the input x, constrain_element_matrix does C^T*K*C just as before, and constrain_element_vector does C^T*(F-K*H) to get the output heterogeneously constrained variables correct. Hooray??? So we basically need *two* constrain_element_vector (and constrain_element_matrix_and_vector) versions: one to do exactly what we did before, and one to handle cases where a single linear solve for a heterogeneous problem is being done rather than a series of linear solves for homogeneous updates. Old user apps for linear problems will need to change their assembly routines to handle true Dirichlet boundaries. I don't see that as a concern, since old user apps would need to change their assembly and initialization routines anyway just to *specify* true Dirichlet boundaries. Anything I'm missing? I'm going to start coding and testing everything this weekend, bearing in mind that I may tear it all up and start over if Ben's got better ideas based on his recent FIN-S true-Dirichlet experience. --- Roy ------------------------------------------------------------------------------ Virtualization & Cloud Management Using Capacity Planning Cloud computing makes use of virtualization - but cloud computing also focuses on allowing computing to be delivered as a service. http://www.accelacomm.com/jaw/sfnl/114/51521223/ _______________________________________________ Libmesh-devel mailing list Libmesh-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-devel
