Alpine3D

Alpine3D Svn Source Tree

Root/trunk/doc/manual/alpine3d_manual.tex

1\documentclass[12pt]{report}
2\usepackage{latexsym}
3\usepackage{times}
4\usepackage[T1]{fontenc}
5\usepackage{pstcol}
6\usepackage{bm}
7\usepackage{psfrag}
8\usepackage{epsfig}
9\usepackage{amsfonts,amssymb,amsmath}
10\usepackage{a4,isolatin1}
11\usepackage{listings} % Quellcode
12
13
14\lstset{language=[latex]tex}
15
16
17%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18% ------------------------------------------------------------
19% macros
20% ------------------------------------------------------------
21%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22\def\figuredir{./figures}
23\newcommand{\noi}{\noindent}
24
25
26
27
28
29
30%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31% ------------------------------------------------------------
32\title{PLEASE NOTE THAT THIS IS NOW OBSOLETE! DO NOT USE THIS MANUAL}
33% ------------------------------------------------------------
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35\author{
36 Alpine3D Contributers
37}
38\date{\today}
39
40
41
42
43
44
45
46
47
48\begin{document}
49\maketitle
50\tableofcontents
51%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53% ------------------------------------------------------------
54\chapter{Introduction}\label{ch:intro}
55% ------------------------------------------------------------
56%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58
59
60
61%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62% ------------------------------------------------------------
63\section{Overview}
64% ------------------------------------------------------------
65%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66Alpine3D is a spatially distributed, land-surface type model
67comprising several phyically and conceptually submodules: a snowpack
68model (SNOWPACK), a radiation balance model, a runoff model, and a
69snow transport model. Additionally, it provides some routines for
70interpolation of meteorological data. A brief overview of Alpine3D and
71its application to hydrological issues can be found in in
72\cite{lehning_2006}.
73
74Alpine3D is written in C (Snowpack), C++(energy balance, Snowdrift) and
75Fortran (runoff, interpolation) on Unix/Linux. It exploits the GRID
76technology POPC++ \cite{popc} for parallelizing individual submodules.
77To use the parallel capabilities of Alpine3D it must be installed
78together with POPC++. In a multiuser environment it is recommended to
79install additionally a GRID resource management system like the Globus
80Toolkit \cite{gt}.
81
82
83
84
85
86
87
88%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90% ------------------------------------------------------------
91\chapter{Installation}\label{ch:installation}
92% ------------------------------------------------------------
93%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95
96
97
98%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99% ------------------------------------------------------------
100\section{System Requirements}
101% ------------------------------------------------------------
102%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
103
104Alpine3D is a Linux application and in its sequential mode it can be
105installed on every (?) Unix/Linux system. Presently it can be used in
106its parallel mode only on the Linux cluster at SLF and on the HPC
107cluster ``Zeus'' at WSL.
108
109
110%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111% ------------------------------------------------------------
112\section{Source code}
113% ------------------------------------------------------------
114%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115
116The Alpine3D source code is maintained with SVN. For a good online
117documentation see \cite{svn}, and the basic usage of the \verb+svn+ is
118also briefly described in Ch.\ \ref{ch:svn}. The repository can be
119accessed (only locally at SLF) via the URL
120\verb+svn://svn.slf.local/alpine3d+ after receiving a password from
121the IT team. To check out a working copy of Alpine3D use the svn
122command
123\begin{verbatim}
124 svn co svn://svn.slf.local/alpine3d/alpine3d/trunk
125\end{verbatim}
126in a terminal. This will create a subdirectory \verb+trunk+ in your
127{\em present} working directory. In the \verb+trunk+ directory you will find
128the following files and subdirectories of Alpine3D:
129
130\verb+./main/+ (main application class)
131
132\verb+./common/+ (array classes, data marshalling functions from POPC++)
133
134\verb+./ebalance/+ (energy balance class, radiation model, view factors
135etc)
136
137\verb+./snowdrift/+ (snowdrift class: saltation, suspension)
138
139\verb+./snowdrift_par/+ (parallel snowdrift implementation)
140
141\verb+./snowpack/+ (snowpack class, core files of SNOWPACK)
142
143\verb+./inout/+ (i/o classes)
144
145\verb+./runoff/+ (runoff functions, FORTRAN!)
146
147\verb+./interpol/+ (interpolation routines for 2d meteo input, FORTRAN!)
148
149\verb+./deploy/+ (directory for popc objects )
150
151\verb+./Interface/+ (a gui for Alpine3D)
152
153\verb+./doc/+ (this manual, documentation)
154
155\verb+./tools/+ (some auxiliary scripts)
156
157\verb+./current_snow/+ (input directory, should be deleted from svn!)
158
159\verb+run_dischma_seq.sh+ (an example start script for sequential runs)
160
161\verb+run_dischma_par.sh+ (an example start script for parallel runs)
162
163\verb+Makefile+ (the top level Makefile)
164
165\verb+Makefile.par+ (the top level Makefile only for parallel drift)
166\\[1ex]
167
168\noindent Note that the program is continously improved and corrected.
169Make sure that you are up to date by checking the status of your
170working copy and update if necessary by typing
171\begin{verbatim}
172 svn up .
173\end{verbatim}
174in your \verb+trunk+-directory (see also Ch.\ \ref{ch:svn} for details)
175
176
177%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
178% ------------------------------------------------------------
179\section{Environment}
180% ------------------------------------------------------------
181%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182
183Using compilers, POPC++, and GT requires you to have your environment
184set up correctly. At SLF, add the following lines to an appropriate
185dot-file, i.e. to \verb+.profile+ or \verb+.bash_profile+ in your home
186directory.
187
188\begin{verbatim}
189 #-- globus --------------------------------------------------
190 export GLOBUS_LOCATION=/usr/local/gt3/gt3.2
191 source $GLOBUS_LOCATION/etc/globus-user-env.sh
192 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/gt3/gt3.2/lib
193
194 #-- popc ---------------------------------------------------
195 # Note : Var names changed from PAROC_ to POPC_ in POP-C++-1.3
196 export POPC_LOCATION=/usr/local/popc
197 export PATH=$PATH:/usr/local/popc/bin
198 export POPC_JOBSERVICE=grid1.slf.local:2711
199
200 #-- fortran -------------------------------------------------
201 export PATH=$PATH:/usr/local/intel_fc_81/bin
202 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/intel_fc_81/lib
203\end{verbatim}
204
205and source the file. Before complaining, make sure that everything is
206set correctly (i.e. verify your environment by typing \verb+set+ in a
207terminal).
208
209
210%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211% ------------------------------------------------------------
212\section{Compilation}
213% ------------------------------------------------------------
214%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215
216% ------------------------------------------------------------
217\subsection{Making at SLF}
218% ------------------------------------------------------------
219Compiling Alpine3D requires a C, C++ and a FORTRAN90 compiler. At SLF,
220gcc, g++ and ifort (Intel compiler) is used. Compilation is done using
221\verb+make+ with the top level Makefile:
222Use
223\begin{verbatim}
224make all
225\end{verbatim}
226in your trunk directory to generate the \verb+Alpine3D+ executable for
227sequential computation. In order to generate the executable
228\verb+Alpine3D.popc+ for parallel computing use
229\begin{verbatim}
230make all_par
231\end{verbatim}
232and subsequently
233\begin{verbatim}
234make deploy_par
235\end{verbatim}
236to move the popc-modules to the \verb+./deploy+ directory and create
237the object description file \verb+./deploy/objectmap.conf+ required by
238POPC++. For making the parallel snowdrift implementation see Ch.\
239\ref{ch:snowdrift}.
240
241
242% ------------------------------------------------------------
243\subsection{Specifics on Zeus}
244% ------------------------------------------------------------
245
246A brief summary of the HPC cluster Zeus at WSL is given below. You
247can login to zeus@wsl.ch after receiving a password from Markus
248Reinhard \\(\verb+markus.reinhardt@wsl.ch+)
249
250\begin{itemize}
251\item CPUs: AMD Opteron Dual Core 64 Bit processors. Each unit has two
252 Dual CPUs and the cluster comprises 16 nodes (ip-adresses
253 172.16.2.1-16) which makes 16 x 2 x 2 = 64 cores
254
255\item Memory:
256 Each unit has 8 GByte RAM.
257
258 \item Netork: The units are connnected by a Gigabyte Ethernet via
259 switch. Additionally, compute nodes are connected by Myrinet.
260 The transfer rate is ~900 MBytes/s (r+w).
261
262 \item OS:
263 Suse Linux 9.3
264
265 \item Software: Sun Grid Engine (sge), ls-dyna, globus, mpich, petsc,
266 pathscale (Compiler), R, IDL, ganglia, nagios
267\end{itemize}
268
269Moving the program to Zeus or any other computer requires to modify
270compilers and libraries in the top level makefile \verb+Makefile+. On
271Zeus, gcc, g++ and the pathscale Fortran compiler pathf90 can be used.
272Besides updating your environment accordingly the top level makefile
273\verb+Makefile+ has to be changed: Make sure that the respective
274linker flag
275\begin{verbatim}
276FLIBS=-L/usr/local/intel_fc_80/lib -limf -lifcore
277\end{verbatim}
278valid at SLF is replaced by the appropriate pathscale flag
279\begin{verbatim}
280FLIBS=-L/opt/pathscale/lib/2.2 -lmpath -lpathfortran
281\end{verbatim}
282In addition, make sure to add the flag \verb+-fno-second-underscore+
283to \verb+FFLAGS+. Then go on with making as described in the previous
284Section.
285
286Note: Sometimes nothing happens when Alpine3D.popc is started on
287Zeus. Try to open the \verb+./deploy/objectmap.conf+ file and remove
288the \verb+exports+ directory in the object description, i.e. replace
289\verb+/exports/home/...+ by \verb+/home/...+.
290
291It is also possible that the job manager unexpectedly dies on Zeus. In
292order to restart it, login as popc and run the command
293\verb+popc-1.0-inst/sbin/SXXpopc start+.
294
295%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
296%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
297% ------------------------------------------------------------
298\chapter{Running Alpine3D}\label{ch:running}
299% ------------------------------------------------------------
300%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
302
303
304
305%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
306% ------------------------------------------------------------
307\section{Testing Alpine3D}
308% ------------------------------------------------------------
309%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310
311After successful compilation, you can immediately test the sequential
312and the parallel version. By simply typing \verb+Alpine3D+ it will
313print a usage-message and you will get a list of possible command line
314switches with explanations.
315
316An example call to Alpine3D with appropriately set command line
317switches can be found in the example shell scripts
318\verb+run_dischma_*.sh+ in your \verb+trunk+ directory. These scripts start
319examples start your local version of Alpine3D with input data for the
320Dischma domain in the
321directory\\
322{\small \verb+/usr/local/org/FE/SCHNEEP/SMM/projects/alpine3d-testcases/dischma+}
323The input data sets up a simulation of the Dischma catchment with 2d
324meteorological input. Unfortunately, this directory is presently only accessible for the
325team ``Snowcover and Micrometeorology''.
326
327The output and the logfile \verb+stdouterr.log+ is written to\\
328{\small
329 \verb+/usr/local/org/FE/SCHNEEP/SMM/projects/alpine3d-testcases/dischma/output/test+
330}
331Output fields are compatible with the Ascii-ArcInfo grid format. All files
332follow the naming convention \verb+<julian-day>.<extension>+ where
333extensions are used as given in Tab.~\ref{tab1}
334
335\begin{table}[b]
336 \label{tab1}
337 \begin{center}
338 \begin{tabular}{cl}
339 .lwr & long wave radiation\\
340 .swr & short wave radiation\\
341 .sdp & snowdepth\\
342 .alb & albedo\\
343 .swe & snow water equivalent\\
344 .tss & snow surface temperature\\
345 .sca & snow covered area\\
346 \end{tabular}
347 \caption{File name extension in Alpine3D}
348 \end{center}
349\end{table}
350You can open these output files, e.g. with the Alpine3D GUI by switching to
351\verb+Interface+ subdirectory, typing \verb+view.bat+ and opening one of the files in that directory.
352
353
354
355%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356% ------------------------------------------------------------
357\section{Setting up a simulation}
358% ------------------------------------------------------------
359%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360Setting up your own simulation can be done by preparing input data and
361adjusting the parameters and paths in the example start scripts.
362Note, presently there are two (three) different files where parameters and paths have to
363be adjusted if you want to do simulations with 1d (2d) meteorological input. These are
364\begin{itemize}
365\item The start script.
366\item The snowpack parameter file specified by the -snowparam switch in the
367 start script
368\item The runoff/interpolation parameter file specified by the -meteopath switch in the
369 start script (only for 2d meteo input)
370\end{itemize}
371Make sure to have all paths in all files set correctly.
372
373
374%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
375% ------------------------------------------------------------
376\chapter{Input Data}\label{ch:input}
377% ------------------------------------------------------------
378%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
379
380%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
381% ------------------------------------------------------------
382\section{General}
383% ------------------------------------------------------------
384%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
385To be filled: Directories, formatting, time stamps, etc
386
387%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
388% ------------------------------------------------------------
389\section{Meteorological input}
390% ------------------------------------------------------------
391%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
392To be filled. Also ask for \cite{jon07}.
393Even in case 2d meteorological input is applied, some attention has to be paid concerning the remaining 1d meteorological input. Therefore, better read the 1d meteorological input section too.
394\subsection{1d meteorological input}
395\begin{itemize}
396\item No time gaps are possible; make sure the gaps are eliminated by extrapolation before starting Alpine3d since the sun position in the radiation balance module is always computed from the expected date following the main imposed Alpine3D time step.
397\item Starting time and time steps should be consistent between 1d- and 2d meteorological input.
398\item The radiation input (global shortwave and longwave radiation or cloud cover fraction) \textbf{has} to be measured at an \textbf{exposed} measurement site.
399\item Latitude, Longitude in the header of 1d meteorological file should be the coordinates of the center of the model domain; Swiss Coordinates and height in the header of 1d meteorological file should be the coordinates and height of the radiation input station.
400\item All meteorological input has to be in winter time (e.g.Davos: UTC+1) (i.e. no summer time is accounted for).
401\end{itemize}
402\subsection{2d meteorological input}
403%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
404%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
405% ------------------------------------------------------------
406\chapter{ToDo}\label{ch:todo}
407% ------------------------------------------------------------
408%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
409%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
410
411
412%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
413% ------------------------------------------------------------
414\section{Model problems}
415% ------------------------------------------------------------
416%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
417\begin{itemize}
418\item The model overestimates snowdepth in low elevations and
419 underestimates it in higher elevations.
420\end{itemize}
421
422%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
423% ------------------------------------------------------------
424\section{Implementation problems}
425% ------------------------------------------------------------
426%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427
428\begin{itemize}
429\item Reorganization (i.e. simplification) of the specification of
430 parameters and paths (e.g.: each module should come up with its own
431 parameter class \verb+<module>Param.cc/h/ph+ and paths can be set
432 only in the start script)
433\item Improve portability (i.e. GNUization with
434 the configure/autoconf/automake mechanism)
435\item Some binary output is generated by one of the Fortran codes and
436 written to stdout or stderr and hence to the logfile. This prevents
437 to search the logfile with grep. ->Incovenient.
438\end{itemize}
439
440
441
442%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
443%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
444% ------------------------------------------------------------
445\chapter{Subversion in a nutshell}\label{ch:svn}
446% ------------------------------------------------------------
447%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
448%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
449
450%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451% ------------------------------------------------------------
452\section{RTFM}
453% ------------------------------------------------------------
454%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
455Documentation can be found on the website
456
457\verb+http://subversion.tigris.org/+
458
459\noi and in the online book (sectioned html)
460
461\verb+http://svnbook.red-bean.com/+
462
463\noi which nicely tells you what it is all about:
464
465\noi {\em If C gives you enough rope to hang yourself, think of subversion as a
466sort of rope storage facility} (from the preface)
467
468%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
469% ------------------------------------------------------------
470\section{ Before you start }
471% ------------------------------------------------------------
472%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
473\begin{itemize}
474\item The repository at SLF has the (url) location
475
476 \verb+svn://svn.slf.local/alpine3d/alpine3d/trunk+.
477
478 \noi If you prefer access to the svn via command line it might be
479 convenient to to set
480
481 \verb+export SVN_ROOT="svn://svn.slf.local/alpine3d/alpine3d/trunk"+
482
483 There are also nice and more convenient x-applications for accessing
484 the repository. Check it out.
485
486\item If you even fail at {\em closing} old-school vi without the reference
487 card (like me) better
488
489 \verb+export EDITOR=emacs+
490
491 in your shell before committing something via command line.
492\end{itemize}
493
494%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
495% ------------------------------------------------------------
496\section{ Basic work cycle}
497% ------------------------------------------------------------
498%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
499\begin{itemize}
500\item Initial checkout:
501
502 \verb+svn checkout $SVN_ROOT+ (alias: \verb+svn co $SVN_ROOT+)
503
504 \noindent creates a working copy of the repository by copying the
505 complete directory tree (including the directory \verb+trunk+) from
506 the repository to your present working directory. Additionally
507 \verb+.svn+ subdirectories are created in each subdirectory which
508 hold status information about your working copy.
509\item Updating the working copy with the most recent version. Type
510
511 \verb+svn update+ (alias: \verb+svn up+)
512
513 \noindent in your \verb+trunk+ directory. If you want to update your
514 working copy with an older (less buggy?) version use
515
516 \verb+svn up --revision < revision no >+
517
518 \noindent You can also update only particular directories/files of
519 your working copy by using
520
521 \verb+svn up <filename>+
522
523\item Making changes to your working copy. Change your working copy by
524 by moving, copying, deleting or adding files. Use\\
525 \verb+svn add, svn delete, svn copy, svn move+ These changes take
526 place immediately in your working copy and the repository is changed
527 after your next \verb+commit+. Note: if you simply invoke system
528 commands, e.g. \verb+rm+, to delete a file, the deletion is not
529 scheduled as a future change in the repository and you might receive
530 an error at next \verb+commit+.
531
532\item Checking the status of the working copy:
533
534 \verb+svn status+
535
536 reports {\em if} there are changes in your files after your last
537 update. If you want to know {\em what} has changed use
538
539 \verb+svn diff+
540
541 \noindent If you want to suppress status information about your
542 messy working with bunch of additional files use the \verb+-q+
543 option. If you want to additionally have information if there are
544 more recent versions in the repository use
545
546 \verb+svn status -qu+
547
548 The first column of the output contains flags indicating the status
549 of the file, the most important is a capital ``C'': Conflict, then
550 your version won't compile. ``M'' means modified, ``A'' added and
551 ``D'' deleted. If a file has a ``*'' in the second column, a more
552 recent version of that file is available in the repository.
553
554\item Undo:
555
556 \verb+svn revert alpine3dFile.cc+
557
558 recovers the state of the file at previous update.
559\item Commit your changes:
560
561 \verb+svn commit+
562
563 \noindent If not invoked with the "-m" switch the default editor
564 will open for entering comments on your changes.
565
566\item Examining the revision history of the project:
567
568 \verb+svn log+
569
570 By passing a filename, the history of only this file is shown.
571 If you want to make sure to read the complete history of
572 alpine3d-changes then use
573
574 \verb+svn log svn://svn.slf.local/alpine3d/alpine3d/trunk+
575
576 All missing revision increments stem from changes in snowpack which
577 is maintained in the same repository in
578
579 \verb+svn://svn.slf.local/alpine3d/snowpack/+
580
581 If you really want to read {\em every commit message of the whole
582 repository}, type
583
584 \verb+svn log svn://svn.slf.local/alpine3d+
585
586\end{itemize}
587
588
589
590%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
591% ------------------------------------------------------------
592\section{Rules}
593% ------------------------------------------------------------
594%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
595
596Needless to say:
597\begin{enumerate}
598\item Don't commit without comment.
599\item Don't commit if it doesn't compile.
600\end{enumerate}
601
602
603
604
605
606
607%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
608%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
609% ------------------------------------------------------------
610\chapter{Snowdrift module}\label{ch:snowdrift}
611% ------------------------------------------------------------
612%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
613%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
614
615%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
616% ------------------------------------------------------------
617\section{Overview}
618% ------------------------------------------------------------
619%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
620The following description refers to the parallel version of the
621snowdrift module located in \verb+snowdrift_par+. From the
622computational point of view it is identical to the sequential version
623at the date where this directory was committed (20. December 2007).
624But in the parallel implementation in \verb+snowdrift_par+ all
625redundant routines have been removed and the naming of functions and
626variables has been made more intuitive.
627
628The snowdrift module comprises three different classes. The base class
629\verb+SnowDrift+ (Interface in \verb+SnowDrift.h+ and \verb+SnowDrift.ph+ and the derived classes \\
630\verb+SnowDriftParallel+ (Interface in \verb+SnowDriftParallel.h+
631\\and \verb+SnowDriftParallel.ph+) and \verb+SnowDriftWorker+
632(Interface also in \\ \verb+SnowDriftParallel.h+ and
633\verb+SnowDriftParallel.ph+) The implementation of the
634\verb+SnowDriftParallel+ class can be found exclusively in
635\verb+SnowDriftParallel.cc+ whereas the implementation of
636\verb+SnowDrift+ and \verb+SnowDriftWorker+ is distributed over the
637files
638\begin{verbatim}
639SnowDrift.cc
640Suspension.cc (contains solely the SnowDriftWorker::SubSuspension method)
641SnowDriftFEInit (contains initialization routines for the FE method)
642SnowDriftFENumerics (contains numerical element routines for the FE method)
643SnowDriftFEControl (contains additional FE routines for SnowDrift and
644SnowDriftWorker)
645\end{verbatim}
646
647The remaining files
648\begin{verbatim}
649PackSnowDrift.cc
650PackSnowDriftWorker.cc
651checksum.cc
652marshal_drift.cc
653\end{verbatim}
654are required for PopC++ and finally \verb+clock.h/cc+ contains a
655simple class for time measurements.
656
657For a basic understanding of the computation flow of the snowdrift
658model check the method \verb+SnowDrift::Compute()+ in the file
659\verb+SnowDrift.cc+ which consists basically of a call to the
660Saltation and the Suspension routine (besides some initialization and
661io-operations).
662
663For the Saltation you have to dig into the file \verb+Saltation.cc+.
664Some of the routines required there are still located in the
665\verb+snowpack_core+ subdirectory of the repository in the file
666\verb+./snowpack/snowpack_core/Saltation.c+
667
668For the suspension the main sequence of actions can be found in
669\\\verb+SnowDriftWorker::SubSuspension()+ in the file
670\verb+Suspension.cc+ with self-explenatory naming of the routines.
671
672
673%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
674% ------------------------------------------------------------
675\section{Details of the implementation of the Finite Element method}
676% ------------------------------------------------------------
677%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
678All theoretical preliminaries can be found in Marcs report
679\cite{ryser_2004}. As explained there, linear, hexahedral elements are
680used for the FE method. In the following I give further details about
681the location of individual procedures of the FE method.
682
683% ------------------------------------------------------------
684% ------------------------------------------------------------
685\subsection{Isoperimetric transformation}
686% ------------------------------------------------------------
687% ------------------------------------------------------------
688\begin{figure}[t]
689 \begin{center}
690 \noindent\includegraphics[angle=0,width=\textwidth]{\figuredir/elementMapping.eps}
691 \end{center}
692 \caption{Enumeration and coordinate systems in the reference and
693 original element respectively}
694 \label{iso}
695\end{figure}
696
697The nodes in the reference Element $\tilde{K}$ have coordinates
698\begin{align}
699 {P}^{(1)}&=(1,-1,-1)&\qquad{P}^{(5)}&=(1,-1,1) \\\nonumber
700 {P}^{(2)}&=(1, 1,-1)&\qquad{P}^{(6)}&=(1,1,1)\\\nonumber
701 {P}^{(3)}&=(-1,1,-1)&\qquad{P}^{(7)}&=(-1,1,1)\\\nonumber
702 {P}^{(4)}&=(-1,-1,-1)&\qquad{P}^{(8)}&=(-1,-1,1)\\\nonumber
703\end{align}
704Note that the coordinate system in the reference element is different
705from global coordinate system (cf. Fig.\ \ref{iso}. The 8-node
706trilinear hexahedron basisfunctions defined on the reference element
707$K$ are explicitely given by
708\begin{equation}
709 \phi^{(k)}(\mathbf{\xi}) =\frac{1}{8}\prod_{\alpha=1}^3
710 (1+\xi_{\alpha}P^{(k)}_{\alpha}),\qquad k=1\ldots 8
711\end{equation}
712The implementation can be found in \verb+computePhi()+ in
713\verb+SnowDriftFEnumerics.cc+.
714
715In terms of the $\phi^k$, the isoperimetric transformation
716$F:\tilde{K}\to K$ reads
717\begin{equation}\label{isoperimetric}
718 F_{\alpha}(\xi) = \sum_{i=1}^8 Q_{\alpha}^{(i)}\, \phi^{(i)}(\xi)
719\end{equation}
720and hence the Jacobian of the transformation is given by
721\begin{equation}
722 J_{\alpha,\beta}(\xi)
723 = \frac{\partial F_{\alpha}(\xi)}{\partial \xi_{\beta}}
724 = \sum_{k=1}^8 Q_{\alpha}^{(k)}\, \frac{\partial \phi^{(k)}(\xi)}{\partial \xi_{\beta}}
725\end{equation}
726which is implemented in \verb+computeJacobian+ in
727\verb+SnowDriftFEnumerics.cc+. The gradients of $\phi^k$ required for
728the jacobian are computed in \verb+computeGradPhi()+ in
729\verb+SnowDriftFEnumerics.cc+.
730
731
732% ------------------------------------------------------------
733% ------------------------------------------------------------
734\subsection{Integrals over elements}
735% ------------------------------------------------------------
736% ------------------------------------------------------------
737The integrals over elements are performed by integrating over
738reference elements. The weak solution of involves integrals over the
739original elmenents $K$. The integration is performed by transforming
740to the reference element $\tilde{K}$ via $F^{-1}$. Accordingly,
741integrals are computed over $\tilde{K}=F^{-1}(K)$ with the
742transformation formula
743\begin{equation}
744 \int_K \mathrm{d}^3 x\: f(x)
745 = \int_{\tilde{K}} \mathrm{d}^3 \xi\:
746 |\det( J(\xi) )|\,f(F(\xi))
747\end{equation}
748
749The integrals over reference elements are done by a two-point (per
750dimension) Gaussian quadrature with points $\pm 1/\sqrt{3}\approx
7510.57735$. These points are set in \verb+setQuadarturePoints+ in
752\verb+SnowDriftFEnumerics.cc+.
753
754Since some integrals over elements involve the inverse of the jacobian
755$J$ (see Eq.\ in \cite{ryser_2004}) it is advantageous it in terms of
756the adjugate of $J$ (denoted by $J_0$ in \cite{ryser_2004}. The
757inverse of $J$ and the adjugate $\mathrm{adj}(J)=$ are related by
758\begin{equation}
759 J^{-1}=\det(J)^{-1}\mathrm{adj}(J)
760\end{equation}
761This matrix is computed in \verb+computeAdjugateMatrix+ in
762\verb+SnowDriftFEnumerics.cc+.
763
764
765
766% ------------------------------------------------------------
767% ------------------------------------------------------------
768\subsection{Boundary conditions}
769% ------------------------------------------------------------
770% ------------------------------------------------------------
771As opposed to \cite{ryser_2004} boundary conditions are now
772implemented as weak Robin boundary conditions which has the advantage
773of higher flexibility.
774
775% ------------------------------------------------------------
776\subsubsection{Additional element contributions}
777% ------------------------------------------------------------
778The derivation of the weak formulation of Eq.\ (25) in
779\cite{ryser_2004} (i.e. by multiplying with a test function $v$ and
780integrating by parts) in general leads to an additional boundary term,
781viz
782\begin{eqnarray}\label{intbyparts}
783 -\int_{\Omega} \mathrm{d}x\ \bigl[ \nabla\cdot(K(x)\nabla c(x)) \bigr] v(x)
784 &=& \int_{\Omega} \mathrm{d}x\ \bigl[K(x)\nabla c(x)\bigr]\cdot\nabla v(x)\\
785 &\phantom{=}&-\int_{\partial\Omega} \mathrm{d}a\ n(x)\cdot\bigl[K(x)\nabla c(x)\bigr] v(x)\;,
786\end{eqnarray}
787where $n(x)$ is the \emph{outward} normal vector field on the boundary
788$\partial\Omega$ of the domain $\Omega$.
789
790The most general Robin boundary condition (also: generalized Neumann
791or flux boundary condition) can be defined in the form
792\begin{eqnarray}\label{rbc}
793 -n(x)\cdot\bigl[K(x)\nabla c(x)\bigr]=\gamma(x)\big(c(x)-g_{\mathrm{D}}(x)\big)+g_{\mathrm{N}}(x)
794\end{eqnarray}
795in terms of functions $\gamma,g_{\mathrm{D}},g_{\mathrm{N}}$ which are
796defined on the boundary $\partial\Omega$. The limiting case of pure
797Dirichlet conditions $c(x)=g_{\mathrm{D}}(x)$ can be formally obtained
798by $\gamma(x)\to\infty$ and pure Neumann conditions
799$-n(x)\cdot\bigl[K(x)\nabla c(x)\bigr]=g_{\mathrm{N}}$ by setting
800$\gamma(x)=0$.
801
802Consequently, the weak imposition of the Robin boundary condition
803(\ref{rbc}) in Eq. 25 of \cite{ryser_2004} can be achieved by
804inserting (\ref{rbc}) into (\ref{intbyparts}) which leads to two
805additional terms in the weak formulation (27) in \cite{ryser_2004}: An
806additional contribution of the load on the rhs of Eq.\ (27) in
807\cite{ryser_2004} is given by (counting minus-signs correctly)
808\begin{eqnarray}\label{add1}
809 b_r(v) := \int_{\partial\Omega} \mathrm{d}a\ \bigl[\gamma(x)g_{\mathrm{D}}(x)-g_{\mathrm{N}}(x)\bigr]v(x)
810\end{eqnarray}
811And the additional contribution from the Robin condition to the
812element matrix on the lhs of Eq.\ (27) in \cite{ryser_2004} reads
813\begin{eqnarray}\label{add2}
814 b_l(c,v) := \int_{\partial\Omega} \mathrm{d}a\ \gamma(x)c(x)v(x)
815\end{eqnarray}
816Note, to obtain the consistently stabilized conterpart of Eq.\ (31) in
817\cite{ryser_2004} with the Robin boundary conditions, both additional
818terms (\ref{add1},\ref{add2}) have to be evaluated with stabilized test
819functions.
820
821
822% ------------------------------------------------------------
823\subsubsection{Integrals over element surfaces}
824% ------------------------------------------------------------
825The respective element contributions of the surface integrals
826(\ref{add1},\ref{add2}) are also computed by transformation on the
827reference element $\tilde{K}$. The isoperimetric transformation $F$ in
828\ref{isoperimetric} induces six mappings $G_{\alpha,\sigma}$ from the
829surfaces
830\begin{equation}
831 \tilde{A}_{\alpha,\sigma}=
832 \{\xi \in \tilde{K} : \xi_{\alpha}=\sigma,\,
833 \sigma= \pm 1,\,
834 \alpha= 1,2,3\}
835\end{equation}
836of the reference element $\tilde{K}$ to the surfaces of the
837${A}_{\alpha,\sigma}$ of $K$ given by the restriction of $F$ on the
838respective surface $G_{\alpha,\sigma}:=F|_{\tilde{A}_{\alpha,\sigma}}$
839
840Then, the usual transformation formula for surface integrals applies
841to the surfaces of the element $K$ via
842
843\begin{equation}
844 \int_{{A}_{\alpha,\sigma}}\mathrm{d}a f(x,y,z)
845 =
846 \int_{\tilde{A}_{\alpha,\sigma}}
847 \mathrm{d}\xi_{\alpha}
848 \mathrm{d}\xi_{\gamma}\
849 \|\partial_{\xi_{\beta}} F(\xi) \times \partial_{\xi_{\gamma}}
850 F(\xi)\|\Big|_{\xi_{\alpha}=\sigma}
851 f \bigl(F(\xi)\bigr)\Big|_{\xi_{\alpha}=\sigma}
852\end{equation}
853and $\alpha,\beta,\gamma$ chosen cyclic in $1,2,3$. The surface
854jacobian can be expressed in terms of the adjugate $\tilde{J}$ of the
855jacobian via
856\begin{equation}
857 \|\partial_{\xi_{\beta}} F(\xi) \times \partial_{\xi_{\gamma}} F(\xi)\|
858 = \Big( \sum_{\beta} \left[\tilde{J}_{\beta,\alpha}\right]^2 \Big)^{1/2}
859\end{equation}
860The application of Robin boundary conditions is done in the method\\
861\verb+SnowDriftWorker::applyRobinBC+ in \verb+SnowDriftFEControl.cc+.
862The functions $\gamma,g_{\mathrm{D}},g_{\mathrm{N}}$ can be specified
863in \verb+SnowDrift::InitializeSystem+ in \\\verb+SnowDriftFEInit.cc+.
864
865% ------------------------------------------------------------
866% ------------------------------------------------------------
867\subsection{Testing the numerics}
868% ------------------------------------------------------------
869% ------------------------------------------------------------
870By choosing the diffusion $K=\mathrm{diag}(0,0,K_z)$, the flow
871$\mathbf{u}=(0,0,w)$ and the source term as $f(z)=mz$ with
872constants $m,K_z,w$ the stationary problem reduces to the one-dimensional
873advection diffusion equation
874\begin{equation}
875 \label{eq:1ddiff}
876 -K_z\ c''(z)+uc'(z)=mz
877\end{equation}
878which can be solved analytically. It can be easiliy verified that with
879Dirichlet boundary condition $c(L)=c_L$ at the top of column and a
880Neumann boundary condition $c'(0)=\zeta$ the solution can be written
881as
882\begin{equation}
883 \label{eq:1ddiff}
884 c(z)=AK/u\ \exp\{uz/K\}+B+m(Kz+uz^2/2)/u^2
885\end{equation}
886where the constants $A,B$ are determined by the boundary conditions
887via
888\begin{equation}
889 \label{eq:1ddiff}
890 A=\zeta-km/u^2,\quad B=c_L-m(KL+uL^2/2)/u^2-AK/u\exp\{uL/K\}
891\end{equation}
892The comparison between the exact solution and the numerical solution
893for values $K_z=100,u=1,m=1,c_L=5,\zeta=1,\Delta z=1,L=32$ is given in
894Fig.\ \ref{test1}.
895
896\begin{figure}[t]
897 \begin{minipage}[c]{0.48\textwidth}
898% \begin{center}
899 \includegraphics[angle=-90,width=7cm]{\figuredir/compHighK.eps}
900% \end{center}
901 \caption{Comparison between exact and numerical solution for an
902 diffusion dominated case}
903 \label{test1}
904 \end{minipage}
905 \hfill
906 \begin{minipage}[c]{0.48\textwidth}
907% \begin{center}
908 \includegraphics[angle=-90,width=7cm]{\figuredir/compLowK.eps}
909% \end{center}
910 \caption{Comparison between exact and numerical solution for an
911 advection dominated case}
912 \label{test2}
913 \end{minipage}
914\end{figure}
915The imposition of the Dirichlet condition has been done by setting
916$\gamma(0)=10^7$. Note, that different values of $K,u,etc$ might
917require to adjust $\gamma(0)$ to a higher or lower value in order get
918a suitably well conditioned linear system of equations.
919
920For an advection dominated case we choose
921$K_z=1,u=10,m=1,c_L=5,\zeta=-1,\Delta z=0.01,L=32$ and obtain Fig.\
922\ref{test2} Note, that the solution remains stable (no oscillations)
923also for quite larger grid spacings $\Delta z$ whereas the error is
924drastically increased. It is recommended to apply {\em a posteriori}
925error bounds to the advection diffusion problem for the particular
926cases of interest to better judge the quality of the solution.
927
928
929
930%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
931% ------------------------------------------------------------
932\section{Details of the implementation of parallelism}
933% ------------------------------------------------------------
934%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
935
936% ------------------------------------------------------------
937% ------------------------------------------------------------
938\subsection{Background}
939% ------------------------------------------------------------
940% ------------------------------------------------------------
941
942% ------------------------------------------------------------
943\subsubsection{General}
944% ------------------------------------------------------------
945
946As detailed elsewhere, the snowdrift model is basically an advection
947diffusion equation for a passive tracer. The old implementation of the
948model was an explicit Euler scheme for the time integration of the
949transient advection diffusion equation. Its parallelization was easily
950accomplished by domain decomposition. Due to the advection dominated
951behaviour of the equation the scheme was however unstable and an
952implicit solution was desired. Initiated by Marc Ryser (cf
953\cite{ryser_2004}) an semi-implicit scheme (SUPG with Crank Nicolson
954time stepping) was implemented which requires the solution of a large,
955unsymmetric, sparse linear system. For an efficient parallel
956implementation of the solution the whole workflow of the finite
957element algorithm should be parallelized, i.e. parallel data (vectors,
958matrices) parallel assembling (domain decomposition) and, most
959important for efficiency, a parallel solver (stabilized bi-conjugate
960gradient)
961
962
963
964% ------------------------------------------------------------
965\subsubsection{MPI/PETSc}
966% ------------------------------------------------------------
967Since the GRID middleware POPC++ (cf \cite{popc}) does not come with
968parallel libraries for standard numerical problems (linear algebra) we
969was either forced to re-invent the wheel or to combine
970well-established libraries with POPC++.
971
972An appropriate starting point is PETSc (cf. \cite{petsc}) which is an
973MPI based, parallel-numerics library. Standalone example programs are
974included in the documentation of PETSc such that a standalone parallel
975snowdrift model can be implemented in a straightforward manner. Here
976the main difficulty remains to unify the parallelization in such a
977manner that the inter-module GRID parallelism achieved by POPC++
978(ebalance, snowpack, snowdrift, runoff, etc) can be maintained while
979introducing additional intra-module parallelism within the snowdrift
980module with MPI/PETSc.
981
982
983% ------------------------------------------------------------
984\subsubsection{POPC++ version 1.1}
985% ------------------------------------------------------------
986A new (temporary) version of POPC++ was released which regards certain
987modules as MPI/PETSc processes. In this way, ``workers'' of a
988particular modules (such as in snowpack, snowdrift) can be initialized
989as MPI/PETSc processes and standard MPI/PETSc-code can be simply used
990in the worker class. Basically, this can be achieved by a new
991construction mechanism of the SnowDriftWorker class and eventually
992using the \verb+-''object=petsc''+-flag during linking of the
993\verb+snowdriftworker.module+ In this way the GRID-parallel class
994\verb+SnowDrift+ simply spawns its workers as MPI/PETSc processes.
995For details I refer to the source code (interface of
996\verb+SnowDriftParallel+ and the construction of the workers in
997\verb+SnowDriftParallel.cc+)
998
999
1000% ------------------------------------------------------------
1001% ------------------------------------------------------------
1002\subsection{Implementation details}
1003% ------------------------------------------------------------
1004% ------------------------------------------------------------
1005
1006% ------------------------------------------------------------
1007\subsubsection{Node and element enumeration}
1008% ------------------------------------------------------------
1009Ideally the complete parallelization should be done by using solely
1010PETSc which support appropriate data structures with ghost nodes and
1011mapping between different enumeration schemes.
1012
1013However, here the old structure of the domain decomposition has been
1014maintained, where the whole simulation domain, is cut (by PopC++) into
1015approximately equal chunks along planes parallel to the $yz$
1016coordinate plane. The implementation always assumes a regular grid
1017where $N_x,N_y,N_z$ are the {\em global} number of nodes in each
1018coordinate direction. Within such a domain decomposition the
1019(processor-) local numbers of nodes in each coordinate direction,
1020denoted by $n_x,n_y,n_z$ are related to the global values by $N_y=n_y$
1021and $N_z=n_z$ and only the local $n_x$ is different from $N_x$. The
1022illustration of the domain decomposition is depicted in Fig.\
1023\ref{decomp}.
1024\begin{figure}[t]
1025 \label{decomp}
1026 \begin{center}
1027 \noindent\includegraphics[angle=0,width=\textwidth]{\figuredir/domainDecomp.eps}
1028 \end{center}
1029 \caption{Domain decomposition, here for an example with three
1030 SnowDriftWorkers and $N_x=24$. Only the bottom layer of nodes of
1031 the whole simulation domain in the $xy$ plane is shown. Some
1032 examples values of the variables used in the code are shown.}
1033\end{figure}
1034The overlap (or ghost nodes) between adjacent domains is determined by
1035PopC and contains two layers of nodes. Thus, the global $N_x$ is not
1036simply the sum over the local $n_x$ due to the overlap.
1037
1038For the present decomposition it is advantageous to use a node
1039enumeration such that consecutive blocks of node-numbers are hold by
1040one processor. This leads to a global node enumeration as shown in
1041Fig.\ \ref{globalNode}: The global mapping of a node with global
1042coordinates $I_x,I_y,I_z$ onto an global node index $I$ is given by
1043$I=I_xN_zN_y+I_zN_y+I_y$, i.e. the lattice is traversed along
1044coordinate directions in the order $y,z,x$. Global node indices for an
1045element with local element number are available on each processor via
1046the \verb+globalNodeMap+ array which is defined in
1047\verb+SnowDriftWorker.cc+ It maps the node numbers $\in\{0,1..8\}$
1048within an element with local element number
1049$\in\{0,1..(n_x-1)(n_y-1)(n_z-1)-1\}$ to the global node index which
1050is e.g. $\in\{5N_yN_z,..18N_yN_z-1\}$ for Worker no 2 in the three
1051processor example given above. Note, that the elements are enumerated
1052by traversing the lattice in along coordinate directions in the order
1053$x,y,z$.
1054
1055In contrast, locally the nodes are enmuerated also in the
1056$x,y,z$-scheme which was the old enumeration scheme and which was left
1057for data structure which survived in the present parallel
1058implementation (i.e. the \verb+nodes+-array). More precisely the
1059enumeration is according to $i=i_x+i_yn_x+i_zn_yn_x$ for a node with
1060local coordinates $i_x,i_y,i_z$ and local index $i$. The
1061\verb+nodeMap+-array is defined in \verb+SnowDrift.cc+ such that node
1062numbers $\in\{0,1..8\}$ within an element with a local element number
1063$\in\{0,1..(n_x-1)(n_y-1)(n_z-1)-1\}$ is mapped onto the the local
1064node index $\in\{0,1..n_xn_yn_z-1\}$. Again, the elements
1065are enumerated in the $x,y,z$-scheme.
1066
1067\begin{figure}[t]
1068 \begin{minipage}[c]{0.48\textwidth}
1069 \includegraphics[angle=0,width=7cm]{\figuredir/globalNodeEnumeration.eps}
1070 \caption{Global node enumeration on the example of Worker 2 from
1071 the domain decomposition in Fig.\ \ref{decomp}}
1072 \label{globalNode}
1073 \end{minipage}
1074 \hfill
1075 \begin{minipage}[c]{0.48\textwidth}
1076 \includegraphics[angle=0,width=7cm]{\figuredir/localNodeEnumeration.eps}
1077 \caption{Local node enumeration on the example of Worker 2 from
1078 the domain decomposition in Fig.\ \ref{decomp}}
1079 \label{localNode}
1080 \end{minipage}
1081\end{figure}
1082
1083
1084% ------------------------------------------------------------
1085\subsubsection{Communication}
1086% ------------------------------------------------------------
1087During the assembling and the solution of the linear system the
1088communication is automatically done by PETSc. The values of the
1089solution vector on the processor domain boundary are afterwards
1090distributed to neighboring processors by POPC++ in the
1091\verb+ExchangeBoundSuspension+ method in \verb+SnowDriftWorker.cc+.
1092This is necessary to compute the deposition flux in
1093\verb+SnowDriftFEControl.cc+ consistently.
1094
1095
1096% ------------------------------------------------------------
1097% ------------------------------------------------------------
1098\subsection{Problems with the present parallel drift}
1099% ------------------------------------------------------------
1100% ------------------------------------------------------------
1101
1102% ------------------------------------------------------------
1103\subsubsection{No unifying source code for sequential and parallel mode}
1104% ------------------------------------------------------------
1105The main problem of the present parallel implementation is that it
1106cannot be run sequentially. To explain this difficulty it is necessary
1107to understand the following details of the implementation: Every
1108MPI/PETSc program requires the initialization of the parallel
1109communicator via \verb+PetscInitialize()+ which in turn calls
1110\verb+MPI_Init()+. This is implicitely done by POPC++ during
1111construction of the SnowDriftWorkers. This implies i) all PETSc
1112variables (vectors,matrices) have to be members of the
1113\verb+SnowDriftWorker+-class and ii) all workers are automatically
1114ready to execute MPI/PETSc code.
1115
1116This reveals the two problems when one wants to take over this
1117implementation to a sequential mode: Presently, in sequential mode no
1118workers are constructed. Consequently i) no PETSc variables are
1119available at all ii) and no call to \verb+PetscInitialize()+ is done.
1120
1121A workaround might be possible by (a lot of?) conditional compilation:
1122First, in sequential mode \verb+PetscInitialize()+ must be invoked
1123manually somewhere else (where? probably in AlpineMain). Second, in
1124sequential mode all PETSc variables must be members of the class
1125\verb+SnowDrift+ and {\em not} members of the class
1126\verb+SnowDriftWorker+. Vice versa in parallel mode.
1127
1128Therefore, the parallel version is included in an additional directory
1129\verb+snowdrift_par+ in the repository.
1130
1131
1132% ------------------------------------------------------------
1133\subsubsection{Dynamic library completion}
1134% ------------------------------------------------------------
1135It is necessary to additionally hack the snowdriftworker.module after
1136running \verb+make deploy_par+ in order to appropriately include the
1137dynamic libraries. This is achieved by the helper application
1138\verb+parocexe+ which is located in the \verb+tools+ directory located
1139in your \verb+trunk+ directory. The call to parocexe is included in
1140the \verb+Makefile.par+.
1141
1142
1143% ------------------------------------------------------------
1144% ------------------------------------------------------------
1145\subsection{If you want to use it}
1146% ------------------------------------------------------------
1147% ------------------------------------------------------------
1148\begin{enumerate}
1149\item adjust your environment:
1150 \verb+source ./tools/popc-petsc.env+
1151\item \verb+make all_par -f Makefile.par+
1152\item \verb+make deploy_par -f Makefile.par+
1153\item \verb+./sc_drift_par.sh+
1154\end{enumerate}
1155
1156With the present ``installation'' of popc only the MPI processes are
1157computed in parallel, the remaining modules are running on your local
1158machine.
1159
1160
1161% ------------------------------------------------------------
1162% ------------------------------------------------------------
1163\section{TODO}
1164% ------------------------------------------------------------
1165% ------------------------------------------------------------
1166\begin{enumerate}
1167\item install a proper version of popc-1.1 in /usr/local against mpi
1168 (should be the version \verb+/usr/local/mpich-1.2.5.2+ since
1169 petsc-2.3.0 is compiled against it) and \verb+petsc-2.3.0-popc+
1170\end{enumerate}
1171
1172
1173
1174%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1175%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1176% ------------------------------------------------------------
1177\chapter{Data Assimilation Module}\label{ch:assimilation}
1178% ------------------------------------------------------------
1179%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1180%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1181
1182%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1183% ------------------------------------------------------------
1184\section{Files}
1185% ------------------------------------------------------------
1186%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1187The data assimilation module in \verb+./assimilation+
1188contains the files
1189\begin{verbatim}
1190 DataAssimilation.h
1191 DataAssimilation.ph
1192 DataAssimilation.cc
1193 PackDataAssimilation.cc
1194\end{verbatim}
1195
1196
1197%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1198% ------------------------------------------------------------
1199\section{Description}
1200% ------------------------------------------------------------
1201%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1202
1203The assimilation module (often abbreviated by ``DA'' or ``da'' in the
1204code) is responsible for any data assimilation specific computations.
1205Presently, it basically hosts memory for the assimilation data which
1206could have been included also in some of the other modules since only
1207direct insertion like assimilation without any computation is
1208implemented. It has been decided to add an additional module though
1209for the sake of generality since for more sophisticated assimilation
1210schemes (like a Kalman Filter) extensive numerics might be required
1211(i.e. for the computation of the Kalman gain matrix).
1212
1213The module can be enabled by the switch \verb+-enable-da+ in the
1214command line options of Alpine3d limilar to the other modules.
1215Additionally, a path to the assimilation data has to be provided
1216by adding the switch
1217
1218\verb+-dadir=<path to da-data, without terminating slash>+
1219
1220Having enabled the data assimilation the basic working steps of
1221Alpine3d which involve the assimilation are given below:
1222
1223\begin{enumerate}
1224\item At each time step the \verb+Compute()+ method of the
1225 DataAssimlation class is invoked in \verb+AlpineControl::Run+ in
1226 \verb+AlpineControl.cc+ which forces its input member to try to read
1227 assimilation data from the specified directory \verb+dadir+. In that
1228 directory, the program expects data in files with names
1229 \verb+<YYYYMMDDHH>.sca+ which must contain integer data and should
1230 be formatted according to the usual Ascii ArcInfo grid format.
1231\item If data is available in that directory with the correct
1232 timestamp it is read into the memory of the assimilation class and
1233 send to snowpack, which itself distributes the data to its workers
1234 (if run in parallel). If no DA-data is available at that timestep an
1235 exception is thrown and everything continues regularly.
1236\item Within snowpack, \verb+Snowpack::Assimlate()+ is invoked and the
1237 DA data is used to control actions in Snowpack. For the present
1238 example in the case of SCA data cf \verb+Snowpack::Assimilate()+ in
1239 \verb+SnowInterface.cc+ for the details of the actions.
1240\end{enumerate}
1241
1242
1243
1244%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1245% ------------------------------------------------------------
1246\section{Modifications in other modules}
1247% ------------------------------------------------------------
1248%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1249The modification in other files which were necessary to build in a new
1250class (parclass) are listed below:
1251\begin{itemize}
1252\item add ./assimlation directory
1253\item modify all Makefiles (for obvious reasons)
1254\item in ./main (creating the DA module in AlpineControl)
1255\item in ./common (add new exceptions, etc)
1256\item in ./snowpack: add pointer to assimilation module, add
1257 assimilation routines, add distribution of DA data to workers
1258\item in ./inout: add a method for reading DA-data and a method for
1259 getting the grid dimensions (GetGridSize(...))
1260\end{itemize}
1261
1262%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1263%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1264\begin{thebibliography}{99}
1265%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1266%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1267\bibitem{popc}
1268 See documentation on \\
1269 (\verb+http://gridgroup.tic.hefr.ch/popc/index.php/Main_Page+)
1270
1271\bibitem{lehning_2006}
1272 Lehning, M., V\"olksch, I., Gustafsson, D., Nguyen, T.A., St\"ahli, M.,
1273 Zappa, M., (2006)
1274 ALPINE3D: A detailed model of mountain surface processes and its
1275 application to snow hydrology,
1276 \textit{Hydrol. Processes, 20}, 2111-2128.
1277
1278\bibitem{petsc}
1279 See documentation on \\
1280 (\verb+http://www-unix.mcs.anl.gov/petsc/petsc-as/+)
1281
1282\bibitem{ryser_2004}
1283 Ryser, M.
1284 Numerical Simulation of Snow Transport Above Alpine Terrain
1285 Internship Report
1286 see \verb+MarcRyser_Drift_Final.pdf+
1287
1288\bibitem{gt}
1289 See documentation on \\
1290 (\verb+http://www.globus.org/toolkit/+)
1291
1292\bibitem{svn}
1293 Online book on\\
1294 (\verb+http://svnbook.red-bean.com+)
1295
1296\bibitem{jon07}
1297 Jonas, T. (SLF)
1298 ``Alpine3d Crash Course for the Snow Hydrology Research Group''
1299
1300\end{thebibliography}
1301
1302
1303
1304
1305
1306\end{document}

Archive Download this file

Revision: HEAD