This assignment considers an interesting 2D shape: the "U" of MUDE! We will use the diffusion equation to compute the distribution of temperature in the U subject to specific boundary conditions and initial conditions.
This assignment contains three parts:
1. expressing the diffusion equation for triangle volumes and formulating algebraic equations with the finite volume method,
2. implementing and solving this method in an unstructured orthogonal mesh,
3. analyzing potential downsides to the discretized approach used if applied to non-orthogonal meshes, and considering how they might be corrected.
Remember, even though the problem is 2D, the numerical scheme will be fromulated by treating _volumes_ (as opposed to dimensions $x$ and $y$ directly).
%% Cell type:markdown id:4fa20681 tags:
## Part 1: Using Just Your Hands
The diffusion equation expressed in its reduced form is:
$$
\frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \phi
$$
Note that for triangle shape volumes, the fluxes are not necessarily directed in the $x$ direction. Rather, in the normal direction to the surfaces of each volume, as the flux direction is given by the $\phi$ gradients and the information propagates in all directions.
Over the next tasks, you will write step by step its derivation.
Approximate the derivatives! Use central differences in space and Forward Euler in time. The subindices can be defined by arbitrarily naming the volume of interest and the adjacent volumes.
The equation you derived above describes volumes with <b>three</b> interior sides. For a volume with <b>one</b> side being an exterior one, modify the equation above to implement the following Neumann condition (which replaces the discretized gradient term):
+ D\frac{\Delta t \Delta L}{\Delta A} \left( \frac{\phi^n_{v1}-\phi^n_i}{\Delta d_c} \right)
+ D\frac{\Delta t \Delta L}{\Delta A} \left( \frac{\phi^n_{v2}-\phi^n_i}{\Delta d_c} \right)
+ D\frac{\Delta t \Delta L}{\Delta A} \cdot 10
$$
</p>
</div>
%% Cell type:markdown id:fb25fa6b tags:
## Part 2: Implementation!
Below, the coordinates of the triangle vertices that cover the domain of interest are defined. The boundary conditions are also specified and incorporated in the resulting object `mesh` once the class `Mesh` is instantiated. The class `Mesh` is defined in `utilities.py` and has an identical structure and similar functions as provided in your PA for this week, which define key characteristics of the mesh and volumes, as well as providing useful plotting and solving methods.
Note: the `Mesh` class and it's usage is illustrated extensively in the companion notebook `mesh_tips.ipynb`.
Execute the cell below and study the resulting geometry until you are comfortable recognizing the problem that will be solved. In particular, check that you can identify which boundary conditions are applied, and where, before moving on.
If you are not sure how to interpret the code or figures, refer to <code>mesh_tips.ipynb</code>.
</p>
</div>
%% Cell type:code id:e10227a9 tags:
``` python
length = 10
coordinates = [[ 0., 0. ],
[ 10., 0. ],
[ 5., -8.660],
[ 15., -8.660],
[ 10., -17.320],
[ 20., -17.320],
[ 25., -8.660],
[ 30., -17.320],
[ 35., -8.660],
[ 30., 0. ],
[ 40., 0. ]]
coordinates = np.array(coordinates)
boundaries = [[['Neumann', +1], [0, 1]],
[['Neumann', 0], [1, 3, 6, 9]],
[['Neumann', -1], [9, 10]],
[['Neumann', 0], [10, 7, 4, 0]]]
mesh = Mesh(coordinates, length, boundaries)
mesh.plot_triangles();
mesh.plot_boundaries();
```
%% Output
%% Cell type:markdown id:942e3065 tags:
Before we continue to finding the solution, let's explore the relationship of today's mesh with the bars and kapsalon shops we considered in the PA for this week.
Write down the similarities between the housing-bar-kapsalon problem (PA 2.1) and the FVM representation of the diffusion equation here. Specifically, the locations and quantities associated with each of the three things.
- **Houses**: these are the vertices of the volumes; we don't actually compute anything at these locations, they are only used to define the triangle/mesh geometry.
- **Bars**: these are the centroids of the triangles and are the locations where the unknowns $\phi(x, y, t)$ are computed. Note however, we compute something more like $\phi(i,t)$ where $i$ refers to the triangle number.
- **Kapsalon shops**: not only are these places where we can get tasty food, but also they are the centers of the side of each triangle (but only when the triangles are equal size and equilateral: orthogonal). This is the midpoint of the surface integral, an important observation that you probably would not have recognized in this task, but is considered in the third part of this notebook.
Suppose **just for this single task** that the FDM would be used, sketch the potential grid points and draw the boundaries that would represent the same domain in $x$ and $y$ as defined by <code>coordinates</code> in the previous task.
This is to illustrate one of the contrasts between FDM and FVM.
Sketch several sets of grid points to illustrate the shape of the U; you will need several. Depending on how "wide" the U is (number of grid points), the geometry has trouble capturing the angles. Especially the boundary is not a smooth line, as "steps" are introduced to capture the boundary.
It is also possible to reshape the U into a rectangular shape; this would be much easier to solve with a simple FDM scheme, however, the problem is very different and may not be representative of the situation in reality that must be modelled.
Main takeaway is to see that FDM has challenges modelling non-rectangular geometry.
</p>
</div>
%% Cell type:markdown id:3c7ee41a tags:
Now we will **continue with code implementation** by adding some information to the object that defines our problem, `mesh`. Remember to refer to `mesh_tips.ipynb` if you are not sure how to use it!
First, use the method <code>set_initial_temperature()</code> to define the initial conditions for the volumes. The initial temperature should be 20 degrees everywhere except at the volume in the middle, where the temperature is 40 degrees.
</p>
</div>
%% Cell type:code id:6ae17c54 tags:
``` python
# mesh.plot_triangles(YOUR_CODE_HERE) # useful for identifying the triangle id
Using the <code>solve</code> method, solve the problem for conditions where $t_{final}=20$, $N_t=100$ and $D=50$.
Use the resulting plot to see if you reach the solution you expect. Remember to consider the boundary conditions that were defined for you as well, not only the initial conditions.
In the previous task, you should have realized that the solution is _not_ what we expect.
It turns out we screwed up! In the massive file <code>utilities.py</code>, there is </b>one</b> line of code that needs to be fixed in order to solve this problem. Can you find it and fix it?
The first group to find the solution will win a prize!
<em>Hint: it's not in the plotting method, and is related to something you derived in the tasks above.</em>
<b>Note:</b> you should expect to spend some time reading the code (at least one very specific part of the code); not only will this help you fix the problem, but it will be very useful for answering the first few questions in the Report.
The process of reading the code and solving this problem is directly related to the first few questions of the Report.
The first thing to recognize is that the "solution" using the incorrect code looks exactly like the initial conditions. So this means that we are probably calculating something incorrectly in the solution algorithm, which is contained in the method <code>solve</code>.
It is challenging to read through the code in this method because there are a lot of loops and conditional statements. However, one should be able to recognize the time and space discretization steps, as well as the conditional statements which identify neighbor triangles to perform the surface flux calculations and apply the boundary conditions, when necessary. Note in particular that the "space" integration actually is a loop over all triangles (the volumes), and for each triangle there is a loop over each face.
Note that Python method <code>enumerate</code> returns the index of a list/array, as well as the item at that index; this is useful to more easily define values in the loop.
Once the structure of the code is recognized, it should be straightforward to check that the algebraic equations defined above are implemented correctly (note variables like <code>unknowns</code>, <code>phi</code>, <code>constant</code> and <code>flux</code>).
At the end of this method, we can see how the Euler scheme is applied:
After close inspection, it should be obvious that this simply takes the solution from the previous time step for the triangle with given id, $\phi_{i}^{n}$ and assigns the value to the next time step, $\phi_{i}^{n+1}$. Because this happens for all time steps and all triangles, this results in a "solution" that simply repeats the initial conditions at every time step.
To fix the issue, we need to implement the algebraic expression
The solution for the refined mesh did not work---this time it is not a problem with the code, but another issue. See if you can make some adjustments and fix the solution!
Note that you will be asked about stability in the Report, so you might as well calculate this and record the values now.
<em>Tip: even though the mesh has smaller volumes (triangles), the solution should look similar to the previous (unrefined) mesh.</em>
First, note that the unstable solution is recognized by the large values of the solution indicated by the color bar in the figure: order 10 to the 87th power! Clearly the solution is on its way to infinity (and beyond?!).
Instead of solving with these parameters:
<pre>
<code>
mesh.solve(20, 100, 50)
</code>
</pre>
Increasing the number of time steps by a factor of 10 works:
<pre>
<code>
mesh.solve(20, 100, 50)
</code>
</pre>
This reduces the time step size $\Delta t$ by a factor of 10, and the reason it results in a stable solution can be explained using the stability criteria...
</p>
</div>
%% Cell type:code id:a3c2b83a tags:
``` python
side_length= 5.0
t_final=20
D=50
Nt_vec=100 * np.arange(1, 11)
Nt_vec=np.arange(200, 650, 50, dtype=int)
Stability = []
max_temp = []
for Nt in Nt_vec:
mesh.solve(20, Nt=Nt, D=D)
_,stability,temp=mesh.solve(20, Nt=Nt, D=D)
Stability.append(stability)
max_temp.append(temp.max())
print("==================================")
print()
plt.figure(figsize=(10, 5))
plt.plot(Nt_vec, Stability, 'o-')
for x, y, temp in zip(Nt_vec, Stability, max_temp):
## Part 3: Evaluating the implementation for non-equilateral triangles
Computations in meshes with non-equilateral triangles have added error sources that would need to be corrected to have an accurate solution. In this section you will analyse and reflect on the potential downsides of your implementation for non-equilateral triangle volumes by looking only at the fluxes between two volumes. The file `utilities.py` and rest of the code in this notebook is completely irrelevant to this Part (except the `numpy` and `matplotlib` import).
The vertex coordinates of the first triangle to be analyzed are: [0,0] , [0,1] , [1,0]. The second triangle share the first two vertices and have coordinates [0,0] , [0,1] , [-1,0].
Plot the centroid of each triangle, its connecting line and the middle distance. For the scheme you implemented in the previous task, write down what this implies for the **numerical integration approximation of the flux at the surface**.
Specifically, consider the location where the flux is calculated, as implied by the numerical scheme, as well as the integration that is implied by the analytic surface integral over the surface of the shared triangle side.
When you approximate the integral with a numerical one, it is assumed that one value of the flux is representative of the fluxes at the entire surface (recall numerical integration of Q1). For equilateral triangles, using the centroids to approximate the gradient implies that we are using a midpoint rule that is second order accurate. However, as the figure shows, for this situation the evaluation is _not_ done at the middle point; rather, it is representative of another location on the surface; thus the accuracy is not second order accurate.
Now lets analyze a different situation and change the coordinates of the second triangle to [0,0] , [0,1] , [-1,1].
Plot the centroid of each triangle, its connecting line and the middle distance. For the numerical scheme you implemented, write the implication for the **numerical derivative approximation of the flux at the surface**.
Specifically, consider the _direction_ in which the flux is calculated, as implied by the numerical scheme, as well as the quantity demanded by the analytic expression (i.e., the surface integral).
Now the evaluation of the flux is done at the midpoint of the surface (in contrast to the previous task) but its direction evaluation is not normal to the surface!
Thus, the gradient and flux are computed incorrectly if we applied the scheme implemented in our code.
Now lets analyze one last situation and change the coordinates of the second triangle to [0,0] , [0,1] , [-1.5,1].
Plot the centroid of each triangle, its connecting line and the middle distance. Write down what this implies for the **numerical derivative approximation of the flux normal to the surface**.
Now the evaluation of the flux (as implemented in task 2) is not done at the midpoint of the surface, not even at the surface, and its direction is not normal to the surface!
Thus, there are two sources of error present, when it comes to the numerical scheme in the code of this assignment, applied to the volumes illustrated here.
