Changeset 7086
- Timestamp:
- May 26, 2009, 11:33:53 AM (16 years ago)
- Location:
- anuga_core/documentation/user_manual
- Files:
-
- 4 edited
-
anuga_user_manual.tex (modified) (30 diffs)
-
demos/cairns/runcairns.py (modified) (1 diff)
-
demos/channel1.py (modified) (2 diffs)
-
demos/runup.py (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/documentation/user_manual/anuga_user_manual.tex
r7071 r7086 29 29 \usepackage[english]{babel} 30 30 \usepackage{datetime} 31 \usepackage[hang,small,bf]{caption} 31 32 32 33 \input{definitions} … … 158 159 159 160 The latest installation instructions may be found at: 160 \url{http://datamining.anu.edu.au/ ~ole/anuga/user_manual/anuga_installation_guide.pdf}.161 \url{http://datamining.anu.edu.au/\~{}ole/anuga/user_manual/anuga_installation_guide.pdf}. 161 162 162 163 \section{Audience} … … 320 321 \subsection{The Code} 321 322 322 %FIXME: we are using the \code function here.323 %This should be used wherever possible324 323 For reference we include below the complete code listing for 325 324 \file{runup.py}. Subsequent paragraphs provide a … … 327 326 significance. 328 327 328 \label{ref:runup_py_code} 329 329 \verbatiminput{demos/runup.py} 330 330 … … 384 384 \end{verbatim} 385 385 386 In addition, the following statement now specifies that the386 In addition, the following statement could be used to state that 387 387 quantities \code{stage}, \code{xmomentum} and \code{ymomentum} are 388 388 to be stored: … … 391 391 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 392 392 \end{verbatim} 393 394 However, this is not necessary, as the above quantities are always stored by default. 393 395 394 396 \subsection{Initial Conditions} … … 416 418 417 419 \begin{verbatim} 418 def f(x,y):420 def topography(x, y): 419 421 return -x/2 420 422 \end{verbatim} … … 425 427 the \code{y} direction. 426 428 427 Once the function \function{ f} is specified, the quantity429 Once the function \function{topography} is specified, the quantity 428 430 \code{elevation} is assigned through the simple statement: 429 431 430 432 \begin{verbatim} 431 domain.set_quantity('elevation', f)433 domain.set_quantity('elevation', topography) 432 434 \end{verbatim} 433 435 434 436 NOTE: If using function to set \code{elevation} it must be vector 435 compatible. For example square root will not work.437 compatible. For example, using square root will not work. 436 438 437 439 \subsubsection{Friction} … … 465 467 466 468 which specifies that the surface level is set to a height of $-0.4$, 467 i.e. 0.4 units (m ) below the zero level.469 i.e. 0.4 units (metres) below the zero level. 468 470 469 471 Although it is not necessary for this example, it may be useful to … … 496 498 \begin{verbatim} 497 499 Br = Reflective_boundary(domain) 498 499 500 Bt = Transmissive_boundary(domain) 500 501 501 Bd = Dirichlet_boundary([0.2, 0.0, 0.0]) 502 503 Bw = Time_boundary(domain=domain, f=lambda t: [(0.1*sin(t*2*pi)-0.3), 0.0, 0.0])502 Bw = Time_boundary(domain=domain, 503 f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0]) 504 504 \end{verbatim} 505 505 … … 511 511 follows: 512 512 \begin{itemize} 513 \item \textbf{Reflective boundary}\label{def:reflective boundary} Returns same \code{stage} as514 as present in its neighbour volume but momentum vector515 reversed 180 degrees (reflected).513 \item \textbf{Reflective boundary}\label{def:reflective boundary} 514 Returns same \code{stage} as in its neighbour volume but momentum 515 vector reversed 180 degrees (reflected). 516 516 Specific to the shallow water equation as it works with the 517 517 momentum quantities assumed to be the second and third conserved … … 615 615 616 616 \begin{verbatim} 617 for t in domain.evolve(yieldstep=0.1, duration= 4.0):617 for t in domain.evolve(yieldstep=0.1, duration=10.0): 618 618 print domain.timestepping_statistics() 619 619 \end{verbatim} … … 629 629 The output is a NetCDF file with the extension \code{.sww}. It 630 630 contains stage and momentum information and can be used with the 631 ANUGA viewer \code{animate} visualisation packageto generate a visual631 ANUGA viewer \code{animate} to generate a visual 632 632 display (see Section \ref{sec:animate}). See Section \ref{sec:file formats} 633 633 (page \pageref{sec:file formats}) for more on NetCDF and other file … … 661 661 on the previously dry bed. 662 662 663 All figures are screenshots from an interactive animation tool called \code{animate} 664 which is part of \anuga and distributed as in the package anuga\_viewer. 665 \code{animate} is described in more detail is Section \ref{sec:animate}. 663 \code{animate} is described in more detail in Section \ref{sec:animate}. 666 664 667 665 \begin{figure}[htp] … … 713 711 714 712 Here is the code for the first version of the channel flow \file{channel1.py}: 713 715 714 \verbatiminput{demos/channel1.py} 715 716 716 In discussing the details of this example, we follow the outline 717 717 given above, discussing each major step of the code in turn. … … 747 747 as we will later in this example, one merely decreases the values of \code{dx} and \code{dy}. 748 748 749 The rest of this script is as in the previous example.749 The rest of this script is similar to the previous example on page \pageref{ref:runup_py_code}. 750 750 % except for an application of the 'expression' form of \code{set\_quantity} where we use 751 751 % the value of \code{elevation} to define the (dry) initial condition for \code{stage}: … … 784 784 785 785 \begin{verbatim} 786 for t in domain.evolve(yieldstep = 0.2, finaltime =40.0):786 for t in domain.evolve(yieldstep=0.2, finaltime=40.0): 787 787 domain.write_time() 788 788 … … 793 793 794 794 \label{sec:change boundary code} 795 The \code{if} statement in the timestepping loop ( evolve) gets the quantity795 The \code{if} statement in the timestepping loop (\code{evolve}) gets the quantity 796 796 \code{stage} and obtains the interpolated value at the point (10m, 797 797 2.5m) which is on the right boundary. If the stage exceeds 0m a … … 923 923 using a mesh-generator. 924 924 925 In its simplest form, the mesh-generator creates the mesh within a single925 The mesh-generator creates the mesh within a single 926 926 polygon whose vertices are at geographical locations specified by 927 927 the user. The user specifies the \emph{resolution} -- that is, the 928 maximal area of a triangle used for triangulation R-- and a triangular928 maximal area of a triangle used for triangulation -- and a triangular 929 929 mesh is created inside the polygon using a mesh generation engine. 930 930 On any given platform, the same mesh will be returned. 931 %Figure932 %\ref{fig:pentagon} shows a simple example of this, in which the933 %triangulation is carried out within a pentagon.934 935 %\begin{figure}[htp]936 % \caption{Mesh points are created inside the polygon}937 %\label{fig:pentagon}938 %\end{figure}939 931 940 932 Boundary tags are not restricted to \code{'left'}, \code{'bottom'}, … … 951 943 polygons in which no triangulation is required. 952 944 953 %\begin{figure}[htp]954 % \caption{Interior meshes with individual resolution}955 % \label{fig:interior meshes}956 %\end{figure}957 958 945 In its general form, the mesh-generator takes for its input a bounding 959 946 polygon and (optionally) a list of interior polygons. The user … … 963 950 964 951 The function used to implement this process is 965 \function{create\_mesh\_from\_regions}. Its arguments include the 952 \function{create\_domain\_from\_regions} which creates a Domain object as 953 well as a mesh file. Its arguments include the 966 954 bounding polygon and its resolution, a list of boundary tags, and a 967 list of pairs \code{[polygon, resolution]} ,specifying the interior955 list of pairs \code{[polygon, resolution]} specifying the interior 968 956 polygons and their resolutions. 969 957 … … 1033 1021 \begin{verbatim} 1034 1022 remainder_res = 10000000 1035 create_mesh_from_regions(project.bounding_polygon,1036 boundary_tags={'top': [0],1037 'ocean_east': [1],1038 'bottom': [2],1039 'onshore': [3]},1040 maximum_triangle_area=remainder_res,1041 filename=meshname,1042 interior_regions=interior_regions,1043 use_cache=True,1044 verbose=True)1023 domain = create_domain_from_regions(project.bounding_polygon, 1024 boundary_tags={'top': [0], 1025 'ocean_east': [1], 1026 'bottom': [2], 1027 'onshore': [3]}, 1028 maximum_triangle_area=project.default_res, 1029 mesh_filename=project.meshname, 1030 interior_regions=project.interior_regions, 1031 use_cache=True, 1032 verbose=True) 1045 1033 \end{verbatim} 1046 1034 … … 1056 1044 (Here, the values associated with each boundary tag are one-element lists, but they can have as many indices as there are edges) 1057 1045 If polygons intersect, or edges coincide (or are even very close) the resolution may be undefined in some regions. 1058 Use the underlying mesh interface for such cases .1059 See Section \ref{sec:mesh interface}.1046 Use the underlying mesh interface for such cases 1047 (see Chapter \ref{sec:mesh interface}). 1060 1048 If a segment is omitted in the tags definition an Exception is raised. 1061 1049 … … 1067 1055 \subsection{Initialising the Domain} 1068 1056 1069 As with \file{runup.py}, once we have created the mesh, the next 1070 step is to create the data structure \code{domain}. We did this for 1071 \file{runup.py} by inputting lists of points and triangles and 1072 specifying the boundary tags directly. However, in the present case, 1073 we use a method that works directly with the mesh file 1074 \code{meshname}, as follows: 1057 Since we used \code{create_domain_from_regions} to create the mesh file, we do not need to 1058 create the domain explicitly, as the above function also does that. 1059 1060 The following statements specify a basename and data directory, and 1061 sets a minimum storable height, which helps with visualisation. 1075 1062 1076 1063 \begin{verbatim} 1077 domain = Domain(meshname, use_cache=True, verbose=True) 1078 \end{verbatim} 1079 1080 Providing a filename instead of the lists used in \file{runup.py} 1081 above causes \code{Domain} to convert a mesh file \code{meshname} 1082 into an instance of \code{Domain}, allowing us to use methods like 1083 \method{set\_quantity} to set quantities and to apply other 1084 operations. 1085 1086 %(In principle, the 1087 %second argument of \function{pmesh\_to\_domain\_instance} can be any 1088 %subclass of \class{Domain}, but for applications involving the 1089 %shallow-water wave equation, the second argument of 1090 %\function{pmesh\_to\_domain\_instance} can always be set simply to 1091 %\class{Domain}.) 1092 1093 The following statements specify a basename and data directory, and 1094 identify quantities to be stored. For the first two, values are 1095 taken from \file{project.py}. 1096 1097 \begin{verbatim} 1098 domain.set_name(project.basename) 1099 domain.set_datadir(project.outputdir) 1100 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 1064 domain.set_name('cairns_' + project.scenario) # Name of sww file 1065 domain.set_datadir('.') # Store sww output here 1066 domain.set_minimum_storable_height(0.01) # Store only depth > 1cm 1101 1067 \end{verbatim} 1102 1068 … … 1111 1077 \subsubsection{Stage} 1112 1078 1113 For the scenario we are modelling in this case, we use a callable 1114 object \code{tsunami\_source}, assigned by means of a function 1115 \function{slide\_tsunami}. This is similar to how we set elevation in 1116 \file{runup.py} using a function -- however, in this case the 1117 function is both more complex and more interesting. 1118 1119 The function returns the water displacement for all \code{x} and 1120 \code{y} in the domain. The water displacement is a double Gaussian 1121 function that depends on the characteristics of the slide (length, 1122 width, thickness, slope, etc), its location (origin) and the depth at that 1123 location. For this example, we choose to apply the slide function 1124 at a specified time into the simulation. {\bf Note, the parameters used 1125 in this example have been deliberately chosen to generate a suitably 1126 large amplitude tsunami which would inundate the Cairns region.} 1079 The stage is initially set to 0.0 by the following statements: 1080 1081 \begin{verbatim} 1082 tide = 0.0 1083 domain.set_quantity('stage', tide) 1084 \end{verbatim} 1085 1086 %For the scenario we are modelling in this case, we use a callable 1087 %object \code{tsunami_source}, assigned by means of a function 1088 %\function{slide\_tsunami}. This is similar to how we set elevation in 1089 %\file{runup.py} using a function -- however, in this case the 1090 %function is both more complex and more interesting. 1091 1092 %The function returns the water displacement for all \code{x} and 1093 %\code{y} in the domain. The water displacement is a double Gaussian 1094 %function that depends on the characteristics of the slide (length, 1095 %width, thickness, slope, etc), its location (origin) and the depth at that 1096 %location. For this example, we choose to apply the slide function 1097 %at a specified time into the simulation. {\bf Note, the parameters used 1098 %in this example have been deliberately chosen to generate a suitably 1099 %large amplitude tsunami which would inundate the Cairns region.} 1127 1100 1128 1101 \subsubsection{Friction} … … 1140 1113 \begin{verbatim} 1141 1114 domain.set_quantity('elevation', 1142 filename=project.dem _name+'.pts',1115 filename=project.demname + '.pts', 1143 1116 use_cache=True, 1144 verbose=True) 1117 verbose=True, 1118 alpha=0.1) 1145 1119 \end{verbatim} 1146 1147 %However, before this step can be executed, some preliminary steps1148 %are needed to prepare the file from which the data is taken. Two1149 %source files are used for this data -- their names are specified in1150 %the file \file{project.py}, in the variables \code{coarsedemname}1151 %and \code{finedemname}. They contain 'coarse' and 'fine' data,1152 %respectively -- that is, data sampled at widely spaced points over a1153 %large region and data sampled at closely spaced points over a1154 %smaller subregion. The data in these files is combined through the1155 %statement1156 1157 %\begin{verbatim}1158 %combine_rectangular_points_files(project.finedemname + '.pts',1159 % project.coarsedemname + '.pts',1160 % project.combineddemname + '.pts')1161 %\end{verbatim}1162 %1163 %The effect of this is simply to combine the datasets by eliminating1164 %any coarse data associated with points inside the smaller region1165 %common to both datasets. The name to be assigned to the resulting1166 %dataset is also derived from the name stored in the variable1167 %\code{combinedname} in the file \file{project.py}.1168 1120 1169 1121 \subsection{Boundary Conditions}\index{boundary conditions} … … 1174 1126 boundary tag names introduced when we established the mesh. In place of the four 1175 1127 boundary types introduced for \file{runup.py}, we use the reflective 1176 boundary for each of the eight tagged segments defined by \code{create_mesh_from_regions}:1128 boundary for each of the tagged segments defined by \code{create_domain_from_regions}: 1177 1129 1178 1130 \begin{verbatim} 1179 Bd = Dirichlet_boundary([0.0,0.0,0.0]) 1180 domain.set_boundary({'ocean_east': Bd, 'bottom': Bd, 'onshore': Bd, 'top': Bd}) 1131 Bd = Dirichlet_boundary([tide,0,0]) # Mean water level 1132 Bs = Transmissive_stage_zero_momentum_boundary(domain) # Neutral boundary 1133 1134 if project.scenario == 'fixed_wave': 1135 # Huge 50m wave starting after 60 seconds and lasting 1 hour. 1136 Bw = Time_boundary(domain=domain, 1137 function=lambda t: [(60<t<3660)*50, 0, 0]) 1138 domain.set_boundary({'ocean_east': Bw, 1139 'bottom': Bs, 1140 'onshore': Bd, 1141 'top': Bs}) 1142 1143 if project.scenario == 'slide': 1144 # Boundary conditions for slide scenario 1145 domain.set_boundary({'ocean_east': Bd, 1146 'bottom': Bd, 1147 'onshore': Bd, 1148 'top': Bd}) 1181 1149 \end{verbatim} 1182 1150 1151 Note that we use different boundary conditions depending on the \code{scenario} 1152 defined in \file{project.py}. 1153 1183 1154 \subsection{Evolution} 1184 1155 1185 1156 With the basics established, the running of the 'evolve' step is 1186 very similar to the corresponding step in \file{runup.py}. For the slide 1187 scenario, the simulation is run for 5000 seconds with the output stored every ten seconds. 1188 For this example, we choose to apply the slide at 60 seconds into the simulation: 1157 very similar to the corresponding step in \file{runup.py}, except we have different \code{evolve} 1158 loops for the two scenarios. 1159 1160 For the slide scenario, the simulation is run for an intial 60 seconds, at which time 1161 the slide occurs. We use the function \function{tsunami_source} to adjust \code{stage} 1162 values. We then run the simulation until 5000 seconds with the output stored 1163 every ten seconds. 1189 1164 1190 1165 \begin{verbatim} 1191 import time 1192 1193 t0 = time.time() 1194 1195 for t in domain.evolve(yieldstep=10, finaltime=60): 1196 domain.write_time() 1197 domain.write_boundary_statistics(tags = 'ocean_east') 1198 1199 # add slide 1166 if project.scenario == 'slide': 1167 for t in domain.evolve(yieldstep=10, finaltime=60): 1168 print domain.timestepping_statistics() 1169 print domain.boundary_statistics(tags='ocean_east') 1170 1171 # Add slide 1200 1172 thisstagestep = domain.get_quantity('stage') 1201 1173 if allclose(t, 60): … … 1205 1177 1206 1178 for t in domain.evolve(yieldstep=10, finaltime=5000, 1207 skip_initial_ste=True): 1208 domain.write_time() 1209 1210 domain.write_boundary_statistics(tags = 'ocean_east') 1179 skip_initial_step=True): 1180 print domain.timestepping_statistics() 1181 print domain.boundary_statistics(tags='ocean_east') 1182 1183 if project.scenario == 'fixed_wave': 1184 # Save every two mins leading up to wave approaching land 1185 for t in domain.evolve(yieldstep=120, finaltime=5000): 1186 print domain.timestepping_statistics() 1187 print domain.boundary_statistics(tags='ocean_east') 1188 1189 # Save every 30 secs as wave starts inundating ashore 1190 for t in domain.evolve(yieldstep=10, finaltime=10000, 1191 skip_initial_step=True): 1192 print domain.timestepping_statistics() 1193 print domain.boundary_statistics(tags='ocean_east') 1211 1194 \end{verbatim} 1212 1195 … … 1216 1199 This functionality is especially convenient as it allows the detailed 1217 1200 parts of the simulation to be viewed at higher time resolution. 1218 1219 \begin{verbatim}1220 # save every two mins leading up to wave approaching land1221 for t in domain.evolve(yieldstep=120, finaltime=5000):1222 domain.write_time()1223 domain.write_boundary_statistics(tags='ocean_east')1224 1225 # save every 30 secs as wave starts inundating ashore1226 for t in domain.evolve(yieldstep=10, finaltime=10000, skip_initial_step=True):1227 domain.write_time()1228 domain.write_boundary_statistics(tags='ocean_east')1229 \end{verbatim}1230 1231 1201 1232 1202 \section{Exploring the Model Output} … … 3730 3700 3731 3701 3732 /pagebreak3702 \pagebreak 3733 3703 \section{utilities/polygons} 3734 3704 -
anuga_core/documentation/user_manual/demos/cairns/runcairns.py
r7078 r7086 146 146 # Continue propagating wave 147 147 for t in domain.evolve(yieldstep=10, finaltime=5000, 148 skip_initial_step =True):148 skip_initial_step=True): 149 149 print domain.timestepping_statistics() 150 150 print domain.boundary_statistics(tags='ocean_east') -
anuga_core/documentation/user_manual/demos/channel1.py
r7064 r7086 24 24 # Setup initial conditions 25 25 #------------------------------------------------------------------------------ 26 def topography(x, y):26 def topography(x, y): 27 27 return -x/10 # linear bed slope 28 28 … … 30 30 domain.set_quantity('friction', 0.01) # Constant friction 31 31 domain.set_quantity('stage', # Dry bed 32 expression='elevation + 0.0')32 expression='elevation') 33 33 34 34 #------------------------------------------------------------------------------ -
anuga_core/documentation/user_manual/demos/runup.py
r7064 r7086 29 29 # Setup initial conditions 30 30 #------------------------------------------------------------------------------ 31 def topography(x, y):31 def topography(x, y): 32 32 return -x/2 # linear bed slope 33 33 #return x*(-(2.0-x)*.5) # curved bed slope … … 35 35 domain.set_quantity('elevation', topography) # Use function for elevation 36 36 domain.set_quantity('friction', 0.1) # Constant friction 37 domain.set_quantity('stage', - .4)# Constant negative initial stage37 domain.set_quantity('stage', -0.4) # Constant negative initial stage 38 38 39 39 #------------------------------------------------------------------------------ … … 44 44 Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values 45 45 Bw = Time_boundary(domain=domain, # Time dependent boundary 46 f=lambda t: [( .1*sin(t*2*pi)-0.3) *exp(-t), 0.0, 0.0])46 f=lambda t: [(0.1*sin(t*2*pi)-0.3)*exp(-t), 0.0, 0.0]) 47 47 48 48 # Associate boundary tags with boundary objects … … 52 52 # Evolve system through time 53 53 #------------------------------------------------------------------------------ 54 for t in domain.evolve(yieldstep = 0.1, finaltime =10.0):54 for t in domain.evolve(yieldstep=0.1, finaltime=10.0): 55 55 print domain.timestepping_statistics()
Note: See TracChangeset
for help on using the changeset viewer.
