mirror of https://github.com/Askill/claude.git
Merge pull request #12 from TechWizzart/master
Manual update: Potential Temperature and Multidimensional FFTs
This commit is contained in:
Simon Clark
GitHub
commit
2a361a63ad
BIN
CLAuDE NOM.pdf
BIN
CLAuDE NOM.pdf
Binary file not shown.
|
|
@ -10,7 +10,12 @@
|
|||
\usepackage{float}
|
||||
\usepackage{verbatim}
|
||||
\usepackage[normalem]{ulem}
|
||||
\usepackage{graphicx}
|
||||
\usepackage{siunitx}
|
||||
|
||||
\setcounter{section}{-1}
|
||||
\newtheorem{lemma}{Lemma}
|
||||
\newcommand{\lemmaautorefname}{Lemma}
|
||||
|
||||
\title{CLimate Analysis using Digital Estimations Non-Offical Manual (CLAuDE NOM)}
|
||||
\author{Sam "TechWizard" Baggen}
|
||||
|
|
@ -22,6 +27,8 @@
|
|||
|
||||
\tableofcontents
|
||||
|
||||
\newpage
|
||||
|
||||
\section{Introduction}
|
||||
|
||||
The CLimate Analysis using Digital Estimations model is a simplified planetary climate model. It will be used to educate people on how climate physics works and to experiment with different
|
||||
|
|
|
|||
|
|
@ -55,4 +55,5 @@ only applies to variables in code, every symbol in equations are explained at th
|
|||
\item $smooth_u$: The smoothing parameter for the $u$ component of the velocity.
|
||||
\item $smooth_v$: The smoothing parameter for the $v$ component of the velocity.
|
||||
\item $smooth_w$: The smoothing parameter for the $w$ component of the velocity.
|
||||
\item $smooth_{vert}$: The smoothing parameter for the vertical part of the velocities.
|
||||
\end{itemize}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 202 KiB |
|
|
@ -149,6 +149,18 @@ chapter={1},
|
|||
pages={12}
|
||||
}
|
||||
|
||||
@inbook{thermalPotential,
|
||||
place={Cambridge},
|
||||
edition={2nd},
|
||||
DOI={10.1017/9781107588417.003},
|
||||
title={Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation},
|
||||
publisher={Cambridge University Press},
|
||||
author={Vallis, Geoffrey K.},
|
||||
year={2017},
|
||||
chapter={1},
|
||||
pages={24}
|
||||
}
|
||||
|
||||
%Calculus
|
||||
@inbook{areaCircle,
|
||||
place={Don Mills, Ont.},
|
||||
|
|
@ -384,3 +396,13 @@ month={Aug}
|
|||
year={2011},
|
||||
publisher={IEEE Computer Society}
|
||||
}
|
||||
|
||||
@inbook{cancellation,
|
||||
author = {Thomas H. Cormen and Charles E and Leiserson and Ronald L. Rivest and Clifford Stein},
|
||||
title = {Introduction to Algorithms},
|
||||
edition = {3rd},
|
||||
chapter = {30},
|
||||
pages = {907},
|
||||
publisher = {MIT Press},
|
||||
year = {2009}
|
||||
}
|
||||
|
|
@ -36,9 +36,9 @@ With thermal diffusion in place, the temperature will spread out a bit, however
|
|||
going to change that. The advection equation is shown in \autoref{eq:advection}. The symbols are:
|
||||
|
||||
\begin{itemize}
|
||||
\item $\psi$: What is carried along (in our case temperature, $K$).
|
||||
\item $t$: The time ($s$).
|
||||
\item $u$: The fluid velocity vector ($ms^{-1}$).
|
||||
\item $\psi$: What is carried along (in our case temperature, \si{K}).
|
||||
\item $t$: The time (\si{s}).
|
||||
\item $u$: The fluid velocity vector (\si{ms^{-1}}).
|
||||
\item $\nabla$: The divergence operator (as explained in \autoref{sec:laplace}).
|
||||
\end{itemize}
|
||||
|
||||
|
|
@ -54,7 +54,7 @@ With the divergence functon defined in \autoref{alg:divergence}, we now need to
|
|||
$T_{add} \leftarrow T_a + \delta t \alpha_a \nabla^2(T_a) + \nabla(T_a)$ \;
|
||||
$T_a \leftarrow T_a + T_{add}[5:-5, :] \text{ //Only add } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + 5, nlat - 5]$. \;
|
||||
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
|
||||
\caption{The main calculations for calculating the effects of diffusion}
|
||||
\caption{The main calculations for calculating the effects of advection}
|
||||
\label{alg:advection}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -78,7 +78,7 @@ to calculate $\rho$ after the movement of air has taken place, so we need to cha
|
|||
$T_a \leftarrow T_a + T_{add}[5:-5, :] \text{ //Only add } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + 5, nlat - 5]$. \;
|
||||
$\rho \leftarrow \rho + \delta t \nabla \rho$ \;
|
||||
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
|
||||
\caption{The main calculations for calculating the effects of diffusion}
|
||||
\caption{The main calculations for calculating the effects of advection}
|
||||
\label{alg:advectionv2}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -92,7 +92,7 @@ boundary.
|
|||
$T_a \leftarrow T_a - 0.5T_{add}[adv\_bound:-adv\_boun, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$. \;
|
||||
$\rho[adv\_boun: -adv\_boun, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$ \;
|
||||
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
|
||||
\caption{The main calculations for calculating the effects of diffusion}
|
||||
\caption{The main calculations for calculating the effects of advection}
|
||||
\label{alg:advectionfix}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -105,7 +105,7 @@ add it, with \autoref{alg:advection layer} as a result. Here the ':' means all i
|
|||
$T_a \leftarrow T_a - 0.5T_{add}[adv\_boun:-adv\_boun, :, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$. \;
|
||||
$\rho[adv\_boun: -adv\_boun, :, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$ \;
|
||||
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
|
||||
\caption{The main calculations for calculating the effects of diffusion}
|
||||
\caption{The main calculations for calculating the effects of advection}
|
||||
\label{alg:advection layer}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -113,4 +113,112 @@ First thing to mention is that vertical advection is still broken. Why? Because
|
|||
we calculate the difference from one layer to the other is by differencing them (subtracting) which is always finite. Therefore we always get some inaccuracies. Usually that is fine, but with an
|
||||
exponential function the differences, you guessed it, become exponentially wrong. As such, the function would eventually be so far off that the model would blow up. So we need to fix it. To
|
||||
prevent a blow up, we have disabled the call to the gradient $z$ function in \autoref{alg:divergence}. This ensures that the horizontal bits still work, but the vertical stuff does not.
|
||||
As always, we will try to fix this in a future stream.
|
||||
As always, we will try to fix this in the future.
|
||||
|
||||
\subsection{Adiabatic Motion}
|
||||
Up until now, we have been moving air and the density of the atmosphere. However we have not transported any temperature yet. What this means is that if we have a packet of air (see it as a box
|
||||
filled with air, but the box is invisible) $P$ which is at the planet surface. There $P$ has temperature $T_1$. If we then move $P$ to a layer higher up in the atmosphere, $P$ will still have
|
||||
the same temperature, which is wrong because the density differs. Due to this difference, there is either more or less pressure applied to the box of air which means that $P$ will contract or
|
||||
expand. This obviously changes it's temperature as the air molecules are closer together/further away from each other. Therefore the energy spread is different which affects the temperature.
|
||||
For a visual representation, please consult \autoref{fig:thermal potential}.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=0.9\textwidth]{figures/potential_temperature.jpg}
|
||||
\caption{Visual representation of why we need Thermal Potential \cite{simon}.}
|
||||
\label{fig:thermal potential}
|
||||
\end{figure}
|
||||
|
||||
As seen in \autoref{fig:thermal potential}, the packet of air (next to $T_1$) moves up into a higher layer of the atmosphere and ends up at the top (next to $T_2$). The packet has grown bigger,
|
||||
as the density in that atmospheric layer is lower and hence the air expands. Because the same energy is in that packet, but it has expanded, the temperature of that packet drops which gives us
|
||||
$T_2 < T_1$. The light blue graph in the background shows how the density behaves the higher you go, meaning at the far right (when we are at the highest point) the density is quite low and at
|
||||
the far left (when we are at the planet surface) the density is quite high. The yellow graph shows the temperature of the air at that height, which behaves similarly to the density graph.
|
||||
|
||||
Thermal potential or potential temperature, is the temperature that the packet of air would be if it is at the planet surface. So take the packet at $T_2$ in \autoref{fig:thermal potential}. If we move that down to the
|
||||
surface, then it would be temperature $T$ which is then equal to the thermal potential. The equation corresponding to thermal potential is shown in \autoref{eq:thermal potential}
|
||||
\cite{thermalPotential}. The symbols in \autoref{eq:thermal potential} mean:
|
||||
|
||||
\begin{itemize}
|
||||
\item $T$: The temperature of the packet of air (\si{K}).
|
||||
\item $p_0$: Reference pressure, usually the pressure at the planet surface (\si{Pa}).
|
||||
\item $p$: Pressure of the packet of air (\si{Pa}).
|
||||
\item $R$: Gas constant as defined in \autoref{sec:gas constant} (\si{JK^{-1}mol^{-1}}).
|
||||
\item $C_a$: Specific heat capacity of air at a constant pressure (\si{Jkg^{-1}K^{-1}}).
|
||||
\end{itemize}
|
||||
|
||||
\begin{subequations}
|
||||
\begin{equation}
|
||||
\theta = T(\frac{p_0}{p})^{\frac{R}{C_a}}
|
||||
\label{eq:thermal potential}
|
||||
\end{equation}
|
||||
\begin{equation}
|
||||
T = \theta(\frac{p}{p_0})^{\frac{R}{C_a}}
|
||||
\label{eq:potential temp}
|
||||
\end{equation}
|
||||
\end{subequations}
|
||||
|
||||
If we now re-arrange \autoref{eq:thermal potential} so that we have $T$ on one side and the rest on the other side, we get \autoref{eq:potential temp}. With this we can convert temperature into
|
||||
potential temperature and vice versa. The whole process of moving temperature around is called adiabatic motion. Now it is time to get this into code. For this to work we need to translate both
|
||||
\autoref{eq:thermal potential} and \autoref{eq:potential temp} into code and create an overarching function that calls the previously mentioned equations. Let us start with
|
||||
\autoref{eq:thermal potential} which is described in \autoref{alg:temp to pot}. Note that $\frac{R}{C_a}$ does not change as they are constants, therefore we can precompute them which saves quite
|
||||
a bit of time (namely $O(n^3)$ divisions, where $n$ is the length of each dimension of $T_a$). Also note that we can inverse the process by inserting a minus in the exponent and swapping $T_a$
|
||||
and $\theta$, which can easily be done in the call to \autoref{alg:temp to pot}. To avoid confusion we rename $p_0$ to $p_r$ as we talk about a reference pressure, instead of the already defined
|
||||
pressure from a previous calculation round.
|
||||
|
||||
\begin{algorithm}
|
||||
\SetKwInOut{Input}{Input}
|
||||
\SetKwInOut{Output}{Output}
|
||||
\Input{temperature of the atmosphere $T_a$, air pressure $p$, boolean $back$}
|
||||
\Output{potential temperature $\theta$}
|
||||
\uIf{$back$}{
|
||||
$\kappa \leftarrow -\frac{R}{C_a}$ \;
|
||||
}\uElse{
|
||||
$\kappa \leftarrow \frac{R}{C_a}$ \;
|
||||
}
|
||||
\For{$i \leftarrow 0$ \KwTo $T_a.length$}{
|
||||
\For{$j \leftarrow 0$ \KwTo $T_a[i].length$}{
|
||||
$p_r \leftarrow p[i, j 0]$ \;
|
||||
\For{$k \leftarrow 0$ \KwTo $T_a[i, j].length$}{
|
||||
$\theta[i, j, k] \leftarrow T_a[i, j, k] (\frac{p[i, j, k]}{p_r})^{\kappa}$
|
||||
}
|
||||
}
|
||||
}
|
||||
\Return{$\theta$}
|
||||
\caption{Converting temperature into potential temperature}
|
||||
\label{alg:temp to pot}
|
||||
\end{algorithm}
|
||||
|
||||
With that done we now only need to create the overarching function that will call \autoref{alg:temp to pot}. This function will be based on \autoref{alg:divergence} as it is quite similar for
|
||||
what we need to do, however we need to replace a call so we do need to make it a seperate function. Let us do that in \autoref{alg:thermal pot}. Here \texttt{TToTheta} represents the call to
|
||||
\autoref{alg:temp to pot}.
|
||||
|
||||
\begin{algorithm}[!hbt]
|
||||
\SetKwInOut{Input}{Input}
|
||||
\SetKwInOut{Output}{Output}
|
||||
\Input{temperature of the atmosphere $T_a$, pressure $p$}
|
||||
\Output{temperature of the atmosphere after advection $output$}
|
||||
$\theta \leftarrow $ \texttt{TToTheta}($T_a, p,$ \texttt{FALSE}) \;
|
||||
\For{$i \leftarrow 0$ \KwTo $a.length$}{
|
||||
\For{$j \leftarrow 0$ \KwTo $a[i].length$}{
|
||||
\For{$k \leftarrow 0$ \KwTo $a[i, j].length$}{
|
||||
$output[i, j, k] \leftarrow \Delta_x(\theta u, i, j, k) + \Delta_y(\theta v, i, j, k) + \Delta_z(\theta w, i, j, k)$ \;
|
||||
}
|
||||
}
|
||||
}
|
||||
$output \leftarrow $ \texttt{TToTheta}($output, p,$ \texttt{TRUE}) \;
|
||||
\Return{$output$} \;
|
||||
\caption{Calculate the result of the thermal advection}
|
||||
\label{alg:thermal pot}
|
||||
\end{algorithm}
|
||||
|
||||
Now we only need to do one more thing, replace the algorithm that calculates $T_a$ as a result of advection. Let us do that in \autoref{alg:advection pot}, where \texttt{ThermalAdv} is the call
|
||||
to \autoref{alg:thermal pot}.
|
||||
|
||||
\begin{algorithm}
|
||||
$T_{add} \leftarrow T_a + $ \texttt{ThermalAdv}$(T_a, p)$\;
|
||||
$T_a \leftarrow T_a - 0.5T_{add}[adv\_boun:-adv\_boun, :, :] \text{ //Only subtract } T_{add} \text{ to } T_a \text{ for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$. \;
|
||||
$\rho[adv\_boun: -adv\_boun, :, :] \leftarrow \rho - 0.5(\delta t \nabla \rho) \text{ //Only change the density for indices in the interval } [-nlat + adv\_boun, nlat - adv\_boun]$ \;
|
||||
$T_p \leftarrow T_p + \delta t \alpha_p \nabla^2(T_p)$ \;
|
||||
\caption{The main calculations for calculating the effects of advection}
|
||||
\label{alg:advection pot}
|
||||
\end{algorithm}
|
||||
|
|
@ -11,7 +11,7 @@ do that we start with the fixed part of the Control Panel, the physical constant
|
|||
gravity is, how fast they spin around their axis and many many more. What does not change are the physical constants, well because they are constant. The Stefan-Boltzmann constant for instance
|
||||
does not change. Whether you are on Earth, in space or on Jupiter, the value of the Stefan-Boltzmann constant will remain the same.
|
||||
|
||||
The Stefan-Boltzmann constant is denoted by $\sigma$ and has a value of $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}. The $\sigma$ is a greek letter called sigma. Greek
|
||||
The Stefan-Boltzmann constant is denoted by $\sigma$ and has a value of $5.670373 \cdot 10^-8$ (\si{Wm^{-2}K^{-4}}) \cite{stefan-boltzmann}. The $\sigma$ is a greek letter called sigma. Greek
|
||||
letters are often used in mathematics, as well as in physics or any other discipline that relies on maths (spoiler alert, quite a lot). Treat it like a normal letter in maths, representing a
|
||||
number that you either do not know yet or is too long or cumbersome to write down every time. The Stefan-Boltzmann constant is denoted in scientific notation, a number followed by the order of
|
||||
magnitude. It is denoted as a multiplication, because that is what you have to do to get the real number. An example: $4.3 \cdot 10^2 = 430$ and $4.3 \cdot 10^{-2} = 0.043$. The
|
||||
|
|
@ -21,29 +21,29 @@ letters you will find the following: [number]. This is a citation, a reference t
|
|||
citation to show that I am not making these numbers up. This is what scientists use to back up their claims if they do not want to redo the work that others have done. I mean what is the point
|
||||
of re-inventing the wheel if there is a tyre company next door? That is why scientists citate.
|
||||
|
||||
So with that out of the way, let us write down some constants. Why do I do this here? Because a lot of constants are used everywhere and I am too lazy to relicate them every time. If you see a
|
||||
So with that out of the way, let us write down some constants. Why do I do this here? Because a lot of constants are used everywhere and I am too lazy to replicate them every time. If you see a
|
||||
letter or symbol that is not explicitly explained, then it is most likely a constant that we discuss here in the control panel.
|
||||
|
||||
\subsection{Physical Constants}
|
||||
As mentioned before, physical constants do not change based on where you are in the universe. Below you will find an overview of all the relevant constants together with their units. And a short
|
||||
explanation where they are used or what they represent. To see them in action, consult the other sections of this manual, you will find them in equations all throughout this document.
|
||||
|
||||
\subsubsection{The Gas Constant}
|
||||
The Gas constant, $R = 8.3144621$ ($J(mol)^{-1}K$) \cite{idealGas} is the constant used to relate the temperature of the gas to the pressure and the volume. One would expect this constant to be
|
||||
\subsubsection{The Gas Constant} \label{sec:gas constant}
|
||||
The Gas constant, $R = 8.3144621$ (\si{JK^{-1}mol^{-1}}) \cite{idealGas} is the constant used to relate the temperature of the gas to the pressure and the volume. One would expect this constant to be
|
||||
different per gas, but under high enough temperatures and low enough pressure the gas constant is the same for all gases.
|
||||
|
||||
\subsubsection{The Specific Heat Capacity}
|
||||
The specific heat capacity $c$ depicts how much energy is required to heat the object by one degree Kelvin per unit mass ($\frac{J}{Kg \cdot K}$) \cite{specificHeat}. This varies per material
|
||||
and is usually indicated by a subscript. The specific heat capacity for water for instance is $c_w = 4190 JKg^{-1}K^{-1}$. Specific heat capacities also exist in the form of $Jg^{-1}K^{-1}$,
|
||||
$Jmol^{-1}K^{-1}$ and $Jcm^{-3}K^{-1}$ which you can use in various circumstances, depending on what information you have.
|
||||
The specific heat capacity $c$ depicts how much energy is required to heat the object by one degree Kelvin per unit mass (\si{Jkg^{-1}K^{-1}}) \cite{specificHeat}. This varies per material
|
||||
and is usually indicated by a subscript. The specific heat capacity for water for instance is $c_w = 4190$ \si{Jkg^{-1}K^{-1}}. Specific heat capacities also exist in the form of
|
||||
\si{Jkg^{-1}K^{-1}}, \si{Jmol^{-1}K^{-1}} and \si{Jcm^{-3}K^{-1}} which you can use in various circumstances, depending on what information you have.
|
||||
|
||||
\subsubsection{Mole}
|
||||
Mole is the amount of particles ($6.02214076 \cdot 10^{23}$) in a substance, where the average weight of one mole of particles in grams is about the same as the weight of one particle in atomic
|
||||
mass units ($u$)\cite{mole}. This is not a physical constant perse, but more like a unit ($mol$). Though it is still important enough to be added here for future reference. All other units are
|
||||
mass units (\si{u})\cite{mole}. This is not a physical constant perse, but more like a unit (\si{mol}). Though it is still important enough to be added here for future reference. All other units are
|
||||
way more intuitive and are assumed to be known.
|
||||
|
||||
\subsubsection{The Stefan-Boltzmann Constant}
|
||||
The Stefan-Boltzmann constant, $\sigma = 5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann} is used in the Stefan-Boltzmann law (more on that in \autoref{sec:first thermolaw}).
|
||||
The Stefan-Boltzmann constant, $\sigma = 5.670373 \cdot 10^-8$ \si{Wm^{-2}K^{-4})} \cite{stefan-boltzmann} is used in the Stefan-Boltzmann law (more on that in \autoref{sec:first thermolaw}).
|
||||
|
||||
\subsection{Planet Specific Variables}
|
||||
The following set of variables vary per planet, that's why we call them variables since they vary. Makes sense right? We add them here as we will use them throughout the manual. The advantage
|
||||
|
|
@ -62,9 +62,9 @@ means less work for the computer. Yes computers are lazy too.
|
|||
\caption{Definition of how much time it takes for a day and a year on a planet and how much time on the planet passes before we start another calculation run}
|
||||
\label{alg:time constants}
|
||||
\SetKwComment{Comment}{//}{}
|
||||
$day \leftarrow 60*60*24$ \Comment*[l]{Length of one day in seconds ($s$)}
|
||||
$year \leftarrow 365*day$ \Comment*[l]{Length of one year in seconds ($s$)}
|
||||
$\delta t \leftarrow 60 * 9$ \Comment*[l]{How much time is between each calculation run in seconds ($s$)}
|
||||
$day \leftarrow 60*60*24$ \Comment*[l]{Length of one day in seconds (\si{s})}
|
||||
$year \leftarrow 365*day$ \Comment*[l]{Length of one year in seconds (\si{s})}
|
||||
$\delta t \leftarrow 60 * 9$ \Comment*[l]{How much time is between each calculation run in seconds (\si{s})}
|
||||
\end{algorithm*}
|
||||
|
||||
\subsubsection{The Planet Passport}
|
||||
|
|
@ -75,12 +75,12 @@ all the relevant variables that are unique to a planet, or well not necessarily
|
|||
\caption{Defining the constants that are specific to a planet}
|
||||
\label{alg:planet constants}
|
||||
\SetKwComment{Comment}{//}{}
|
||||
$g \leftarrow 9.81$ \Comment*[l]{Magnitude of gravity on the planet in $ms^{-2}$}
|
||||
$g \leftarrow 9.81$ \Comment*[l]{Magnitude of gravity on the planet in \si{ms^{-2}}}
|
||||
$\alpha \leftarrow -23.5$ \Comment*[l]{By how many degrees the planet is tilted with respect to the star's plane}
|
||||
$top \leftarrow 50*10^3$ \Comment*[l]{How high the top of the atmosphere is with respect to the planet surface in meters ($m$)}
|
||||
$ins \leftarrow 1370$ \Comment*[l]{Amount of energy from the star that reaches the planet per unit area ($Jm^{-2}$)}
|
||||
$top \leftarrow 50*10^3$ \Comment*[l]{How high the top of the atmosphere is with respect to the planet surface in meters (\si{m})}
|
||||
$ins \leftarrow 1370$ \Comment*[l]{Amount of energy from the star that reaches the planet per unit area (\si{Jm^{-2}})}
|
||||
$\epsilon \leftarrow 0.75$ \Comment*[l]{Absorbtivity of the atmosphere, fraction of how much of the total energy is absorbed (unitless)}
|
||||
$r \leftarrow 6.4*10^6$ \Comment*[l]{The radius of the planet in meters ($m$)}
|
||||
$r \leftarrow 6.4*10^6$ \Comment*[l]{The radius of the planet in meters (\si{m})}
|
||||
\end{algorithm}
|
||||
|
||||
\subsubsection{Model Specific Parameters}
|
||||
|
|
@ -97,12 +97,12 @@ definitions can be found in \autoref{alg:model constants}. What the $adv$ boolea
|
|||
$resolution \leftarrow 3$ \Comment*[l]{The amount of degrees on the latitude longitude grid that each cell has, with this setting each cell is 3 degrees latitude high and 3 degrees
|
||||
longitude wide}
|
||||
$nlevels \leftarrow 10$ \Comment*[l]{The amount of layers in the atmosphere}
|
||||
$\delta t_s \leftarrow 60*137$ \Comment*[l]{The time between calculation rounds during the spin up period in seconds ($s$)}
|
||||
$t_s \leftarrow 5*day$ \Comment*[l]{How long we let the planet spin up in seconds ($s$)}
|
||||
$\delta t_s \leftarrow 60*137$ \Comment*[l]{The time between calculation rounds during the spin up period in seconds (\si{s})}
|
||||
$t_s \leftarrow 5*day$ \Comment*[l]{How long we let the planet spin up in seconds (\si{s})}
|
||||
$adv \leftarrow \texttt{FALSE}$ \Comment*[l]{Whether we want to enable advection or not}
|
||||
$adv\_boun \leftarrow 8$ \Comment*[l]{How many cells away from the poles where we want to stop calculating the effects of advection}
|
||||
$C_a \leftarrow 287$ \Comment*[l]{Heat capacity of the atmosphere in $JKg^{-1}K^{-1}$}
|
||||
$C_p \leftarrow 1 \cdot 10^6$ \Comment*[l]{Heat capacity of the planet in $JKg^{-1}K^{-1}$}
|
||||
$C_a \leftarrow 287$ \Comment*[l]{Heat capacity of the atmosphere in \si{Jkg^{-1}K^{-1}}}
|
||||
$C_p \leftarrow 1 \cdot 10^6$ \Comment*[l]{Heat capacity of the planet in \si{Jkg^{-1}K^{-1}}}
|
||||
$\delta y \leftarrow \frac{2\pi r}{nlat}$ \Comment*[l]{How far apart the gridpoints in the y direction are (degrees latitude)}
|
||||
$\alpha_a \leftarrow 2 \cdot 10^{-5}$ \Comment*[l]{The diffusivity constant for the atmosphere}
|
||||
$\alpha_p \leftarrow 1.5 \cdot 10^{-6}$ \Comment*[l]{The diffusivity constant for the planet surface}
|
||||
|
|
@ -110,9 +110,10 @@ definitions can be found in \autoref{alg:model constants}. What the $adv$ boolea
|
|||
$smooth_u \leftarrow 0.8$ \Comment*[l]{The smoothing parameter for the $u$ component of the velocity}
|
||||
$smooth_v \leftarrow 0.8$ \Comment*[l]{The smoothing parameter for the $v$ component of the velocity}
|
||||
$smooth_w \leftarrow 0.3$ \Comment*[l]{The smoothing parameter for the $w$ component of the velocity}
|
||||
$smooth_{vert} \leftarrow 0.3$ \Comment*[l]{The smoothing parameter for the vertical part of the velocities}
|
||||
$count \leftarrow 0$ \;
|
||||
\For{$j \in [0, top]$}{
|
||||
$heights[j] \leftarrow count$ \Comment*[l]{The height of a layer}
|
||||
$heights[j] \leftarrow count$ \Comment*[l]{The height of a layer (\si{m})}
|
||||
$count \leftarrow count + \frac{top}{nlevels}$ \;
|
||||
}
|
||||
|
||||
|
|
@ -121,6 +122,6 @@ definitions can be found in \autoref{alg:model constants}. What the $adv$ boolea
|
|||
}
|
||||
|
||||
\For{$k \in [0, nlevels - 1]$}{
|
||||
$\delta z[k] \leftarrow heights[k + 1] - heights$ \Comment*[l]{How far apart the gridpoints in the z direction are ($m$)}
|
||||
$\delta z[k] \leftarrow heights[k + 1] - heights$ \Comment*[l]{How far apart the gridpoints in the z direction are (\si{m})}
|
||||
}
|
||||
\end{algorithm}
|
||||
|
|
@ -21,7 +21,7 @@ calculations in \autoref{alg:temperature with density} we would calculate the ve
|
|||
}
|
||||
}
|
||||
}
|
||||
\caption{Main loop that can simulate flow and advection conditionally}
|
||||
\caption{Main loop that can simulate velocities and advection conditionally}
|
||||
\label{alg:stream4v1}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -150,7 +150,7 @@ do not consider the data that occurs so infrequently that it means nothing. We d
|
|||
\begin{algorithm}
|
||||
$u \leftarrow \texttt{Smooth}(u, smooth_u)$ \;
|
||||
$v \leftarrow \texttt{Smooth}(v, smooth_v)$ \;
|
||||
$w \leftarrow \texttt{Smooth}(w, smooth_w)$ \;
|
||||
$w \leftarrow \texttt{Smooth}(w, smooth_w, smooth_{vert})$ \;
|
||||
\caption{Smoothing the velocity}
|
||||
\label{alg:smoothv}
|
||||
\end{algorithm}
|
||||
|
|
@ -11,23 +11,23 @@ body, i.e. it absorbs all radiation on all wavelengths. This is captured in Stef
|
|||
In \autoref{eq:stefan-boltzmann} the symbols are:
|
||||
|
||||
\begin{itemize}
|
||||
\item $S$: The energy that reaches the top of the atmosphere, coming from the sun or a similar star, per second per meters squared $Wm^{-2}$. This is also called the insolation.
|
||||
\item $\sigma$: The Stefan-Boltzmann constant, $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}.
|
||||
\item $T$: The temperature of the planet ($K$).
|
||||
\item $S$: The energy that reaches the top of the atmosphere, coming from the sun or a similar star, per second per meters squared (\si{Wm^{-2}}). This is also called the insolation.
|
||||
\item $\sigma$: The Stefan-Boltzmann constant, $5.670373 \cdot 10^-8$ (\si{Wm^{-2}K^{-4}}) \cite{stefan-boltzmann}.
|
||||
\item $T$: The temperature of the planet (\si{K}).
|
||||
\end{itemize}
|
||||
|
||||
The energy difference between the energy that reaches the atmosphere and the temperature of the planet must be equal to the change in temperature of the planet, which is written in
|
||||
\autoref{eq:sb rewritten}. The symbols on the right hand side are:
|
||||
|
||||
\begin{itemize}
|
||||
\item $\Delta U$: The change of internal energy ($J$) \cite{thermo1}.
|
||||
\item $C$: The specific heat capacity of the object, i.e. how much energy is required to heat the object by one degree Kelvin ($\frac{J}{K}$).
|
||||
\item $\Delta T$: The change in temperature ($K$).
|
||||
\item $\Delta U$: The change of internal energy (\si{J}) \cite{thermo1}.
|
||||
\item $C$: The specific heat capacity of the object, i.e. how much energy is required to heat the object by one degree Kelvin (\si{JK^{-1}}).
|
||||
\item $\Delta T$: The change in temperature (\si{K}).
|
||||
\end{itemize}
|
||||
|
||||
We want to know the change of temperature $\Delta T$, so we rewrite the equation into \autoref{eq:sb rewritten2}. Here we added the $\delta t$ term to account for the time difference (or time
|
||||
step). This is needed as we need an interval to calculate the difference in temperature over. Also we need to get the energy that we get ($J$) and not the energy per second ($W$), and by adding
|
||||
this time step the units all match up perfectly.
|
||||
step). This is needed as we need an interval to calculate the difference in temperature over. Also we need to get the energy that we get (\si{J}) and not the energy per second (\si{W}), and by
|
||||
adding this time step the units all match up perfectly.
|
||||
|
||||
\begin{subequations}
|
||||
\begin{equation}
|
||||
|
|
@ -146,13 +146,13 @@ we change that) amd the current time is obviously not constant. All constants ca
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{time $t$, amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
\For{$lat \in [-90, 90]$}{
|
||||
\For{$lon \in [0, 360]$}{
|
||||
\For{$lat \leftarrow -90$ \KwTo $90$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $360$}{
|
||||
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
|
||||
$T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
|
||||
}
|
||||
}
|
||||
\caption{The main function of the temperature calculations}
|
||||
\caption{The main function for the temperature calculations}
|
||||
\label{alg:stream1v1}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -197,13 +197,13 @@ that $S$ is defined as the call to \autoref{alg:solar}.
|
|||
\Input{amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lon \in [0, nlot]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlot$}{
|
||||
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
|
||||
$T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
|
||||
}
|
||||
}
|
||||
\caption{The main function of the temperature calculations}
|
||||
\caption{The main function for the temperature calculations}
|
||||
\label{alg:stream1v2}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -229,13 +229,13 @@ reflected instead of absorbed, where we need the amount that is absorbed which i
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lon \in [0, nlot]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlot$}{
|
||||
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t ((1 - a[lat, lon])S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p[lat, lon]}$ \;
|
||||
$T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
|
||||
}
|
||||
}
|
||||
\caption{The main loop of the temperature calculations}
|
||||
\caption{The main function for the temperature calculations}
|
||||
\label{alg:stream2v3}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -249,13 +249,13 @@ The air density is not at all points exactly the same. This may be due to the wi
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lon \in [0, nlot]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlot$}{
|
||||
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t ((1 - a[lat, lon])S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{\rho[lat, lon]C_p[lat, lon]}$ \;
|
||||
$T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{\rho[lat, lon]C_a}$ \;
|
||||
}
|
||||
}
|
||||
\caption{The main loop of the temperature calculations}
|
||||
\caption{The main function for the temperature calculations}
|
||||
\label{alg:temperature with density}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -294,15 +294,15 @@ regards to the storage of the values) by the addition of multiple atmospheric la
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lon \in [0, nlot]$}{
|
||||
\For{$layer \in [0, nlevels]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlot$}{
|
||||
\For{$layer \leftarrow 0$ \KwTo $nlevels$}{
|
||||
$T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t ((1 - a[lat, lon])S + \sigma(4\epsilon[0](T_a[lat, lon, 0])^4 - 4(T_p[lat, lon])^4))}
|
||||
{4C_p[lat, lon]}$ \;
|
||||
\uIf{$layer == 0$}{
|
||||
\uIf{$layer = 0$}{
|
||||
$T_a[lat, lon, layer] \leftarrow T_a[lat, lon, layer] + \frac{\delta t \sigma((T_p[lat, lon])^4 - 2\epsilon[layer](T_a[lat, lon, layer])^4)}
|
||||
{\rho[lat, lon, layer]C_a\delta z[layer]}$ \;
|
||||
}\uElseIf{$layer == nlevels - 1$}{
|
||||
}\uElseIf{$layer = nlevels - 1$}{
|
||||
$T_a[lat, lon, layer] \leftarrow T_a[lat, lon, layer] + \frac{\delta t \sigma(\epsilon[layer - 1](T_a[lat, lon, layer - 1])^4 - 2\epsilon[layer](T_a[lat, lon, layer])^4)}
|
||||
{\rho[lat, lon, layer]C_a\delta z[layer]}$ \;
|
||||
}\uElse{
|
||||
|
|
@ -312,7 +312,7 @@ regards to the storage of the values) by the addition of multiple atmospheric la
|
|||
}
|
||||
}
|
||||
}
|
||||
\caption{The main loop of the temperature calculations}
|
||||
\caption{The main function for the temperature calculations}
|
||||
\label{alg:temperature layer}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -320,7 +320,7 @@ We also need to initialise the $\epsilon$ value for each layer. We do that in \a
|
|||
|
||||
\begin{algorithm}
|
||||
$\epsilon[0] \leftarrow 0.75$ \;
|
||||
\For{$i \in [1, nlevels]$}{
|
||||
\For{$i \leftarrow 1$ \KwTo $nlevels$}{
|
||||
$\epsilon[i] \leftarrow 0.5\epsilon[i - 1]$
|
||||
}
|
||||
\caption{Intialisation of the insulation of each layer (also known as $\epsilon$)}
|
||||
|
|
@ -379,8 +379,8 @@ is described in \autoref{eq:optical depth}. The symbols in the equation mean:
|
|||
|
||||
\begin{itemize}
|
||||
\item $\tau_0$: Optical depth at the surface.
|
||||
\item $p$: Atmospheric pressure ($Pa$).
|
||||
\item $p_s$: Atmospheric pressure at the surface ($Pa$).
|
||||
\item $p$: Atmospheric pressure (\si{Pa}).
|
||||
\item $p_s$: Atmospheric pressure at the surface (\si{Pa}).
|
||||
\item $f_l$: The linear optical depth parameter, with a value of 0.1.
|
||||
\end{itemize}
|
||||
|
||||
|
|
@ -406,23 +406,23 @@ $S_z$ represents the call to \autoref{alg:gradient z}. \texttt{solar} represents
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{amount of energy that hits the planet $S$}
|
||||
\Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lon \in [0, nlon]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlot$}{
|
||||
$pressureProfile \leftarrow p[lat, lon, :]$ \;
|
||||
$\tau = \tau_0(lat)f_l\frac{pressureProfile}{pressureProfile[0]} + (1 - f_l)(\frac{pressureProfile}{pressureProfile[0]})^4)$ \;
|
||||
$U[0] \leftarrow \sigma T_p[lat, lon]^4$ \;
|
||||
|
||||
\For{$level \in [1, nlevels]$}{
|
||||
\For{$level \leftarrow 1$ \KwTo $nlevels$}{
|
||||
$U[level] \leftarrow U[level - 1] - \frac{(\tau[level] - \tau[level - 1])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level - 1] - \tau[level])}$ \;
|
||||
}
|
||||
|
||||
$D[nlevels - 1] \leftarrow 0$ \;
|
||||
|
||||
\For{$level \in [nlevels - 1, 0]$}{
|
||||
\For{$level \leftarrow nleves - 1$ \KwTo $0$}{
|
||||
$D[level] \leftarrow D[level + 1] - \frac{(\tau[level + 1] - \tau[level])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level] - \tau[level + 1])}$ \;
|
||||
}
|
||||
|
||||
\For{$level \in [0, nlevels]$}{
|
||||
\For{$level \leftarrow 0$ \KwTo $nlevels$}{
|
||||
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
|
||||
}
|
||||
|
||||
|
|
@ -433,7 +433,7 @@ $S_z$ represents the call to \autoref{alg:gradient z}. \texttt{solar} represents
|
|||
$T_p[lat, lon] \leftarrow T_p[lat, lon] \frac{\delta t((1 - a[lat, lon]) S + S_z(D, 0, 0, 0) - \sigma T_p[lat, lon]^4)}{C_p[lat ,lon]}$ \;
|
||||
}
|
||||
}
|
||||
\caption{Adding in radiation}
|
||||
\caption{Main function for calculating the temperature using radiation}
|
||||
\label{alg:optical depth}
|
||||
\end{algorithm}
|
||||
|
||||
|
|
@ -450,7 +450,7 @@ $Q$. We add in a check to see if we are currently calculating the radiation in t
|
|||
height.
|
||||
|
||||
\begin{algorithm}
|
||||
\For{$level \in [0, nlevels]$}{
|
||||
\For{$level \leftarrow 0$ \KwTo $nlevels$}{
|
||||
$Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
|
||||
\uIf{$heights[level] > 20 \cdot 10^3$}{
|
||||
$Q[level] \leftarrow Q[level] + \texttt{solar}(5, lat, lon, t) \frac{24 \cdot 60 \cdot 60(\frac{heights[level] - 20 \cdot 10^3}{10^3})^2}{30^2}$ \;
|
||||
|
|
|
|||
|
|
@ -23,7 +23,7 @@ calculations we will do later on.
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
|
||||
\Output{Gradient in the $x$ direction}
|
||||
\eIf{$k == \texttt{NULL}$}{
|
||||
\eIf{$k = \texttt{NULL}$}{
|
||||
$grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
|
||||
}{
|
||||
$grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon, k] - a[i, (j - 1) \text{ mod } nlon, k]}{\delta x[i]}$ \;
|
||||
|
|
@ -38,18 +38,18 @@ calculations we will do later on.
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
|
||||
\Output{Gradient in the $y$ direction}
|
||||
\eIf{$k == \texttt{NULL}$}{
|
||||
\eIf{$k = \texttt{NULL}$}{
|
||||
\uIf{$i == 0$}{
|
||||
$grad \leftarrow 2 \frac{a[i + 1, j] - a[i, j]}{\delta y}$ \;
|
||||
}\uElseIf{$i == nlat - 1$}{
|
||||
}\uElseIf{$i = nlat - 1$}{
|
||||
$grad \leftarrow 2 \frac{a[i, j] - a[i - 1, j]}{\delta y}$ \;
|
||||
}\uElse{
|
||||
$grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
|
||||
}
|
||||
}{
|
||||
\uIf{$i == 0$}{
|
||||
\uIf{$i = 0$}{
|
||||
$grad \leftarrow 2 \frac{a[i + 1, j, k] - a[i, j, k]}{\delta y}$ \;
|
||||
}\uElseIf{$i == nlat - 1$}{
|
||||
}\uElseIf{$i = nlat - 1$}{
|
||||
$grad \leftarrow 2 \frac{a[i, j, k] - a[i - 1, j, k]}{\delta y}$ \;
|
||||
}\uElse{
|
||||
$grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
|
||||
|
|
@ -65,18 +65,18 @@ calculations we will do later on.
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$}
|
||||
\Output{Gradient in the $z$ direction}
|
||||
\uIf{$a.dimensions == 1$}{
|
||||
\uIf{$k == 0$}{
|
||||
\uIf{$a.dimensions = 1$}{
|
||||
\uIf{$k = 0$}{
|
||||
$grad \leftarrow \frac{a[k + 1] - a[k]}{\delta z[k]}$ \;
|
||||
}\uElseIf{$k == nlevels - 1$}{
|
||||
}\uElseIf{$k = nlevels - 1$}{
|
||||
$grad \leftarrow \frac{a[k] - a[k - 1]}{\delta z[k]}$ \;
|
||||
}\uElse{
|
||||
$grad \leftarrow \frac{a[k + 1] - a[k - 1]}{2\delta z[k]}$ \;
|
||||
}
|
||||
} \uElse {
|
||||
\uIf{$k == 0$}{
|
||||
\uIf{$k = 0$}{
|
||||
$grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
|
||||
}\uElseIf{$k == nlevels - 1$}{
|
||||
}\uElseIf{$k = nlevels - 1$}{
|
||||
$grad \leftarrow \frac{a[i, j, k] - a[i, j, k - 1]}{\delta z[k]}$ \;
|
||||
}\uElse{
|
||||
$grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k - 1]}{2\delta z[k]}$ \;
|
||||
|
|
@ -131,17 +131,17 @@ respectively.
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{A matrix (double array) a}
|
||||
\Output{A matrix (double array) with results for the laplacian operator for each element}
|
||||
\eIf{$a.dimensions == 2$}{
|
||||
\For{$lat \in [1, nlat - 1]$}{
|
||||
\For{$lon \in [0, nlon]$}{
|
||||
\eIf{$a.dimensions = 2$}{
|
||||
\For{$lat \leftarrow 1$ \KwTo $nlat - 1$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlon$}{
|
||||
$output[lat, lon] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon)}{\delta x[lat]} + \frac{\Delta_y(a, lat + 1, lon) -
|
||||
\Delta_y(a, lat - 1, lon)}{\delta y}$\;
|
||||
}
|
||||
}
|
||||
}{
|
||||
\For{$lat \in [1, nlat - 1]$}{
|
||||
\For{$lon \in [0, nlon]$}{
|
||||
\For{$k \in [0, nlevels - 1]$}{
|
||||
\For{$lat \leftarrow 1$ \KwTo $nlat - 1$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlon$}{
|
||||
\For{$k \leftarrow 0$ \KwTo $nlevels - 1$}{
|
||||
$output[lat, lon, k] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon, k) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon, k)}{\delta x[lat]} + \frac{\Delta_y(a,
|
||||
lat + 1, lon, k) - \Delta_y(a, lat - 1, lon, k)}{\delta y} + \frac{\Delta_z(a, lat, lon, k + 1) - \Delta_z(a, lat, lon, k + 1)}{2\delta z[k]}$\;
|
||||
}
|
||||
|
|
@ -164,12 +164,9 @@ as we expect that we might use it in combination with the divergence operator mo
|
|||
\SetKwInOut{Output}{Output}
|
||||
\Input{A matrix (double array) $a$}
|
||||
\Output{A matrix (double array) containing the result of the divergence operator taken over that element}
|
||||
$dim_1 \leftarrow \text{ Length of } a \text{ in the first dimension}$ \;
|
||||
\For{$i \in [0, dim_1]$}{
|
||||
$dim_2 \leftarrow \text{ Length of } a \text{ in the second dimension (i.e. the length of the array stored at index } i)$ \;
|
||||
\For{$j \in [0, dim_2]$}{
|
||||
$dim_3 \leftarrow \text{ Length of } a \text{ in the third dimension}$ \;
|
||||
\For{$k \in [0, dim_3]$}{
|
||||
\For{$i \leftarrow 0$ \KwTo $a.length$}{
|
||||
\For{$j \leftarrow 0$ \KwTo $a[i].length$}{
|
||||
\For{$k \leftarrow 0$ \KwTo $a[i, j].length$}{
|
||||
$output[i, j, k] \leftarrow \Delta_x(au, i, j, k) + \Delta_y(av, i, j, k) + \Delta_z(aw, i, j, k)$ \;
|
||||
}
|
||||
}
|
||||
|
|
@ -245,41 +242,47 @@ Fourier Transform. The final equation is given in \autoref{eq:ft}.
|
|||
These Fourier Transforms have the great property that if you add them together you still have the relatively large distances of the center of mass to the origin at the original frequencies. It
|
||||
is this property that enables us to find the frequencies that compose a sound wave.
|
||||
|
||||
Before moving on we first need to discuss some important properties of FTs. Actually, these properties only apply to complex roots of unity (the $e^{2\pi ift}$ part is a complex root of unity).
|
||||
Due to the FT only being a coefficient $g(t)$ in front of a complex root of unity all these properties apply to FTs as well. Let us first start with the cancellation lemma as described in
|
||||
\autoref{lemma:cancellation} \cite{cancellation}. Note that $\omega^k_n = e^{\frac{2\pi k}{n}}$ which is very similar to the form we discussed previously if you realise that $f = \frac{1}{T}$.
|
||||
So replace $k$ by $t$ and $n$ by $T = \frac{1}{f}$ and you get the form we have discussed earlier.
|
||||
|
||||
\begin{lemma}[Cancellation Lemma] \label{lemma:cancellation}
|
||||
$\omega^{dk}_{dn} = \omega^k_n$, for all $k \geq 0, n \geq 0$ and $d > 0$.
|
||||
\end{lemma}
|
||||
|
||||
With that we also need to talk about the halving lemma (\autoref{lemma:halving}). This means that $(\omega^{k + \frac{n}{2}}_n)^2 = (\omega^k_n)^2$.
|
||||
|
||||
\begin{lemma}[Halving Lemma] \label{lemma:halving}
|
||||
If $n > 0$ is even, then the squares of the $n$-complex $n^{\text{th}}$ roots of unity are the $\frac{n}{2}$-complex $\frac{n}{2}^{\text{th}}$ roots of unity.
|
||||
\end{lemma}
|
||||
|
||||
Now that we know what a Fourier Transform is, we need to make it a Fast Fourier Transform, as you can imagine that calculating such a thing is quite difficult. Some smart people have thought
|
||||
about this and they came up with quite a fast algorithm, the Cooley-Tukey algorithm \cite{fft}, named after the people that thought of it. They use something we know as a Discrete Fourier
|
||||
Transform which is described by \autoref{eq:dft}. Here $N$ is the total amount of samples from the continuous sound wave. This means that $0 \leq k \leq N - 1$ and $0 \leq n \leq N - 1$.
|
||||
|
||||
Now with the DFT out of the way we can discuss the algorithm. It makes use of a clever property of the DFT. If you replace $k$ by $(N + k)$, as in \autoref{eq:dftex}, you can split up the
|
||||
exponent which will transform one of the two parts into $1$ due to $e^{-i2\pi n} = 1$ for any integer $n$. This means that $X_{N + k} = X_k \Rightarrow X_{k + iN} = X_k$ for any integer $i$.
|
||||
This symmetry as it is called can be exploited to produce a divide and conquer algorithm, which will recursivley calculate the FT which gives a running time of $N\log(n)$ as shown in
|
||||
\autoref{alg:FFT}.
|
||||
|
||||
\begin{subequations}
|
||||
\begin{equation}
|
||||
X_k = \sum_{n = 0}^{N - 1} x_ne^{-\frac{i2\pi}{N}kn}
|
||||
\label{eq:dft}
|
||||
\end{equation}
|
||||
\begin{equation}
|
||||
X_{N + k} = \sum_{n = 0}^{N - 1} x_ne^{-\frac{i2\pi(N + k)}{N}n} = \sum_{n = 0}^{N - 1} x_ne^{-i2\pi n}x_ne^{-\frac{i2\pi}{N}kn} = \sum_{n = 0}^{N - 1} x_ne^{-\frac{i2\pi}{N}kn}
|
||||
\label{eq:dftex}
|
||||
\end{equation}
|
||||
\end{subequations}
|
||||
Now with the DFT out of the way we can discuss the algorithm. Before we do, we assume that we have an input of $n$ elements, where $n$ is an exact power of $2$. Meaning $n = 2^k$ for some
|
||||
integer $k$.
|
||||
|
||||
\begin{algorithm}
|
||||
\SetKwInOut{Input}{Input}
|
||||
\SetKwInOut{Output}{Output}
|
||||
\Input{array $A$, integer $i$, integer $N$, integer $s$}
|
||||
\Input{array $A$}
|
||||
\Output{array $B$ with length $A.length - 1$ containing the DFT}
|
||||
\uIf{$N = 1$}{
|
||||
$B[0] \leftarrow A[0]$ \;
|
||||
} \uElse{
|
||||
$B[0], \dots, B[\frac{N}{2} - 1] \leftarrow \texttt{FFT}(A, i, \frac{N}{2}, 2s)$ \;
|
||||
$B[\frac{N}{2}], \dots, B[N - 1] \leftarrow \texttt{FFT}(A, i + s, \frac{N}{2}, 2s)$ \;
|
||||
\For{$k = 0$ to $\frac{N}{2} - 1$}{
|
||||
$t \leftarrow B[k]$ \;
|
||||
$B[k] \leftarrow t + e^{-2i\pi \frac{k}{N}} B[k + \frac{N}{2}]$ \;
|
||||
$B[k + \frac{N}{2}] \leftarrow t - e^{-2i\pi \frac{k}{N}} B[k + \frac{N}{2}]$ \;
|
||||
}
|
||||
$n \leftarrow A.length$ \;
|
||||
\uIf{$n = 1$}{
|
||||
\Return{$A$} \;
|
||||
}
|
||||
$\omega_n \leftarrow e^{\frac{2\pi i}{n}}$ \;
|
||||
$\omega \leftarrow 1$ \;
|
||||
$a^0 \leftarrow (A[0], A[2], \dots, A[n - 2])$ \;
|
||||
$a^1 \leftarrow (A[1], A[3], \dots, A[n - 1])$ \;
|
||||
$y^0 \leftarrow $ \texttt{FFT}($a^0$) \;
|
||||
$y^1 \leftarrow $ \texttt{FFT}($a^1$) \;
|
||||
\For{$k \leftarrow 0$ \KwTo $\frac{n}{2} - 1$}{
|
||||
$B[k] \leftarrow y^0[k] + \omega y^1[k]$ \;
|
||||
$B[k + \frac{n}{2}] \leftarrow y^0[k] - \omega y^1[k]$ \;
|
||||
$\omega \leftarrow \omega \omega_n$ \;
|
||||
}
|
||||
\Return{B}
|
||||
\caption{One dimensional Fast Fourier Transformation}
|
||||
|
|
@ -296,25 +299,31 @@ second dimension. This can of course be extended in the same way for a $p$-dimen
|
|||
\end{equation}
|
||||
|
||||
It is at this point that the algorithm becomes very complicated. Therefore I would like to invite you to use a library for these kinds of calculations, like Numpy \cite{numpy} for Python. If you
|
||||
really want to use your own made version, you need to wait for a bit as I try to decode the literature on the algorithm for it. It is one heck of a thing to decode, so I decided to treat that at
|
||||
a later point in the future. It will then be found here, so you will have to stay tuned. Sorry about that, it's quite complex...
|
||||
really want to use your own made version, or if you want to understand how the libraries sort of do it, then you may continue. I still advise you to use a library, as other people have made the
|
||||
code for you.
|
||||
|
||||
With that out of the way (or rather on the TODO list), we need to create a smoothing operation out of it. We do this in \autoref{alg:smooth}. Keep in mind that \texttt{FFT} the call is to the
|
||||
We can calculate the result of such a multidimensional FFT by computing one-dimensional FFTs along each dimension in turn. Meaning we first calculate the result of the FFT along one dimension,
|
||||
then we do that for the second dimension and so forth. Now the order of calculating the FFTs along each dimension does not matter. So you can calculate the third dimension first, and the first
|
||||
dimension after that if you wish. This has the advantage that we can program a generic function without keeping track of the order of the dimensions. Now, the code is quite complicated but it
|
||||
boils down to this: calculate the FFT along one dimension, then repeat it for the second dimension and multiply the two together and keep repeating that for all dimensions.
|
||||
|
||||
With that out of the way, we need to create a smoothing operation out of it. We do this in \autoref{alg:smooth}. Keep in mind that \texttt{FFT} the call is to the
|
||||
multidimensional Fast Fourier Transform algorithm, \texttt{IFFT} the call to the inverse of the multidimensional Fast Fourier Transform algorithm (also on the TODO list) and that the $int()$
|
||||
function ensures that the number in brackets is an integer. Also note that the inverse of the FFT might give complex answers, and we only want real answers which the $.real$ ensures. We only
|
||||
take the real part and return that.
|
||||
take the real part and return that. $v$ is an optional parameter with default value $0.5$ which does absolutely nothing in this algorithm.
|
||||
|
||||
\begin{algorithm}
|
||||
\SetKwInOut{Input}{Input}
|
||||
\SetKwInOut{Output}{Output}
|
||||
\Input{Array $a$, smoothing factor $s$}
|
||||
\Input{Array $a$, smoothing factor $s$, vertical smoothing factor $v \leftarrow 0.5$}
|
||||
\Output{Array $A$ with less variation}
|
||||
$nlat \leftarrow a.length$ \;
|
||||
$nlon \leftarrow a[0].length$ \;
|
||||
$nlevels \leftarrow a[0][0].length$ \;
|
||||
$temp \leftarrow \texttt{FFT}(a)$ \;
|
||||
$temp[int(nlat s):int(nlat(1 - s)),:,:] \leftarrow 0$ \;
|
||||
$temp[:,int(nlon s):int(nlon(1 - s)),:] \leftarrow 0$ \;
|
||||
$temp[:,int(nlon s):int(nlon(1 - s)),:] \leftarrow 0$ \;
|
||||
$temp[:,:,int(nlevels v):int(nlevels(1 - v))] \leftarrow 0$ \;
|
||||
\Return $\texttt{IFFT}(temp).real$ \;
|
||||
\caption{Smoothing function}
|
||||
\label{alg:smooth}
|
||||
|
|
|
|||
|
|
@ -6,17 +6,17 @@ The equation of state relates one or more variables in a dynamical system (like
|
|||
described by \autoref{eq:ideal gas} \cite{idealGas}. The symbols in that equation represent:
|
||||
|
||||
\begin{itemize}
|
||||
\item $p$: The gas pressure ($Pa$).
|
||||
\item $V$: The volume of the gas ($m^3$).
|
||||
\item $n$: The amount of moles\footnote{Mole is the amount of particles ($6.02214076 \cdot 10^{23}$) in a substance, where the average weight of one mole of particles in grams is about the
|
||||
same as the weight of one particle in atomic mass units ($u$)\cite{mole}} in the gas.
|
||||
\item $R$: The Gas constant, $8.3144621$ ($J(mol)^{-1}K$) \cite{idealGas}.
|
||||
\item $p$: The gas pressure (\si{Pa}).
|
||||
\item $V$: The volume of the gas (\si{m^3}).
|
||||
\item $n$: The amount of moles in the gas (\si{mol}).
|
||||
\item $R$: The Gas constant as defined in \autoref{sec:gas constant} (\si{JK^{-1}mol^{-1}}) \cite{idealGas}.
|
||||
\item $T$: The temperature opf the gas ($K$).
|
||||
\end{itemize}
|
||||
|
||||
If we divide everything in \autoref{eq:ideal gas} by $V$ and set it to be unit (in this case, set it to be exactly $1 m^3$) we can add in the molar mass in both the top and bottom parts of the
|
||||
division as show in \autoref{eq:gas unit}. We can then replace $\frac{nm}{V}$ by $\rho$ the density of the gas ($kgm^{-3}$) and $\frac{R}{m}$ by $R_s$ the specific gas constant (gas constant that varies per
|
||||
gas in $J(mol)^{-1}K$) as shown in \autoref{eq:state gas}. the resulting equation is the equation of state that you get that most atmospheric physicists use when talking about the atmosphere \cite{simon}.
|
||||
If we divide everything in \autoref{eq:ideal gas} by $V$ and set it to be unit (in this case, set it to be exactly $1$ \si{m^3}) we can add in the molar mass in both the top and bottom parts of
|
||||
the division as show in \autoref{eq:gas unit}. We can then replace $\frac{nm}{V}$ by $\rho$ the density of the gas (\si{kgm^{-3}}) and $\frac{R}{m}$ by $R_s$ the specific gas constant (gas
|
||||
constant that varies per gas in \si{JK^{-1}mol^{-1}}) as shown in \autoref{eq:state gas}. The resulting equation is the equation of state that you get that most atmospheric physicists use when
|
||||
talking about the atmosphere \cite{simon}.
|
||||
|
||||
\begin{subequations}
|
||||
\begin{equation}
|
||||
|
|
@ -47,19 +47,19 @@ to the f-plane approximation are given in \autoref{eq:x momentum} and \autoref{e
|
|||
effect on the whole flow. All the symbols in \autoref{eq:x momentum} and \autoref{eq:y momentum} mean:
|
||||
|
||||
\begin{itemize}
|
||||
\item $u$: The east to west velocity ($ms^{-1}$).
|
||||
\item $t$: The time ($s$).
|
||||
\item $u$: The east to west velocity (\si{ms^{-1}}).
|
||||
\item $t$: The time (\si{s}).
|
||||
\item $f$: The coriolis parameter as in \autoref{eq:coriolis}.
|
||||
\item $v$: The north to south velocity ($ms^{-1}$).
|
||||
\item $\rho$: The density of the atmosphere ($kgm^{-3}$).
|
||||
\item $p$: The atmospheric pressure ($Pa$).
|
||||
\item $x$: The local longitude coordinate ($m$).
|
||||
\item $y$: The local latitude coordinate ($m$).
|
||||
\item $v$: The north to south velocity (\si{ms^{-1}}).
|
||||
\item $\rho$: The density of the atmosphere (\si{kgm^{-3}}).
|
||||
\item $p$: The atmospheric pressure (\si{Pa}).
|
||||
\item $x$: The local longitude coordinate (\si{m}).
|
||||
\item $y$: The local latitude coordinate (\si{m}).
|
||||
\end{itemize}
|
||||
|
||||
If we then define a vector $\bar{u}$ as $(u, v, 0)$, we can rewrite both \autoref{eq:x momentum} as \autoref{eq:x momentum laplace}. Here $\nabla u$ is the gradient of
|
||||
$u$ in both $x$ and $y$ directions. Then if we write out $\nabla u$ we get \autoref{eq:x momentum final}. Similarly, if we want to get $\delta v$ instead of $\delta u$ we rewrite
|
||||
\autoref{eq:y momentum} to get \autoref{eq:y momentum laplace} and \autoref{eq:y momentum final}.
|
||||
If we then define a vector $\bar{u}$ as $(u, v, 0)$, we can rewrite both \autoref{eq:x momentum} as \autoref{eq:x momentum laplace}. Here $\nabla u$ is the gradient of $u$ in both $x$ and $y$
|
||||
directions. Then if we write out $\nabla u$ we get \autoref{eq:x momentum final}. Similarly, if we want to get $\delta v$ instead of $\delta u$ we rewrite \autoref{eq:y momentum} to get
|
||||
\autoref{eq:y momentum laplace} and \autoref{eq:y momentum final}.
|
||||
|
||||
\begin{subequations}
|
||||
\begin{equation}
|
||||
|
|
@ -98,8 +98,8 @@ With the gradient functions defined in \autoref{alg:gradient x} and \autoref{alg
|
|||
$S_{yv} \leftarrow \texttt{gradient\_y}(v, lan, lon)$ \;
|
||||
$S_{px} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
|
||||
$S_{py} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
|
||||
\For{$lat \in [1, nlat - 1]$}{
|
||||
\For{$lon \in [0, nlon]$}{
|
||||
\For{$lat \leftarrow 1$ \KwTo $nlat - 1$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlon$}{
|
||||
$u[lan, lon] \leftarrow u[lan, lon] + \delta t \frac{-u[lan, lon]S_{xu} - v[lan, lon]S_{yu} + f[lan]v[lan, lon] - S_{px}}{\rho}$ \;
|
||||
$v[lan, lon] \leftarrow v[lan, lon] + \delta t \frac{-u[lan, lon]S_{xv} - v[lan, lon]S_{yv} - f[lan]u[lan, lon] - S_{py}}{\rho}$ \;
|
||||
}
|
||||
|
|
@ -116,7 +116,7 @@ Another change introduced is in the coriolis parameter. Up until now it has been
|
|||
\SetAlgoLined
|
||||
$\Omega \leftarrow 7.2921 \cdot 10^{-5}$ \;
|
||||
|
||||
\For{$lat \in [-nlat, nlat]$}{
|
||||
\For{$lat \leftarrow -nlat$ \KwTo $nlat$}{
|
||||
$f[lat] \leftarrow 2\Omega \sin(lat \frac{\pi}{180})$ \;
|
||||
}
|
||||
\caption{Calculating the coriolis force}
|
||||
|
|
@ -141,9 +141,9 @@ then so does $y$.}
|
|||
$S_{yv} \leftarrow \texttt{gradient\_y}(v, lan, lon)$ \;
|
||||
$S_{px} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
|
||||
$S_{py} \leftarrow \texttt{gradient\_y}(p, lan, lon)$ \;
|
||||
\For{$lat \in [1, nlat - 1]$}{
|
||||
\For{$lon \in [0, nlon]$}{
|
||||
\For{$layer \in [0, nlevels]$}{
|
||||
\For{$lat \leftarrow 1$ \KwTo $nlat - 1$}{
|
||||
\For{$lon \leftarrow 0$ \KwTo $nlon$}{
|
||||
\For{$layer \leftarrow 0$ \KwTo $nlevels$}{
|
||||
$u[lan, lon, layer] \leftarrow u[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xu} - v[lat, lon, layer]S_{yu} + f[lat]v[lat, lon, layer] - S_{px}}{\rho}$ \;
|
||||
$v[lan, lon, layer] \leftarrow v[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xv} - v[lat, lon, layer]S_{yv} - f[lat]u[lat, lon, layer] - S_{py}}{\rho}$ \;
|
||||
$w[lan, lon, layer] \leftarrow w[lat, lon, layer] - \frac{p[lat, lon, layer] - p_o[lat, lon, layer]}{\delta t\rho[lat, lon, layer]g}$ \;
|
||||
|
|
|
|||
Loading…
Reference in New Issue