Discussion Forum

Posted: 1 decade ago 22 mai 2013, 05:54 UTC−4

Hello Daniel, There are a few sections about constraints in the Equation-Based Modeling chapter in the version 4.3b Reference Manual. These may give you some clues as to how the weak constraints are implemented, even if it is not stated explicitly. The section on Boundary Conditions in the Building a COMSOL Model chapter may also provide useful background information. If you understand how the standard pointwise constraints are incorporated into the system of equations, the jump from there to weak constraints is not very long: when using weak constraints, the Lagrange multipliers are included in the system of equations using standard shape functions, which are also used as test functions multiplying the constraint condition. The way I use to think about weak constraints (and constraints in general) is roughly as follows: Say that you want to constrain N degrees of freedom. This is equivalent to adding N extra equations to the total system (rows in the matrix). However, since the system was square before adding the constraints and must stay square in order to be solvable, you must also add N additional degrees of freedom (columns in the matrix) - the Lagrange multipliers. Now the constraint force rows may be assembled either as constraints directly on each degree of freedom (standard pointwise constraints) or by multiplying the constraint by a test function and integrating (weak constraints). The former case is actually a special case of the latter if we allow for test functions of Dirac distribution type. Since the number of constraint equations is, by definition, equal to the number of Lagrange multiplier DOFs necessary to keep the system square, it makes sense to identify the test functions multiplying the constraint equations with the Lagrange multipliers. This, however, does not in itself uniquely determine the contents of the columns that must be added to return the system to square form. In fact, you are rather free to choose the form of this "constraint force Jacobian" block of the overall system - as long as you make sure that the complete system is non-singular. COMSOL actually lets you specify a constraint force expression manually, if you chose to "Apply reaction terms on" "User defined" in the Weak Constraint feature. The default behavior in COMSOL is, however, to apply constraint reaction terms in the unique way which retains symmetry in the system matrix when the unconstrained system is symmetric. This choice (because it is indeed a choice) is related both to conservation of energy, as explained in http://www.softeng.rl.ac.uk/st/projects/felib3/Docs/html/Intro/intro-node68.html, to the KKT conditions for a constrained local minimum in an optimization problem and to the principle of symmetric action and reaction expressed in, e.g., Newton's third law. When comparing to http://www.softeng.rl.ac.uk/st/projects/felib3/Docs/html/Intro/intro-node68.html, note that there is a transpose missing on the upper right "B" component in the system structure stated and that there is a strange row shift (probably a mistake) in the matrix example. When using Lagrange multipliers to compute accurate fluxes or forces, note that the true "constraint force" is the product of the Lagrange multipliers and the constraint force Jacobian. Therefore, the Lagrange multiplier itself is only the accurate flux if the corresponding constraint force Jacobian column sums to 1, which for a symmetric constraint means that the coefficient in front of the dependent variable being constrained must be equal to 1. For example, if you change the Weak Constraint in your example model to "2*u-2", which is completely equivalent from the constraint point of view and gives the same solution for u, the Lagrange multiplier value will be reduced by a factor 2 compared to the original model. The reason for this is that the constraint force Jacobian used now contains "2" instead of "1" in its only non-zero position, so the product of Jacobian and Lagrange multiplier is still the same. Finally, I don't quite understand your remark: [QUOTE] However, when I do this, I get the same answer for flux using either method (directly elimination or lagrange multipliers) once I solve the system of equations. BUT, when I do the same thing in COMSOL with the same problem definition, same constraints, and same number of elements, same basis function order, etc., the lagrange multiplier value I get matches up perfectly with the analytical solution. [/QUOTE] The model you attached is not really a good example for experimenting with weak constraints and comparing different flux evaluation strategies. The reason is essentially that when both shape functions and exact solution are linear, both the flux evaluated from the gradient of u at the boundary and averaged over the nearest element will coincide with the exact expected flux. To get a more interesting situation where the Lagrange multiplier becomes more accurate than direct evaluation, try adding a non-uniform source term in the domain such that the solution is no longer linear. best regards Nils Malm

Posted: 1 decade ago 22 mai 2013, 20:20 UTC−4

Thanks for the detailed response Nils. I think I have a pretty good understanding now of what's going on with the Lagrange Multipliers. Sorry I uploaded the wrong version of the file. The one I was referring to has a source term of 10J/sm (f=-10 in COMSOL) (an example taken from Zimmerman's book "Process Modelling and Simulation with Finite Element Methods" that COMSOL lists on their book references page). When I write out the matrix by hand and do the matrix calculations based on the link provided before, I get a lagrange value of 8.65. This is the same value that you get if you let temperature at node 1 be a constant in the system of equations, and rearrange them for the unknown flux at the dirichlet node (making it no longer symmetric, but that doesn't matter for such a small matrix), which is what you would expect - since when you use the lagrange approach you assume a test function value of zero at dirichlet boundaries and the lagrange multiplier takes up the necessary flux to maintain a constant temperature at that node as far as I understand. However, my quote that you referenced there was to point out that if I set up this same problem that I do on paper, in COMSOL, I get the flux value of 11.25 when using lagrange multipliers, which matches up with the analytical solution (which is quadratic). But this doesn't match up to what I get on paper so COMSOL must be doing something different/additional that I'm not understanding. EDIT: It's possible I'm setting it up wrong of course!

Posted: 1 decade ago 24 mai 2013, 13:36 UTC−4

Was a calculation on my part, everything matches up now. Thanks for the detailed explanation, I have a pretty good idea of what's going on now!