<p><b>mpetersen@lanl.gov</b> 2013-03-01 12:48:24 -0700 (Fri, 01 Mar 2013)</p><p>Addition of full 3D operators to tensor operators design document.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/documents/shared/current_design_doc/tensor_operations/tensor_operations.tex
===================================================================
--- trunk/documents/shared/current_design_doc/tensor_operations/tensor_operations.tex 2013-03-01 19:37:33 UTC (rev 2536)
+++ trunk/documents/shared/current_design_doc/tensor_operations/tensor_operations.tex 2013-03-01 19:48:24 UTC (rev 2537)
@@ -67,7 +67,7 @@
\chapter{Algorithmic Formulations}
-\section{Design Solution: Strain rate tensor}
+\section{Design Solution: Strain rate tensor, projection}
Date last modified: 2/12/2013 \\
Contributors: Mark Petersen \\
@@ -108,7 +108,7 @@
\label{fig:2D domain}
\end{figure}
-\subsection{Discrete strain rate tensor at cell centers}
+\subsection{Discrete 2D strain rate tensor at cell centers}
Here we compute the strain rate tensor at a cell center ${\bf x}_i$ based on the normal and tangential velocities, ${\bf u}_e$ and ${\bf v}_e$, at the edges of that cell (Figure \ref{fig:edge vectors}). The total velocity at edge $e$ is $({\bf u}_e + {\bf v}_e)$. Discretizing the weak form (\ref{weak strain rate tensor}), the strain rate tensor at cell $i$ is
\begin{equation}
@@ -163,20 +163,21 @@
where $(sym)$ indicates that $\dot{\epsilon}_{21}=\dot{\epsilon}_{12}$ on this symmetric matrix.
-\subsection{Discrete strain rate tensor at vertices}
+\subsection{Discrete 2D strain rate tensor at vertices}
The strain rate computation at vertices is similar to that for cells, except that the control area is now a dual-mesh cell centered on vertex at ${\bf x}_v$, bounded by lines of length $d_e$ connecting primal-mesh cell centers \cite[figure 1]{Ringler_ea10jcp}. The normal vector to an edge on the dual-mesh is
\begin{eqnarray}
&&\tilde{\bf n}_e = \left( \begin{array}{c} \tilde{n}_{e1} \\ \tilde{n}_{e2} \end{array} \right)
= \left( \begin{array}{r} \sin{\theta_{e}} \\ -\cos{\theta_{e}} \end{array} \right),
\end{eqnarray}
-so that $\tilde{\bf n}_e$ is rotated clockwise from ${\bf n}_e$ by 90 degrees. The projects of ${\bf u}_e$ and ${\bf v}_e$ remain the same as in (\ref{u_e}--\ref{v_e}). The strain rate tensor at a vertex is then
+so that $\tilde{\bf n}_e$ is rotated clockwise from ${\bf n}_e$ by 90 degrees. The projections of ${\bf u}_e$ and ${\bf v}_e$ remain the same as in (\ref{u_e}--\ref{v_e}). The strain rate tensor at a vertex is then
\begin{eqnarray}
&\left[</font>
<font color="gray">abla_s {\bf u} \right]_v
-= \frac{1}{A_v}
+= \displaystyle \frac{1}{A_v}
\sum_{e\in EV(v)}{\frac{1}{2}\left(
\tilde{\bf n}_e\otimes\left({\bf u}_e +{\bf v}_e \right)
+ \left({\bf u}_e +{\bf v}_e \right)\otimes\tilde{\bf n}_e\right)l_e },
+\label{vertex strain rate tensor2}
\\
& = \displaystyle
\frac{1}{A_v}
@@ -193,13 +194,72 @@
& -v_e \cos^2\theta_e - u_e \cos{\theta_e} \sin{\theta_e} \end{array} \right)l_e
\end{eqnarray}
+\subsection{Discrete 3D strain rate tensor at cell centers}
+Here we present the full three-dimensional calculation for a strain rate tensor. This version avoids the use of $\theta_e$, the \verb|angleEdge| variable, altogether, and so does not have any issues near the pole. However, this is at the expense of going from a 2x2 symmetric tensor to a 3x3 symmetric tensor. It is also not clear how to generalize the 2D tensor operations in the CICE dynamics \cite[section 3.4]{CICE_manual_4.1}.
+Beginning with (\ref{cell strain rate tensor2})
+\begin{equation}
+\label{cell strain rate tensor2}
+\left[</font>
<font color="blue">abla_s {\bf u} \right]_i
+=
+\frac{1}{A_i}
+\sum_{e\in EC(i)}{\frac{1}{2}\left(
+ {\bf n}_e\otimes\left({\bf u}_e +{\bf v}_e \right)
++ \left({\bf u}_e +{\bf v}_e \right)\otimes{\bf n}_e\right)l_e },
+\end{equation}
+we rewrite the normal and tangential velocities in terms of a scalar coefficient and a unit vector,
+\begin{eqnarray}
+{\bf u}_e = u_e {\bf n}_e,
+ && {\bf n}_e = \frac{{\bf x}_{i2}-{\bf x}_{i1}}{\left||{\bf x}_{i2}-{\bf x}_{i1}\right||}
+ = \frac{1}{d_e} \sum_{i\in CE(e)}n_{e,i}{\bf x}_{i}
+\\
+{\bf v}_e = v_e \tilde{\bf n}_e,
+ && \tilde{\bf n}_e = \frac{{\bf x}_{v2}-{\bf x}_{v1}}{\left||{\bf x}_{v2}-{\bf x}_{v1}\right||}
+ = \frac{1}{d_v} \sum_{v\in VE(e)}n_{e,i}{\bf x}_{v}?
+\end{eqnarray}
+I have not settled on the best notation yet, but ${\bf n}_e$ is the unit vector normal to the edge, and $\tilde{\bf n}_e $ is the unit vector tangent to the edge, pointing in the same direction as ${\bf v}_e$. Continuing from (\ref{cell strain rate tensor2}),
+\begin{equation}
+\label{cell strain rate tensor3}
+\left[</font>
<font color="blue">abla_s {\bf u} \right]_i
+=
+\frac{1}{A_i}
+\sum_{e\in EC(i)}{\frac{1}{2}\left(
+ 2u_e{\bf n}_e\otimes{\bf n}_e
++ v_e{\bf n}_e\otimes\tilde{\bf n}_e
++ v_e\tilde{\bf n}_e\otimes{\bf n}_e
+\right)l_e }.
+\end{equation}
+The unit vectors ${\bf n}_e$ and $\tilde{\bf n}_e$ are 3-vectors in $(x,y,z)$ space, making the tensor products 3x3 tensors.
+
+
+\subsection{Discrete 3D strain rate tensor at vertices}
+
+Beginning with (\ref{vertex strain rate tensor3}),
+\begin{eqnarray}
+\label{vertex strain rate tensor3}
+\left[</font>
<font color="blue">abla_s {\bf u} \right]_v
+&=& \displaystyle \frac{1}{A_v}
+\sum_{e\in EV(v)}{\frac{1}{2}\left(
+ \tilde{\bf n}_e\otimes\left({\bf u}_e +{\bf v}_e \right)
++ \left({\bf u}_e +{\bf v}_e \right)\otimes\tilde{\bf n}_e\right)l_e },
+\\
+&=& \displaystyle \frac{1}{A_v}
+\sum_{e\in EV(v)}{\frac{1}{2}\left(
+ 2v_e\tilde{\bf n}_e\otimes\tilde{\bf n}_e
++ u_e{\bf n}_e\otimes\tilde{\bf n}_e
++ u_e\tilde{\bf n}_e\otimes{\bf n}_e
+\right)l_e },
+\end{eqnarray}
+so that a 3x3 strain rate tensor is computed for each vertex.
+
</font>
<font color="blue">ewpage
\section{Design Solution: divergence of a tensor}
Date last modified: 2/12/2013 \\
Contributors: Mark Petersen \\
+\subsection{2D divergence of a tensor}
+
In two-dimensional Cartesian coordinates, the component of the divergence of a tensor $\sigma$ in the direction of a unit vector ${\bf n}$ is
\begin{eqnarray}
\left(</font>
<font color="gray">abla\cdot\sigma\right)\cdot{\bf n}
@@ -221,6 +281,7 @@
In most cases we are interested in projecting the divergence of a tensor onto ${\bf n}_e$, the normal to a cell edge. For example, the momentum equation for sea-ice requires this term when computing the tendencies of $u_e$ \cite[section 3.4]{CICE_manual_4.1}. It is therefore convenient to arrange a local Cartesian coordinate system where the 1-direction is aligned with ${\bf n}_e$ (i.e. pointing from one cell center to the other) and the 2-direction is aligned with $\tilde{\bf n}_e$ (i.e. pointing from one vertex to the other). This is possible because on a Voronoi tesselation grid, the line joining cell centers across a cell edge is always orthogonal to that cell edge. A tensor may be rotated through an angle $\theta$ into a new coordinate system using rotation matrices as
\begin{eqnarray}
+\label{rotate sigma 2D}
\tilde{\sigma} &=& R_\theta \sigma R_\theta^*,\\
&=&
\left(\begin{array}{cc} \cos{\theta}& -\sin{\theta} \\
@@ -273,11 +334,37 @@
{\it Response, Mark:} I only know the sigma tensor at cell centers and vertices. How could I compute \ref{sigma1} for any n in the 1,2 reference frame using strong derivatives, as in \ref{discrete div1}? The only solution I can think of is to rotate the tensor so that the directions to take derivatives aligns parallel and perpendicular to the local edge. This then conforms with our standard methods of finite differences.
+\subsection{3D divergence of a tensor}
+
+The general form required is the divergence of a tensor normal to edge $e$,
+\begin{equation}
+(</font>
<font color="blue">abla\cdot\sigma)\cdot{\bf n}_e,
+\end{equation}
+is coordinate free, and may be computed in three dimensions. The difficulty is that we only have information at the cell centers and vertices, and so must take strong derivatives in those directions. This is the same issue that required rotating the tensor in 2D in (\ref{rotate sigma 2D}). This time we rotate in 3D,
+\begin{eqnarray}
+\tilde{\sigma}_{i,e} &=& R_e\sigma_i R_e^*,
+\\
+\tilde{\sigma}_{v,e} &=& R_e\sigma_v R_e^*,
+\\
+R_e &=& \left[
+\begin{array}{ccc}
+{\bf n}_e & \tilde{\bf n}_e & {\bf n}_e\times\tilde{\bf n}_e
+\end{array}
+\right]
+\\
+\left[(</font>
<font color="gray">abla\cdot\sigma)\cdot{\bf n}_e\right]_e &=& \tilde\sigma_{11,1} + \tilde\sigma_{21,2},
+\\
+ &=& \frac{1}{d_e} \sum_{i\in CE(e)}-n_{e,i}\left[\tilde\sigma_{11}\right]_{i,e}
+ + \frac{1}{l_e} \sum_{v\in VE(e)}-t_{e,v}\left[\tilde\sigma_{12}\right]_{v,e}
+\end{eqnarray}
+where the 1-direction points from cell to cell with ${\bf n}_e$, the 2-direction points from vertex to vertex with $\tilde{\bf n}_e$, and the 3-direction is perpendicular to those directions, using ${\bf n}_e\times\tilde{\bf n}_e$. The rotation matrix $R_e$ is simply composed of those three vectors, and the operation $R_e\sigma_i R_e^*$ rotates the original tensor in an $(x,y,z)$ reference frame to the new reference frame. All components of the rotated stress tensor with a 3 index ($\tilde\sigma_{13}$, $\tilde\sigma_{23}$, $\tilde\sigma_{33}$, etc.) should be small because velocities are tangent to the earth, but may not be exactly zero due to the curvature of the earth from one cell to the next. For domains in the $x-y$ plane, ${\bf n}_e\times\tilde{\bf n}_e=(0,0,1)$ for all edges, and all elements of $\tilde\sigma_{ij}$ with an index value of three will be exactly zero.
+
+
%-----------------------------------------------------------------------
\chapter{Design and Implementation}
-\section{Implementation: strain rate tensor}
+\section{Implementation: 2D strain rate tensor}
Date last modified: 2/12/2013 \\
Contributors: Mark Petersen \\
@@ -358,7 +445,7 @@
</font>
<font color="red">ewpage
-\section{Implementation: divergence of a tensor}
+\section{Implementation: 2D divergence of a tensor}
Date last modified: 2/12/2013 \\
Contributors: Mark Petersen \\
</font>
</pre>