Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

## Re: [Libmesh-devel] [Libmesh-users] interaction between subdomain_id and dof constraints?

 Re: [Libmesh-devel] [Libmesh-users] interaction between subdomain_id and dof constraints? From: Roy Stogner - 2012-02-25 00:53:13 ```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 ```

 Re: [Libmesh-devel] [Libmesh-users] interaction between subdomain_id and dof constraints? From: Roy Stogner - 2011-03-02 16:55:49 ```On Wed, 2 Mar 2011, Boyce Griffith wrote: > On 3/1/11 6:41 PM, Roy Stogner wrote: >> >> On Thu, 2 Dec 2010, Boyce Griffith wrote: >> >>> The patched code does the "right thing" only when >>> asymmetric_constraint_rows=true when calling >>> DofMap::constrain_element_{matrix,vector,matrix_and_vector}(). This is >>> the default, but if asymmetric_constraint_rows=false, the code will >>> silently do the "wrong thing," probably yielding bizarre results. One >>> approach would be to keep track of whether a particular constraint row >>> is provided by the user or not, and to require >>> asymmetric_constraint_rows==true when such constraints are present. >> >> Thanks for the heads-up. I've got some FEMSystem-based codes that >> could probably be adjusted pretty easily to test that case, which we'd >> want to fix before committing. > > Let me know if you want me to think about what to do for this case, or if it > isn't clear what goes wrong if asymmetric_constraint_rows=false. I haven't looked at the problem too closely yet, but it might be better to think about the math before the code. 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? --- Roy ```
 Re: [Libmesh-devel] [Libmesh-users] interaction between subdomain_id and dof constraints? From: Roy Stogner - 2012-02-25 00:53:13 ```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 ```
 Re: [Libmesh-devel] [Libmesh-users] interaction between subdomain_id and dof constraints? From: Roy Stogner - 2012-02-27 22:59:08 Attachments: heterogenous.diff ```On Mon, 27 Feb 2012, Boyce Griffith wrote: > On 2/24/12 7:53 PM, Roy Stogner wrote: > >> 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. > > Let me know if I can help. I have some time this week that I could spend > working on this. That would be great, since I don't have time to work on this. But why let that stop me? Updated patch is posted; thanks also to Derek for offering to test. Assuming this doesn't break anyone's adaptive, periodic, etc. codes, then I'll commit and add a Dirichlet constraints generation routine (probably an API along the lines of the functor version of System::project_boundary_solution unless anyone dislikes that) this week. --- Roy```