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

[Roy S. <roy...@ic...>]

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