Numerical Solution of PDE: Comsol Multiphysics: Example 1, The Poisson Equation On An Ellipse
Numerical Solution of PDE: Comsol Multiphysics: Example 1, The Poisson Equation On An Ellipse
Numerical Solution of PDE: Comsol Multiphysics: Example 1, The Poisson Equation On An Ellipse
ApplicationModes>ComsolMultiphysics>PDEModes>ClassicalPDEs>Poissons Equation. We see that the element type is Lagrange P2-triangles. The solution is approximated
with second degree polynomials, one over each triangle, a globally C0 function. The elements can be changed later, if you wish. But the name of the solution (here u) must be changed here.
The next screen shows the work plane (2D ) and one usually works from left to right in the menu bar: 1. Create geometry (Draw) 2. Define the PDE (Physics) 3. Set boundary conditions (Physics) 4. Make the mesh (Mesh) 5. Solve (Solve) 6. Plot (Postprocessing) Under File we find operations like open, save, save as, etc. Options gives facilities like defining constants and expressions. It is practical to have names for important model parameters like heat transfer coefficients, specific heat, etc.
Geometry The geometry is built in a drawing program, from elementary bodies (2D Rectangles and ellipses), or surfaces bounded by Beziercurves. The bodies are put together by set operations (union, intersection, difference). We make an ellipse, centered at the origin and half axes 1 and 0.2. Note that the objects are snapped to the grid, change the grid under Options>GridSettings if needed. Our domain is only the single ellipse so we move on to Define PDE
Physics>SubdomainSettings
(there is
boundary too, in a minute) The PDE is in the Equation pane, ( c u ) f and we need to define the coefficient c and the source term f. One can give functions of x,y, u, ux, uy and, if time-dependent, t. There is only one subdomain, the ellipse, so click it and accept a c = f = 1 by OK.
Physics>BoundarySettings
Each ellipse gives four boundaries. Choose them all by ctrl-a; They can be selected also by clicking (shift-clicking, ctrl-clicking) on the work pane. The default boundary condition is apparently Dirichlet: h u = r with h = 1 and r = 0 and we change to r = sin(10x) for some variation, OK.
Element mesh The triangle icon gives a mesh. It looks good, OK. The mesh can be controlled in some detail by parameters under Mesh>MeshParameters .
Solution
The generated equations are solved when you click the = icon. Scroll down to Solve>SolverParameters and check the possibilities offered, Our choice of c, f and boundary conditions makes the problem linear so the default selection stationary linear is correct, OK. Details that can be very important for large (and non-linear) problems, like choice of solver for linear systems, are controllable here. It says Unsymmetric about the matrix whereas it is in fact symmetric but never mind, OK, and solve by =.
After a fraction of a second, the system with 613 unknowns has been constructed and solved and painted as colored isosurfaces on the work pane, and the logline (bottom of screen) says
Number of degrees of freedom solved for: 613 Solution time: 0.093 s
Convergence study To see how the solution converges with refinement of the mesh, we compute u in the center. Choose
Postprocessing>DataDisplay>Subdomain
where we give the coordinates for the center (0,0) and see
Value: 0.020201, Expression: u, Position: (0,0)
in the log-line. Make a regular refinement (each triangle split into four) by the four triangle icon, solve again, and print u(center) again:
Number of degrees of freedom solved for: 2345 Solution time: 0.188 s Value: 0.019069, Expression: u, Position: (0,0)
and again, and again , ... try up to, say, a million degrees of freedom when the solver runs out of memory. Convergence:
u(0,0) diff. 1201 69 -1132 231 162 234 3 231 -3
The order of convergence should be 2 for regular solutions in energy norm and 3 in max-norm. The exponent is hard to see from the data, but the convergence is rapid. The
convergence is not as regular as you have seen for e,g, the trapezoidal rule for quadrature, because the mesh is so irregular. Now try piecewise linear basis functions instead. Lagrange P1 elements are chosen under Physics>SubdomainSettings>Element Repeat the experiment: Slower convergence.
#dof 613 2345 9169 36257 144193 time 0.094 0.156 0.547 2.172 11.204 u(0.0) 0.019103 0.019115 0.019205 0.019224 0.019230 diff. 12 90 19 6