User:Tclamb/FEniCS

 equation numbering. backward and forward reference. Egm6322.s12 (talk) 15:10, 20 July 2012 (UTC)

=Newton's Law of Cooling=

Problem Statement


Consider Newton's Law of Cooling. The PDE for the heat problem is given by:

Problem Statement






Problem Statement References




Obtaining the Weak Variational Form

 * $$\nabla \cdot \left( \mathbf \kappa \nabla u \right) + f(\mathbf x,t) = \rho c \frac{\partial u}{\partial t}$$

Substitute in given values:
 * $$\nabla \cdot \left( \mathbf I \nabla u \right) + 0 = \rho c \, 0$$
 * $$\nabla^2 u = 0$$

Multiply by test function:
 * $$\left( \nabla^2 u \right) v = 0$$

Integrate over domain:
 * $$\int_\Omega \left( \nabla^2 u \right) v \, d\Omega = 0$$

Apply Green's first identity:
 * $$\oint_\Gamma \left( \nabla u \cdot \hat \mathbf n \right) v \, d\Gamma - \int_\Omega \nabla u \cdot \nabla v \, d\Omega = 0$$

Split line integral and substitute in Neumann boundary conditions:
 * $$\left( \int_{\Gamma_h} -h v \, d\Gamma_h + \int_{\Gamma_H} -H \left( u - u_\infty \right) v \, d\Gamma_H + \int_{\Gamma_g} 0 v \, d\Gamma_g \right) - \int_\Omega \nabla u \cdot \nabla v \, d\Omega = 0$$

Separate into generic form, a(u,v) = L(v): $$\int_\Omega \nabla u \cdot \nabla v \, d\Omega + \int_{\Gamma_H} H u v \, d\Gamma_H = \int_{\Gamma_H} H u_\infty v \, d\Gamma_H - \int_{\Gamma_h} h v \, d\Gamma_h$$

Plots of Solution
My Solution:
 * [[image:NewtCooling.pdf|480px]]

This solution is in agreement with ccook:

While viewing the plot, the data can be saved as a "vtk" file by pressing "v". This vtk file can be opened in ParaView for further analysis. Below, a ParaView plot with isolines on every multiple of 0.25.


 * [[File:Nloc biunit.png]]]

ParaView can also be used to determine the value at a given point.

=Poisson's Equation with Known Analytic Solution=

Problem Statement

 * $$\Omega = \square$$
 * $$u \left(\mathbf x\right) = \sin \left(2 \pi x\right) \sin\left(2 \pi y\right)$$

Finding the Poisson Equation and Dirichlet Boundary Conditions
From the given solution, we take the Laplacian to find the following Poisson Equation:
 * $$- \nabla^2 u = f \quad \text{where } f = 8 \pi^2 \sin \left(2 \pi x\right) \sin\left(2 \pi y\right) = 8 \pi^2 u$$

Likewise, we can see that on the unit square, the given solution will satisfy the following boundary condition:
 * $$u = 0 \quad \forall \, x \in \partial\Omega$$

Obtaining the Weak Variational Form
Multiply by test function:
 * $$-\left( \nabla^2 u \right) v = f v$$

Integrate over domain:
 * $$-\int_\Omega \left( \nabla^2 u \right) v \, d\Omega = \int_\Omega f v \, d\Omega$$

Apply Green's first identity:
 * $$\int_\Omega \nabla u \cdot \nabla v \, d\Omega - \oint_\Gamma \left( \nabla u \cdot \hat \mathbf n \right) v \, d\Gamma = \int_\Omega f v \, d\Omega$$

Simplify:
 * $$\int_\Omega \nabla u \cdot \nabla v \, d\Omega = \int_\Omega f v \, d\Omega$$

Computed Solution

 * [[File:Poisson_Computed.png]]

Exact Solution

 * [[File:Poisson_Exact.png]]

=Steady, Incompressible Navier-Stokes=

Derivation
Starting from Navier-Stokes:
 * $$\rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot\boldsymbol{\mathsf{T}} + \mathbf{f}$$

By assuming a homogeneous, incompressible, Newtonian fluid, we can simplify the stress tensor, resulting in the following equation:
 * $$\rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}$$

By assuming steady flow, we can eliminate the time derivative of u, resulting in the following equation:
 * $$\rho \left(\mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}$$

Assuming no body forces, we simplify to this form:
 * $$\rho \left(\mathbf{u} \cdot \nabla \mathbf{u} \right) - \mu \nabla^2 \mathbf{u} + \nabla p = 0$$

By dividing through by $$\rho$$, we can eliminate the first term's coefficient and absorb the density into the pressure (effectively scaling the pressure):
 * $$\mathbf{u} \cdot \nabla \mathbf{u} - \frac\mu\rho \nabla^2 \mathbf{u} + \nabla p = 0$$

We now take the dot product of a vector valued test function for velocity, v, and integrate over the domain:
 * $$\int_\Omega \left( \mathbf{u} \cdot \nabla \mathbf{u} - \frac\mu\rho \nabla^2 \mathbf{u} + \nabla p \right) \cdot \mathbf{v} \; dx = 0$$

We can take advantage of the fact that the test function, $$\mathbf{v}$$, evaluates to 0 along the boundary to transform two of the terms into more manageable forms.

First, we manipulate the pressure term into a more useful form involving the divergence of the test function:
 * $$\int_\Omega \nabla p \cdot \mathbf{v} \; dx = \int_\Omega \frac{\partial p}{\partial x_i} v_i \; dx = \left. p v_i \right|_{\partial\Omega} - \int_\Omega p \frac{\partial v_i}{\partial x_i} \; dx = 0 - \int_\Omega p \left(\nabla \cdot \mathbf{v}\right) \; dx$$

Observing that in Cartesian coordinates, the following identity holds:
 * $$\left(\nabla^2 \mathbf{u}\right)_i = \nabla^2 u_i$$

We then manipulate the vector Laplacian term into one involving only first derivatives:
 * $$\int_\Omega \nabla^2 \mathbf{u} \cdot \mathbf{v} \; dx = \int_\Omega \frac{\partial^2 u_i}{\partial x_j^2} v_i \; dx =

\left. \frac{\partial u_i}{\partial x_j} v_i \right|_{\partial\Omega} - \int_\Omega \frac{\partial u_i}{\partial x_j} \frac{\partial v_i}{\partial x_j} \; dx = 0 - \int_\Omega \left(\nabla \mathbf{u}\right) : \left(\nabla \mathbf{v}\right) \; dx$$

Now, we can obtain the following weak form of the solution:
 * $$\int_\Omega \left( \mathbf{u} \cdot \nabla \mathbf{u} \right) \cdot \mathbf{v} + \frac\mu\rho \left(\nabla \mathbf{u}\right) : \left(\nabla \mathbf{v}\right) - p \left( \nabla \cdot \mathbf{v} \right) \; dx = 0$$

Or equivalently:
 * $$\int_\Omega \mathbf{u}\cdot\text{grad}(\mathbf{u})\cdot\mathbf{v} + \frac\mu\rho \, \text{grad}(\mathbf{u}) : \text{grad}(\mathbf{v}) - p \, \text{div}(\mathbf{v}) \; dx = 0$$

Example: 3D Flow in Elbow Pipe
 describe (1) how you found out about gmesh as an open-source mesh generator, (2) how you used it to create a mesh, (3) how you incorporate the mesh data into FENiCS (formatting, etc.) Egm6322.s12 (talk) 14:41, 20 July 2012 (UTC)

 explain the different degree of spatial interpolation for the velocity field and for the pressure field. why do you need to use 2nd order for the velocity field and first order for the pressure field ? the strain rate is proportional to the spatial gradient of the velocity field, and thus if the velocity field is interpolated with 2nd order polynomials, then the strain would be interpolated with 1st order polynomials. so both the strain rate and the pressure would have the same 1st order interpolation. but why do you need to have the same interpolation order for both the strain rate and the pressure ? document what we discussed today. Egm6322.s12 (talk) 15:10, 20 July 2012 (UTC)

Output

 * Ns.split.re4k.u.small.png
 * Ns.split.re4k.p.small.png

=Quantum Harmonic Oscillator=

Problem Statement
The Schrödinger equation for a wavefunction $$\Psi(x,t)$$ is
 * $$-\frac{\hbar^2}{2m}\nabla^2\Psi(\vec{x},t) + V(\vec{x})\Psi(\vec{x},t) = i \hbar \frac{\partial \Psi(\vec{x},t)}{\partial t}$$

Separation of variables can be performed, yielding an equation for the position-dependence of the wavefunction called the time-independent Schrödinger equation:
 * $$-\frac{\hbar^2}{2m}\nabla^2\psi(\vec{x}) + V(\vec{x})\psi(x) = E \psi(\vec{x})$$

Identifying the left-hand side as the Hamiltonian operator, this is an eigenvalue equation for the Hamiltonian with eigenfunctions $$\psi_n(x)$$ and eigenvalues $$E_n$$.

The potential for a harmonic oscillator potential is
 * $$V(x) = \frac12k\vec{x}^2$$

The natural frequency of an oscillator is $$\omega = \sqrt{\tfrac{k}{m}}$$, allowing the potential to be rewritten as
 * $$V(x) = \frac12 m \omega^2 \vec{x}^2$$

Therefore, the time-independent Schrödinger equation within a harmonic oscillator potential is
 * $$-\frac{\hbar^2}{2m}\nabla^2\psi(\vec{x}) + \frac12 m \omega^2 \vec{x}^2 \psi(\vec{x}) = E \psi(\vec{x})$$

This can be rendered unitless by the transformation $$\vec{x} = \sqrt{\tfrac{\hbar}{m\omega}} \vec{\xi}$$:
 * $$-\frac{\hbar \omega}{2}\nabla^2\psi(\vec{\xi}) + \frac12 \hbar \omega \vec{\xi}^2 \psi(\vec{\xi}) = E \psi(\vec{\xi})$$

Defining the unitless constant $$\epsilon = \frac{2 E}{\hbar \omega}$$, the final form of the unitless equation is
 * $$\nabla^2\psi(\vec{\xi}) + (\epsilon - \vec{\xi}^2) \psi(\vec{\xi}) = 0$$

In one dimension, this equation reduces to
 * $$\frac{\partial^2 \psi(\xi)}{\partial \xi^2} + (\epsilon - \xi^2) \psi(\xi) = 0$$

Obtaining the Weak Variational Form
For ease of typing, rename the variables in the one-dimensional equation above to obtain the following equation:
 * $$\frac{\partial^2 u}{\partial x^2} + (k - x^2) u = 0$$

Now, multiply by the test function:
 * $$\frac{\partial^2 u}{\partial x^2}v + (k - x^2) uv = 0$$

Integrate over the domain:
 * $$\int_\Omega \frac{\partial^2 u}{\partial x^2}v + (k - x^2) uv\,dx = 0$$

Integrate by parts, and multiply by negative one. Boundary terms are discarded as the wavefunction drops off exponentially.
 * $$\int_\Omega \frac{\partial u}{\partial x}\frac{\partial v}{\partial x} + x^2 uv\,dx = \int_\Omega k u v\,dx $$

Animations of Solutions
The squares of the computed eigenfunctions generated with FEniCS are animated on the left. Corresponding analytical plots generated in Mathematica are animated on the right.
 * [[Image:U2.gif]][[Image:U2_anal.gif]]