From 76dbed2c3df94696df218718da3bf8dbeb276546 Mon Sep 17 00:00:00 2001 From: Jaime Arriaga Garcia Date: Tue, 12 Nov 2024 14:46:28 +0100 Subject: [PATCH] some corrections --- book/fvm/advection_1.ipynb | 46 ++++++++++++++++-------------- book/fvm/getting_started_new.ipynb | 9 +++++- 2 files changed, 32 insertions(+), 23 deletions(-) diff --git a/book/fvm/advection_1.ipynb b/book/fvm/advection_1.ipynb index bf06c02..a25eb62 100644 --- a/book/fvm/advection_1.ipynb +++ b/book/fvm/advection_1.ipynb @@ -106,7 +106,7 @@ "\n", "The integral equations may now be discretized over each finite volumes with the \"unknown\" at the center, $\\phi_i$, and with dimensions $\\Delta x$, $\\Delta y$ and $\\Delta z$. .\n", "\n", - "For the first term: velocity $\\phi$ is constant over the domain, the integral is applied over the finite volume $i$ and a forward Euler scheme is applied with time-step $\\Delta t$:\n", + "For the first term: velocity $c$ is constant over the domain, the integral is applied over the finite volume $i$ and a forward Euler scheme is applied with time-step $\\Delta t$:\n", "\n", "$$\n", "\\int_{V}\\frac{\\partial \\phi}{\\partial t}dV\n", @@ -147,10 +147,10 @@ "\n", "$$\n", "\\int_{S} c \\, \\phi \\cdot \\mathbf{n} \\, dS\n", - "= -c \\, \\phi_{w} \\, \\Delta y \\Delta z + \\phi_{e} \\, \\Delta y \\Delta z = c\\left(\\phi_{e}-\\phi_{w}\\right)\\Delta y \\Delta z\n", + "= -c \\, \\phi_{w} \\, \\Delta y \\Delta z + c \\phi_{e} \\, \\Delta y \\Delta z = c\\left(\\phi_{e}-\\phi_{w}\\right)\\Delta y \\Delta z\n", "$$\n", "\n", - "where $\\phi_e,w$ are evaluated at the midpoint of the east and west surfaces. As those values are not yet known at the surface (they are known at the volume centroid), an approximation is needed. A linear interpolation can be used, thus only two points per surface value are needed: \n", + "where $\\phi_{e,w}$ are evaluated at the midpoint of the east and west surfaces. As those values are not yet known at the surface (they are known at the volume centroid), an approximation is needed. A linear interpolation can be used, thus only two points per surface value are needed: \n", " \n", "$$\n", "\\phi_{e} = \\frac{\\phi_{i+1} + \\phi_i}{2} \\qquad\n", @@ -319,10 +319,12 @@ "\n", "The initialization of the problem is delicate as only after the second time step the algebraic equation is fully satisfied without additional requirements. \n", "\n", - "The coefficient matrix $A$ only has non-zero elements along the diagonal due to the regular geometry incorporated in the 1D discretization, called a _banded_ matrix. In the derivation of the system of equations, the boundary conditions of end volumes were not presented. Actually, the advection equation has only one boundary condition due to the first spatial derivative! It is defined at the boundary, from where $\\phi$ is being advected. The formulation of the coefficients in the matrix implies that the values of the $0^{\\textrm{th}}$ and $7\\textrm{th}$ volume ($\\phi_{0}$ and $\\phi_7$) were taken as 0, the latter is an artifice to vanish the quantity $\\phi$. Rigorously speaking, it should not be a boundary condition.\n", + "The coefficient matrix $A$ only has non-zero elements along the upper and lower diagonal due to the regular geometry incorporated in the 1D discretization. In the derivation of the system of equations, the boundary conditions of end volumes were not presented. Actually, the advection equation has only one boundary condition due to the first spatial derivative! It is defined at the boundary from where $\\phi$ is being advected. The formulation of the coefficients in the matrix implies that the values of the $0^{\\textrm{th}}$ and $7\\textrm{th}$ volume ($\\phi_{0}$ and $\\phi_7$) were taken as 0, the latter is an artifice to vanish the quantity $\\phi$. Rigorously speaking, it should not be a boundary condition.\n", "\n", "### Solution Techniques\n", "\n", + "This system can be solved without the need for building the matrix as each value of the next time step can be computed from known values of previous steps. So, it can be solved with a for loop as well! \n", + "\n", "Note that the matrix in the above example contains a lot of zeros: 10 out of 36 elements are non-zero (28%). This is called a _sparse_ matrix, and the proportion of non-zero entries grows with the size of the problem. For example, $N=1000$ finite volumes requires $2N-2$ non-zero entries (1998), which is only 0.1% of the matrix! Since matrix calculations are computationally expensive, we would never carry out the matrix calculations as represented here. Instead, various vectorization approaches are implemented in all numerical analysis software packages. Typically, a mapping of volumes and coordinates are stored in a vector format that retains the relationship between volumes and the interpolation points (e.g., neighbor volumes). Then calculations per volume are carried out using this vectorized information. This can be easily implemented in any programming language, for example, with a for loop over each volume. However, note that as geometry of the domain of interest and the discretization scheme changes (e.g., non-square volumes!) more advanced vectorization techniques are required." ] }, @@ -412,6 +414,11 @@ "\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -433,33 +440,28 @@ " \\phi_i^{n+1} = \\phi_i^n + \\Delta t D \\frac{\\phi_{i-1}^n - 2\\phi_i^n + \\phi_{i+1}^n}{\\Delta x^2}\n", " $$\n", "\n", - " Remember the stability criteria of Forward Euler $|1- \\Delta t \\alpha| <1$\n", + " This can be rewritten as \n", + " $$\n", + " \\phi_i^{n+1} = (1-2*\\Delta t D/\\Delta x^2)\\phi_i^n + \\Delta t D \\frac{\\phi_{i-1}^n + \\phi_{i+1}^n}{\\Delta x^2}\n", + " $$\n", "\n", - " For the diffusion equation we get the following:\n", + " Therefore, the amplification factor should be lower than 1:\n", "\n", " $$\n", - " -1 < 1 - \\frac{\\Delta t D}{\\Delta x^2} <0\n", + " -1 < 1-2D \\frac{\\Delta t}{\\Delta x^2} < 1\n", " $$\n", " \n", " $$\n", - " -2 < - \\frac{\\Delta t D}{\\Delta x^2} <0\n", - " $$\n", - "\n", - " $$\n", - " \\Delta t< \\frac{ 2 \\Delta x^2}{D}\n", + " D \\frac{\\Delta t}{\\Delta x^2} <= 1/2\n", " $$\n", - "\n", - "\n", - " - **Backward Euler (Implicit)**:\n", - "\n", - " $$\n", - " \\phi_i^{n+1} = \\phi_i^n + \\Delta t D \\frac{\\phi_{i-1}^{n+1} - 2\\phi_i^{n+1} + \\phi_{i+1}^{n+1}}{\\Delta x^2}\n", - " $$\n", - "\n", - "\n", - "Backward Euler is unconditionally stable.\n" + "\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, { "cell_type": "code", "execution_count": null, diff --git a/book/fvm/getting_started_new.ipynb b/book/fvm/getting_started_new.ipynb index 3e8c9ac..06dc7ce 100644 --- a/book/fvm/getting_started_new.ipynb +++ b/book/fvm/getting_started_new.ipynb @@ -192,9 +192,11 @@ "$$\n", "\\frac{\\partial \\phi}{\\partial t}\n", "+ \\mathbf{u} \\cdot \\nabla \\phi\n", - "= D \\nabla^2 \\phi +F_{\\phi}\n", + "= D^* \\nabla^2 \\phi +F^*_{\\phi}\n", "$$\n", "\n", + "where $D^*$ and $F^*$ are $D$ and $F$ divided by $\\rho$. From now on, the asterisk will be dropped. \n", + "\n", "From left to right, these three terms are explained as follows, relative to the quantity $\\phi$:\n", "1. Transient effect: rate of change of quantity $\\phi$\n", "2. Convection: transport due to velocity field $\\mathbf{u}$\n", @@ -216,6 +218,11 @@ "\n" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, { "cell_type": "markdown", "metadata": {},