Implementing the Weak Form with a COMSOL App

Previously in our weak form series, we discretized the weak form equation to obtain a matrix equation to solve for the unknown coefficients in our simple example problem. Following the same procedure as in this previous blog post, we will implement the equation in the COMSOL Multiphysics® software with additional steps included to examine the matrices. We will find it more convenient to use a COMSOL® software application to display all relevant matrices at once, arranged logically on one screen.

Our Simple Example

Recall our simple example of 1D heat transfer at steady state with no heat source, where the temperature T is a function of the position x in the domain defined by the interval 1 \le x \le 5. We imposed the Neumann boundary condition such that the outgoing flux should be 2 at the left boundary (x=1) and the Dirichlet boundary condition such that the temperature should be 9 at the right boundary (x=5). After discretizing the weak form equation, we obtained this matrix equation:

(1)

where a_1, a_2, \cdots , a_5 are unknown temperature values at the nodal points (x=1, 2, \cdots, 5) and \lambda_2 is an unknown heat flux at the right boundary (x=5). The matrix on the left-hand side is customarily called the stiffness matrix and the vector on the right is called the load vector, due to the application of this technique in structural mechanics.

Viewing the Stiffness Matrix and Load Vector

The steps for implementing the weak form equation in COMSOL Multiphysics have been discussed in this earlier blog entry, thus we will not repeat them here. To view the stiffness matrix and load vector, right-click Solution 1 in the model tree → Other → Assemble, as shown in the screenshot below:

Make sure to drag the Assemble node up to position it above (before) the Stationary Solver node. If this step is not done and the Assemble node is left below (after) the Stationary Solver node, then the content of the load vector may be altered by the solver.

Then, in the settings window for the Assemble 1 node, we can select the matrices of interest by checking the corresponding checkbox for each item. Additionally, set Scaling of Dependent Variables 1 to None. After solving, we can evaluate the matrices by right-clicking Derived Values and then selecting System Matrix, as illustrated below:

In the corresponding settings window, we can select Stiffness matrix from the Matrix drop-down menu and evaluate it in a table. As indicated in the screenshot below, we obtain exactly the same 6×6 matrix as the one in Eq. (1).

The load vector on the right-hand side of Eq. (1) can be evaluated and verified by the same procedure.

Pointwise Constraint

We mentioned before that it is sometimes desirable not to solve for the Lagrange multiplier \lambda_2 — for example, to save computation resources. To do so, we right-click the Weak Form PDE main node → More → Pointwise Constraint, as depicted below:

Then, in the settings window, we enter 9-T for the Constraint expression input field. After solving, we obtain exactly the same solution, but a smaller 5×5 stiffness matrix:

If we look at this matrix closely, we find that it matches the upper left-hand part of the 6×6 matrix that we observed earlier. This should not be too surprising, as we are solving exactly the same problem (represented by Eq. (1)), just in a slightly different way. We will briefly discuss this next.

Alternative Methods for Solving the Matrix Equation

When we implement the Dirichlet boundary condition using the Weak Contribution feature with a Lagrange multiplier \lambda_2 (as shown in this previous blog post), we effectively ask that the entire matrix equation (1) be solved to yield the coefficients a_1, a_2, \cdots , a_5 as well as the multiplier \lambda_2.

On the other hand, when we implement the fixed boundary condition using the Pointwise Constraint feature as shown above, we essentially ask that the same matrix equation (1) be solved without explicitly solving for the Lagrange multiplier \lambda_2. The software effectively segregates Eq. (1) to this form (see below)

(2)

where the stiffness matrix K is the 5×5 matrix shown above, the constraint Jacobian matrix N is (0 \, 0 \, 0 \, 0 \, 1), the constraint force Jacobian matrix N_F is (0 \, 0 \, 0 \, 0 \, 1)^T, the solution vector U is formed by the coefficients a_1, a_2, \cdots , a_5, the (one-element) vector \Lambda is formed by the Lagrange multiplier \lambda_2, the load vector L is (-2 \, 0 \, 0\, 0\, 0)^T, and the (one-element) constraint vector M is (9).

This segregation of Eq. (1) is shown graphically below:

Of course, in reality, when the Pointwise Constraint feature is specified, the software does not bother to assemble the full 6×6 matrix in Eq. (1). Instead, it assembles K, N, and N_F, which requires less computational resources than assembling the full matrix.

Eq. (2) can be written as a system of two matrix equations:

(3)

(4)

where the first one (3) contains the part of the discretized weak form equation with heat flux boundary conditions and the second one (4) contains the constraint imposed by the Dirichlet boundary condition.

The equation system can be solved in two steps. In the first step, the constraint equation (4) is solved for the degree(s) of freedom involved with the Dirichlet boundary condition. In the second step, the solution from the first step is substituted into Eq. (3) to solve for the remaining degrees of freedom.

The resulting solution vector U is then given as the linear combination of the solutions from the two steps:

(5)

The first term U_d is the solution vector of the constraint equation (4) solved in the first step. The second term is obtained from the second step in the form of the product of a matrix Null and a solution vector U_n. The columns of the matrix Null are formed by the basis vectors spanning the null space of the constraint Jacobian matrix N. So, we have

The solution vector U_n is the solution to the eliminated matrix equation

(6)

where the eliminated stiffness matrix K_c and the eliminated load vector L_c are given by

and the columns of the matrix Nullf are formed by the basis vectors spanning the null space of the transpose of the constraint force Jacobian matrix N_F. So, we have

The term eliminated indicates that the degree(s) of freedom involved in the Dirichlet boundary condition have been removed from the matrix equation (6). In our example, the size of the eliminated stiffness matrix K_c is 4×4.

Notice that the solution procedure described above does not involve the Lagrange multiplier vector \Lambda. Indeed, the value of the Lagrange multiplier \lambda_2 is left unsolved by this procedure. The advantage of this method is that the required computational resources are reduced. For our simple example, the size of the matrix is reduced from 6×6 to 4×4 (plus an even smaller one for the constraint equation).

Viewing All Matrices at Once with a COMSOL® App

COMSOL Multiphysics allows us to evaluate and view any of the matrices and vectors mentioned above. All we have to do is to follow the same procedure outlined in the previous section on viewing the stiffness matrix and load vector. However, we quickly find out that we will spend a lot of time clicking on each System Matrix node in the Model Builder and then clicking on the “Evaluate” button in the settings window. Also, we can only see one matrix or vector at a time, making it tedious to examine all the matrices and vectors.

This represents one of the many situations in which a COMSOL application can help enhance modeling experience and productivity (learn about application webinars below). The core package of COMSOL Multiphysics can be used to build a COMSOL app based on any COMSOL model by wrapping a user interface (UI) around it. The UI is completely customizable and can be easily configured to suit individual modeling needs.

For an introduction to COMSOL applications, see the following webinars: Intro to COMSOL Multiphysics® 5.0 and the Application Builder (with a focus on the second half) and How to Build and Run Simulation Apps with COMSOL Server™.

The following screenshots show just one of the essentially infinite number of ways that a UI can be arranged to serve as a convenient tool for investigating the matrices and vectors. The app is set up to switch among several different kinds of boundary conditions via checkboxes. All matrices and vectors are evaluated and displayed by a single click of the “Compute” button. Here is a screenshot of the case in which the full 6×6 matrix is used to solve:

We see that N, N_F, and M are empty and K_c remains the full size of 6×6. In contrast, here is the screenshot for the second case where the eliminated 4×4 matrix is used to solve:

Concluding Remarks

Using our simple heat transfer example, we have explored the two different ways to solve the matrix equation obtained from discretizing the weak form equation. We found that it can be much simpler to evaluate and inspect the various matrices and vectors involved in the solution process by setting up a COMSOL application. In the next blog posts, we will continue to use COMSOL apps to help us investigate more complex examples.

As a general-purpose development environment for multiphysics simulation, the COMSOL Desktop® environment must be structured with a limited arrangement of UI elements. The COMSOL app frees us from such limitations, providing the power for building specialized UI for unique requirements. Additionally, the built-in coding functionality allows much more versatile computations, and COMSOL Server™ license enables the deployment of apps for worldwide access by coworkers as well as clients.

We hope that you will use these powerful features in COMSOL Multiphysics to boost your own work efficiency!