Note: This discussion is about an older version of the COMSOL Multiphysics® software. The information provided may be out of date.
Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.
Incorporating Lagrange Multipliers
Posted 21 mag 2013, 20:52 GMT-4 Results & Visualization, Studies & Solvers Version 4.0a 3 Replies
Please login with a confirmed email address before reporting spam
Hello again,
Thanks to the help of those on the forum I was able to incorporate weak constraints to get accurate flux calculations in my diffusion models. However, I am trying to understand how exactly the lagrange multipliers work in terms of the final matrix construction and there is very little documentation in the comsol manuals and guides. I tried looking up some papers but they are written at a very high level of math that I'm not able to grasp the details of.
I tried writing out a simple heat diffusion problem in 1D with 3 mesh elements and solving with direct elimination, and then incorporating the lagrange multiplier via the approach shown here: www.softeng.rl.ac.uk/st/projects/felib3/Docs/html/Intro/intro-node68.html
I just wrote out the weak form and calculated the matrix by hand which isn't very hard for a 3 element 1D problem.
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. I would like to know how this is implemented, so that I can understand what is going on when I choose these constraints.
If anyone is able to explain it to me or provide me with a link I would be very grateful!
Thanks to the help of those on the forum I was able to incorporate weak constraints to get accurate flux calculations in my diffusion models. However, I am trying to understand how exactly the lagrange multipliers work in terms of the final matrix construction and there is very little documentation in the comsol manuals and guides. I tried looking up some papers but they are written at a very high level of math that I'm not able to grasp the details of.
I tried writing out a simple heat diffusion problem in 1D with 3 mesh elements and solving with direct elimination, and then incorporating the lagrange multiplier via the approach shown here: www.softeng.rl.ac.uk/st/projects/felib3/Docs/html/Intro/intro-node68.html
I just wrote out the weak form and calculated the matrix by hand which isn't very hard for a 3 element 1D problem.
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. I would like to know how this is implemented, so that I can understand what is going on when I choose these constraints.
If anyone is able to explain it to me or provide me with a link I would be very grateful!
Attachments:
3 Replies Last Post 24 mag 2013, 13:36 GMT-4
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
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 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 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:
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
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 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 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:
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.
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
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
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!
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!
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
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!