Matthias K. Gobbert

# Introduction to COMSOL Multiphysics

## Purpose of this Document

COMSOL Multiphysics is an excellent, state-of-the-art software for the solution of many types of partial differential equations (PDEs), both stationary and time-dependent, by numerical techniques based on the finite element method for the spatial discretization.

• Provide step-by-step instructions how to solve one simple PDE using the graphical user interface (GUI) of COMSOL Multiphysics. This example should be useful to get a first idea of how COMSOL's GUI works and can also be used by users (or system administrators) to test whether a given installation of COMSOL works properly.
• The best way to get to the extensive documentation of COMSOL is to go through the graphical user interface (GUI) of COMSOL Multiphysics. Hence, the instructions for the small example are also useful to get into the GUI and to access the documentation from there.
• I have a few files that users of COMSOL might find useful, and I am using this page as a repository of these files. Look near the end of the page.

## Solution of a Simple PDE using the GUI

#### How to Start COMSOL Multiphysics:

Under Windows, double click on the COMSOL desktop icon or look for COMSOL under the Start menu. (This should be similar under Mac OS X.) Under Linux, you can enter comsol -h to see a list of how to start COMSOL; for the GUI version of COMSOL Multiphysics, enter comsol multiphysics or the short-hand form comsol at the Linux command-line. (To run it in the background and thus get your shell back in the terminal window, add an ampersand at the end of either command, such as comsol&.)

#### Start the GUI by Selecting the Mode in the Model Navigator:

COMSOL will bring up a Model Navigator window now. The space dimension is preselected as 2D. Under Application Modes, expand (by clicking on the plus sign) COMSOL Multiphysics, then PDE Modes, then PDE Coefficient Form, and select Stationary analysis. Click the OK button to start the GUI. This may take some time. (If there is a licensing issue or an issue with the connection to the license server, this would typically only show up as an error now, not any earlier.)

Now that you have the GUI on the screen, you can now access the documentation. See below for more pointers on that. The most important parts of the GUI of COMSOL Multiphysics are the large work plane to draw the domain of the PDE and the pull-down menues Draw, Physics, Mesh, Solve, and Postprocessing. These menues correspond to the 5 basic steps necessary to solve a problem in COMSOL. The following list gives an example.

#### Step-by-Step Instructions how to Solve an Example Problem

For many of the following steps, buttons are available in the GUI that can be used as short-hand to access certain features from pull-down menues. It is generally easier for me to explain what to select from a pull-down menu than to describe the short-hand button, so that is why you will see me not reference the short-hand buttons very much, even when they exist for a particular function.

The following example is the Poisson equation -div(grad(u))=f with f(x,y) = (pi^2/2)*(cos(pi*x)*(cos(pi*y/2))^2+(cos(pi*x/2))^2*cos(pi*y)) in the domain (-1,1)x(-1,1) with homogeneous Dirichlet boundary condition (BC) u=0.

1. When the GUI starts, you are automatically in Draw mode, in which you draw the geometry of the domain of the PDE. If you are not in Draw mode for some reason, you can get to it by selecting Draw Mode from the Draw pull-down menu. Now, you have a number of basic geometry shapes on the left side of the work plane. To draw the domain (-1,1)x(-1,1), select from the pull-down menu Draw > Draw Objects > Rectangle/Square (Centered). Place the center at the origin (0,0) by pressing the left mouse button, then continue to hold the button and slide the mouse to the corner point (1,1), then release. You should have the domain (-1,1)x(-1,1) now. If you did not get the right shape (but do have a rectangle on the screen), you can double click on the geometry object to bring up a dialog window, in which you can enter the dimensions and placement of the object numerically. After the domain is drawn, it is convenient to re-center the work plane by selecting Options > Zoom > Zoom Extents.
2. Next, we will specify the coeffients in the PDE and the BC under the Physics pull-down menu. First, select Physics > Subdomain Settings. The top of the window shows the PDE being solved with the coefficient values visible on the right side of the window. Make sure to highlight the subdomain under the Subdomain tab on the left side. Only then will the input boxes become active on the right side. To solve the problem with the source term f given above, enter this expression in the field for f. This is the only coefficient that we needed to change, so click OK now. Second, select Physics > Boundary Settings. The top of the window shows the formulas of the boundary conditions available. To select all segments of the boundary, select Select by group in the lower left corner first, then highlight any of the boundaries on the left side of the window. Dirichlet BC should be preselected, and the coefficients h and r are already set for a homogeneous BC of u=0, so nothing needs to be changed. Click OK to finish.
3. Now mesh the domain by selecting Mesh > Initialize Mesh. The simple triangle is the short-hand button for this step. A mesh of the domain should now apper in the work plane. On my computer, this resulted in a mesh with 521 mesh points, 980 elements, and 2021 degrees of freedom (using the default quadrative Lagrange elements), as seen in the Mesh > Mesh Statistics.
4. Solve the problem by selecting Solve > Solve Problem or clicking on the short-hand button = (the equal sign).
5. After a brief solution process, a two-dimensional colored surface plot should appear. To see a three-dimensional colored surface plot of the solution, click on the second short-hand button on the left side of the work plane. The full capabilities for solution plots are available from the Postprocessing > Plot Parameters dialog window.
Afternote: If you solved exactly the problem discribed in the above steps, the true solution of the problem is (cos(pi*x/2))^2*(cos(pi*y/2))^2. To check that COMSOL obtained the correct solution, you could compute the integral over the domain of the square of the error between the FEM solution u and the formula for the true solution, that is, the integral of (u-(cos(pi*x/2))^2*(cos(pi*y/2))^2)^2. This is available under Postprocessing > Subdomain Integration. In the dialog window, select the subdomain, then enter the above formula for the error in the input field Expression, and click OK. The log field at the bottom of the COMSOL window should show the result, which I found to be 8.883395e-9 on my computer for example. Notice that this integral is the square of the L2-norm of the error, that is, the error norm is sqrt(8.883395e-9)=9.425176e-5.

To confirm that this is an appropriate result, we can refine the mesh once by secting Mesh > Refine Mesh, then solving again by selecting Solve > Solve Problem, and then computing the domain integral of the squre of the error again. This gave 1.340048e-10 on my computer and thus sqrt(1.340048e-10)=1.157604e-5 for the error norm. The ratio of the error on the coarse mesh divided by this error on the fine mesh is 9.425176e-5/1.157604e-5=8.141969. This means that by halving the mesh spacing in the uniform refinement process, we reduced the error by a factor of about 8, which agrees with the theoretical expectation of cubid convergence for quadratic Lagrange elements.

## Some Pointers to the Documentation of COMSOL Multiphysics

I will now detail how to get to the documentation of COMSOL Multiphysics. I am assuming that you have the main GUI on the screen as detailed above. Click on the Help menu on the top-right and select Help Desk (HTML/PDF), which should be the first item; alternatively you could click on the question mark button "?" on the far-right or even press the F1 function key on your keyboard. If you are not presently running a web browser, one will be started with the Documentation. If you had a browser already (possibly not visible), then a new tab in there is created and the documentation displayed. There are two columns of documentation, an HTML based version suitable for interactive searches and PDF files suitable for printing. The PDF files represent individual books for each subject, click on one desired PDF link at a time. The HTML version is all-encompassing with each of all manual chapters listed, and you can click on any of the HTML links to see the full documentation.

## Other Ways to Run COMSOL

The instructions for the sample problem above explained how to use COMSOL Multiphysics by running its graphical user interface (GUI). There are several additional options.
• Running the GUI together with COMSOL Script or MATLAB:
COMSOL Multiphysics provides interfaces between its GUI and the command-line driven environments of COMSOL Script and MATLAB. This is useful for importing and exporting data (e.g., the COMSOL fem struct that holds all information about a simulation) as well as for running scripts and functions (collectively called m-files). On the one hand, if you have already the GUI running, you can start COMSOL Script or MATLAB under the File menu of the GUI; to run either of these, you need to have a separate licence for that software. On the other hand, you can start COMSOL Script by itself first, by selecting it under Windows (and Mac) or by calling comsol script under Linux. You can then start the GUI from the File menu of the Script window.
• Running COMSOL in batch mode under Linux:
Notice that even running COMSOL Script entails opening a new window, including under Linux. Hence, you can only use this, if you are able to open a new window on your screen, which might not be the case if connecting remotely, e.g., from a Windows machine to a Linux machine.
Under the Linux operating system, there is the additional option for running COMSOL in batch mode, which avoids opening any windows (unless you explicitly program graphics commands in your script) and can thus be done within any text-based shell. To use this, you need to have a m-file script of the COMSOL Script commands that you wish to execute. The examples below discuss how to create such a file called driver_getfem.m. (Notice that this must be a script, not a function.) To run this script, enter at the Linux prompt:
```  comsol batch driver_getfem
```
Notice that the argument is the name of the script driver_getfem, which is the file name driver_getfem.m without the extension .m. COMSOL will now run and output to the shell windows whatever it would output in the Script or MATLAB window.
The above implies that you must remain logged in to the shell window, in which you are running COMSOL in batch mode. To truly run COMSOL in the background, you need to capture the screen output to a file. This can be done by re-directing the output. Assuming you are using the tcsh shell of Linux, adding >& driver_getfem.log to the above command will re-direct both stdout and stderr to the file driver_getfem.log. To additionally run COMSOL in the background and allow you to log out, add the ampersand & at the and and put nohup at the beginning of the command. So, the line
```  nohup comsol batch driver_getfem >& driver_getfem.log &
```
runs COMSOL in batch mode in the background and captures all screen output to the file driver_getfem.log. This is suitable for long runs, e.g., overnight of COMSOL.

1. At the COMSOL Conference 2008, I presented a paper with detailed convergence tests with Lagrange finite elements with degrees ranging from 1 to 5 for three smooth and non-smooth test problems.
• Preprint of the paper in PDF-form (6 pages)
• Slides of my presentation in PDF-form (15 slides)
2. At the COMSOL Conference 2007, I presented a paper with detailed instructions how to assess the solution quality of a given FEM solution by performing a convergence study on a sequence of progressively finer meshes. This technique uses the GUI of COMSOL Multiphysics in conjunction with COMSOL Script and might also be interesting just to get a feel for how the two can be used together. I am also posting the two m-files created used for the paper.
3. The following files are analogous to the ones above, but for an example of a problem that has a true solution known in closed form. It generalizes the discussion regarding convergence in the Afternote to the example described in the step-by-step instructions above. The driver function driver_getfem_mod.m is an improved version of the driver script driver_getfem.m above, with more features, such as accepting a variable for the getfem function, formatted output for the table data including additionally the number of elements in the mesh, and computation of true error if desired.
Here is a sample of the COMSOL Script session, in which the inputs to the function are set up first, then driver_getfem_mod is called twice, once with and once without the true solution supplied. This results first in a table of convergence data for the true error and then for the error against the reference solution. In the tables, the columns are the refinement level r = 0, 1, ..., 7, the number of elements N, the number of degrees of freedom D, the (true or estimated, resp.) error E, the ratio R of consecutive errors, and the (true or estimated, resp.) convergence order Q. Notice that the initial mesh with 4 elements created in getfem_cos_mod.m is the coarsest rotationally symmetric mesh possible, thus it takes a couple of refinement levels for the convergence order to reach the value of 2 expected for linear Lagrange elements.
```C>> h_getfem = 'getfem_cos_mod'
h_getfem =
getfem_cos
C>> nrefmax = 8
nrefmax =
6
C>> h_utrue = '(cos(pi*x/2))^2*(cos(pi*y/2))^2'
h_utrue =
(cos(pi*x/2))^2*(cos(pi*y/2))^2
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       4       5  2.9263e-001          NaN          NaN
1       8      13  2.0899e-001  1.4002e+000  4.8562e-001
2      16      41  7.3679e-002  2.8365e+000  1.5041e+000
3      32     145  2.0178e-002  3.6515e+000  1.8685e+000
4      64     545  5.2376e-003  3.8526e+000  1.9458e+000
5     128    2113  1.3264e-003  3.9486e+000  1.9813e+000
6     256    8321  3.3300e-004  3.9833e+000  1.9940e+000
7     512   33025  8.3358e-005  3.9948e+000  1.9981e+000
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax);
r       N       D            E            R            Q
0       4       5  2.7731e-001          NaN          NaN
1       8      13  1.5046e-001  1.8431e+000  8.8210e-001
2      16      41  5.3111e-002  2.8330e+000  1.5024e+000
3      32     145  1.4555e-002  3.6488e+000  1.8674e+000
4      64     545  3.7743e-003  3.8564e+000  1.9473e+000
5     128    2113  9.5361e-004  3.9579e+000  1.9847e+000
6     256    8321  2.3594e-004  4.0418e+000  2.0150e+000
7     512   33025  5.4898e-005  4.2977e+000  2.1036e+000
```
4. This example extends the previous one from linear Lagrange elements, that is, with shape functions of polynomial degree 1, to higher-order elements with shape functions of polynomial degrees 2, 3, 4, and 5. That is accomplished by modififying the file getfem_cos_mod.m from the previous example by successively changing Lag1 to Lag2, Lag3, Lag4, and Lag5 to create files getfem_cos_mod_2.m, ..., getfem_cos_mod_5.m, respectively. The true error is the computed by calling the same driver function driver_getfem_mod.m from the previous example with nrefmax = 5 and h_utrue = (cos(pi*x/2))^2*(cos(pi*y/2))^2 again. This gives the tables
```C>> h_getfem = 'getfem_cos_mod_2'
h_getfem =
getfem_cos_mod_2
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       4      13  1.7392e-001          NaN          NaN
1       8      41  3.1510e-002  5.5196e+000  2.4646e+000
2      16     145  3.8044e-003  8.2825e+000  3.0501e+000
3      32     545  5.1077e-004  7.4484e+000  2.8969e+000
4      64    2113  6.5472e-005  7.8013e+000  2.9637e+000
C>> h_getfem = 'getfem_cos_mod_3'
h_getfem =
getfem_cos_mod_3
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       4      25  7.1091e-002          NaN          NaN
1       8      85  6.5707e-003  1.0819e+001  3.4356e+000
2      16     313  6.4204e-004  1.0234e+001  3.3553e+000
3      32    1201  4.0756e-005  1.5753e+001  3.9776e+000
4      64    4705  2.5597e-006  1.5922e+001  3.9929e+000
C>> h_getfem = 'getfem_cos_mod_4'
h_getfem =
getfem_cos_mod_4
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       4      41  1.2803e-002          NaN          NaN
1       8     145  9.0950e-004  1.4077e+001  3.8153e+000
2      16     545  2.7186e-005  3.3454e+001  5.0641e+000
3      32    2113  8.9823e-007  3.0267e+001  4.9197e+000
4      64    8321  2.8497e-008  3.1520e+001  4.9782e+000
C>> h_getfem = 'getfem_cos_mod_5'
h_getfem =
getfem_cos_mod_5
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       4      61  8.1515e-003          NaN          NaN
1       8     221  9.2938e-005  8.7708e+001  6.4546e+000
2      16     841  2.5066e-006  3.7078e+001  5.2125e+000
3      32    3281  3.9463e-008  6.3517e+001  5.9891e+000
4      64   12961  6.1840e-010  6.3815e+001  5.9958e+000
```
These results show that the L2-norm of the error in the solution of the FEM with Lagrange elements with shape functions of polynomial degree p converges with order p+1 for this test smooth test problem. This is consistent with the theory for higher-order elements, assuming its boundary is smooth.

This demonstration used the fact that the true solution is known for this test problem. It is not readily clear how to estimate the error precisely enough, if the true solution is not known.

5. This example applies the techniques of the previous example to a problem, whose solution does not have the regularity to justify expecting the FEM to converge. This is true, even though the true solution is known in analytical form. The problem is the Poisson equation with a Dirac delta distribution as source term and homogeneous Dirichlet BC on the unit disk. This problem is discussed in the COMSOL documentation in the chapter Modeling Physics and Equations under the section on Specifying Point and Edge Settings as an example for using point constraints to implement a weak term. In two spatial dimensions, we can expect linear Lagrange elements to converge with order 1 (instead of the usual 2). We use the same driver function driver_getfem_mod.m from the previous example.
The following shows the setup of the true solution and the results in the COMSOL Script session, where the files getfem_del_disk_2.m, etc. are identical to getfem_del_disk.m with Lag1 changed to Lag2, etc. These results show that the convergence order does not improve beyond 1 for higher-order elements.
```C>> h_getfem = 'getfem_del_disk'
h_getfem =
getfem_del_disk
C>> nrefmax = 5
nrefmax =
5
C>> h_utrue = '-log(sqrt(x^2+y^2))/(2*pi)'
h_utrue =
-log(sqrt(x^2+y^2))/(2*pi)
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       8       9  4.2659e-002          NaN          NaN
1      16      25  2.1727e-002  1.9634e+000  9.7334e-001
2      32      81  1.0827e-002  2.0067e+000  1.0048e+000
3      64     289  5.4008e-003  2.0048e+000  1.0035e+000
4     128    1089  2.6990e-003  2.0010e+000  1.0008e+000
C>> h_getfem = 'getfem_del_disk_2'
h_getfem =
getfem_del_disk_2
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       8      25  2.6621e-002          NaN          NaN
1      16      81  1.1271e-002  2.3619e+000  1.2400e+000
2      32     289  5.6492e-003  1.9952e+000  9.9652e-001
3      64    1089  2.8240e-003  2.0004e+000  1.0003e+000
4     128    4225  1.4120e-003  2.0000e+000  1.0000e+000
C>> h_getfem = 'getfem_del_disk_3'
h_getfem =
getfem_del_disk_3
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       8      49  2.3227e-002          NaN          NaN
1      16     169  1.2267e-002  1.8934e+000  9.2098e-001
2      32     625  6.1349e-003  1.9996e+000  9.9972e-001
3      64    2401  3.0674e-003  2.0000e+000  1.0000e+000
4     128    9409  1.5337e-003  2.0000e+000  1.0000e+000
C>> h_getfem = 'getfem_del_disk_4'
h_getfem =
getfem_del_disk_4
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       8      81  9.1463e-003          NaN          NaN
1      16     289  3.3492e-003  2.7309e+000  1.4494e+000
2      32    1089  1.6752e-003  1.9993e+000  9.9950e-001
3      64    4225  8.3758e-004  2.0000e+000  1.0000e+000
4     128   16641  4.1879e-004  2.0000e+000  1.0000e+000
C>> h_getfem = 'getfem_del_disk_5'
h_getfem =
getfem_del_disk_5
C>> [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
r       N       D            E            R            Q
0       8     121  1.0465e-002          NaN          NaN
1      16     441  5.2186e-003  2.0053e+000  1.0038e+000
2      32    1681  2.6100e-003  1.9995e+000  9.9962e-001
3      64    6561  1.3050e-003  2.0000e+000  1.0000e+000
4     128   25921  6.5250e-004  2.0000e+000  1.0000e+000
```

As a common courtesy, if any information on this webpage is useful to you and you use it, you should give credit to me and refer your audience to this site (refer to it as the COMSOL area of www.math.umbc.edu/~gobbert).