Changeset 8129


Ignore:
Timestamp:
Mar 8, 2011, 9:47:19 PM (14 years ago)
Author:
steve
Message:

Getting the unit tests for structures to pass when called from test_all.py in directory anuga

Location:
trunk/anuga_core
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/compile_all.py

    r7934 r8129  
    2121os.chdir('..')
    2222os.chdir('advection')
     23execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     24
     25os.chdir('..')
     26os.chdir('operators')
     27execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     28
     29os.chdir('..')
     30os.chdir('structures')
    2331execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
    2432

  • trunk/anuga_core/source/anuga/__init__.py

    r8124 r8129  
    320320        #psyco.background() # attempt to profile code - only compile most used
    321321       
    322        
     322
     323#---------------------------
     324# Operators
     325#---------------------------
     326from anuga.operators.operator import Operator
     327
    323328#---------------------------
    324329# Structures

  • trunk/anuga_core/source/anuga/compile_all.py

    r7865 r8129  
    1212os.chdir('..')
    1313os.chdir('advection')
     14execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     15
     16os.chdir('..')
     17os.chdir('operators')
     18execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     19
     20os.chdir('..')
     21os.chdir('structures')
    1422execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
    1523

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8047 r8129  
    7171    if (dq1>=dq2){
    7272      if (dq1>=0.0)
    73         *qmax=dq0+dq1;
     73    *qmax=dq0+dq1;
    7474      else
    75         *qmax=dq0;
     75    *qmax=dq0;
    7676     
    7777      *qmin=dq0+dq2;
    … …  
    8080    else{// dq1<dq2
    8181      if (dq2>0)
    82         *qmax=dq0+dq2;
     82    *qmax=dq0+dq2;
    8383      else
    84         *qmax=dq0;
     84    *qmax=dq0;
    8585     
    8686      *qmin=dq0+dq1;   
    … …  
    9191    if (dq1<=dq2){
    9292      if (dq1<0.0)
    93         *qmin=dq0+dq1;
     93    *qmin=dq0+dq1;
    9494      else
    95         *qmin=dq0;
     95    *qmin=dq0;
    9696     
    9797      *qmax=dq0+dq2;   
    … …  
    100100    else{// dq1>dq2
    101101      if (dq2<0.0)
    102         *qmin=dq0+dq2;
     102    *qmin=dq0+dq2;
    103103      else
    104         *qmin=dq0;
     104    *qmin=dq0;
    105105     
    106106      *qmax=dq0+dq1;
    … …  
    175175//static inline double _compute_speed(double *uh,
    176176double _compute_speed(double *uh,
    177                       double *h,
    178                       double epsilon,
    179                       double h0,
    180                       double limiting_threshold) {
     177              double *h,
     178              double epsilon,
     179              double h0,
     180              double limiting_threshold) {
    181181 
    182182  double u;
    … …  
    263263                           double n1, double n2,
    264264                           double epsilon,
    265                            double h0,
    266                            double limiting_threshold,
    267                            double g,
     265               double h0,
     266               double limiting_threshold,
     267               double g,
    268268                           double *edgeflux, double *max_speed)
    269269{
    … …  
    304304  z = 0.5*(z_left + z_right); // Average elevation values.
    305305                              // Even though this will nominally allow
    306                               // for discontinuities in the elevation data,
    307                               // there is currently no numerical support for
    308                               // this so results may be strange near
    309                               // jumps in the bed.
     306                  // for discontinuities in the elevation data,
     307                  // there is currently no numerical support for
     308                  // this so results may be strange near
     309                  // jumps in the bed.
    310310
    311311  // Compute speeds in x-direction
    … …  
    314314  uh_left = q_left_rotated[1];
    315315  u_left = _compute_speed(&uh_left, &h_left,
    316                           epsilon, h0, limiting_threshold);
     316              epsilon, h0, limiting_threshold);
    317317
    318318  w_right = q_right_rotated[0];
    … …  
    320320  uh_right = q_right_rotated[1];
    321321  u_right = _compute_speed(&uh_right, &h_right,
    322                            epsilon, h0, limiting_threshold);
     322               epsilon, h0, limiting_threshold);
    323323
    324324  // Momentum in y-direction
    … …  
    330330  // All validation tests pass, so do we really need it anymore?
    331331  _compute_speed(&vh_left, &h_left,
    332                 epsilon, h0, limiting_threshold);
     332        epsilon, h0, limiting_threshold);
    333333  _compute_speed(&vh_right, &h_right,
    334                 epsilon, h0, limiting_threshold);
     334        epsilon, h0, limiting_threshold);
    335335
    336336  // Maximal and minimal wave speeds
    … …  
    459459  uh_left = q_left_rotated[1];
    460460  u_left =_compute_speed(&uh_left, &h_left,
    461                         epsilon, h0, limiting_threshold);
     461            epsilon, h0, limiting_threshold);
    462462
    463463  w_right = q_right_rotated[0];
    … …  
    465465  uh_right = q_right_rotated[1];
    466466  u_right =_compute_speed(&uh_right, &h_right,
    467                           epsilon, h0, limiting_threshold);
     467              epsilon, h0, limiting_threshold);
    468468
    469469
    … …  
    553553
    554554void _manning_friction_sloped(double g, double eps, int N,
    555                            double* x, double* w, double* zv,
    556                            double* uh, double* vh,
    557                            double* eta, double* xmom_update, double* ymom_update) {
     555               double* x, double* w, double* zv,
     556               double* uh, double* vh,
     557               double* eta, double* xmom_update, double* ymom_update) {
    558558
    559559  int k, k3, k6;
    … …  
    602602
    603603void _chezy_friction(double g, double eps, int N,
    604                            double* x, double* w, double* zv,
    605                            double* uh, double* vh,
    606                            double* chezy, double* xmom_update, double* ymom_update) {
     604               double* x, double* w, double* zv,
     605               double* uh, double* vh,
     606               double* chezy, double* xmom_update, double* ymom_update) {
    607607
    608608  int k, k3, k6;
    … …  
    785785      // FIXME: Try with this one precomputed
    786786      for (i=0; i<3; i++) {
    787         dz = max(dz, fabs(zv[k3+i]-zc[k]));
     787    dz = max(dz, fabs(zv[k3+i]-zc[k]));
    788788      }
    789789    }
    … …  
    813813     
    814814      if (dz > 0.0) {
    815         alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     815    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );
    816816      } else {
    817         alpha = 1.0;  // Flat bed
     817    alpha = 1.0;  // Flat bed
    818818      }
    819819      //printf("Using old style limiter\n");
    … …  
    828828   
    829829      if (hmin < H0) {
    830         alpha = 1.0;
    831         for (i=0; i<3; i++) {
    832          
    833           h_diff = hc_k - hv[i];     
    834           if (h_diff <= 0) {
    835             // Deep water triangle is further away from bed than
    836             // shallow water (hbar < h). Any alpha will do
    837      
    838           } else { 
    839             // Denominator is positive which means that we need some of the
    840             // h-limited stage.
    841            
    842             alpha = min(alpha, (hc_k - H0)/h_diff);     
    843           }
    844         }
    845 
    846         // Ensure alpha in [0,1]
    847         if (alpha>1.0) alpha=1.0;
    848         if (alpha<0.0) alpha=0.0;
    849        
     830    alpha = 1.0;
     831    for (i=0; i<3; i++) {
     832
     833      h_diff = hc_k - hv[i];
     834      if (h_diff <= 0) {
     835        // Deep water triangle is further away from bed than
     836        // shallow water (hbar < h). Any alpha will do
     837     
    850838      } else {
    851         // Use w-limited stage exclusively in deeper water.
    852         alpha = 1.0;       
     839        // Denominator is positive which means that we need some of the
     840        // h-limited stage.
     841
     842        alpha = min(alpha, (hc_k - H0)/h_diff);
     843      }
     844    }
     845
     846    // Ensure alpha in [0,1]
     847    if (alpha>1.0) alpha=1.0;
     848    if (alpha<0.0) alpha=0.0;
     849
     850      } else {
     851    // Use w-limited stage exclusively in deeper water.
     852    alpha = 1.0;
    853853      }
    854854    }
    … …  
    875875    if (alpha < 1) {     
    876876      for (i=0; i<3; i++) {
    877        
    878         wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
    879 
    880         // Update momentum at vertices
    881         if (use_centroid_velocities == 1) {
    882           // This is a simple, efficient and robust option
    883           // It uses first order approximation of velocities, but retains
    884           // the order used by stage.
     877
     878    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     879
     880    // Update momentum at vertices
     881    if (use_centroid_velocities == 1) {
     882      // This is a simple, efficient and robust option
     883      // It uses first order approximation of velocities, but retains
     884      // the order used by stage.
    885885   
    886           // Speeds at centroids
    887           if (hc_k > epsilon) {
    888             uc = xmomc[k]/hc_k;
    889             vc = ymomc[k]/hc_k;
    890           } else {
    891             uc = 0.0;
    892             vc = 0.0;
    893           }
    894          
    895           // Vertex momenta guaranteed to be consistent with depth guaranteeing
    896           // controlled speed
    897           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    898           xmomv[k3+i] = uc*hv[i];
    899           ymomv[k3+i] = vc*hv[i];
    900      
    901         } else {
    902           // Update momentum as a linear combination of
    903           // xmomc and ymomc (shallow) and momentum
    904           // from extrapolator xmomv and ymomv (deep).
    905           // This assumes that values from xmomv and ymomv have
    906           // been established e.g. by the gradient limiter.
    907          
    908           // FIXME (Ole): I think this should be used with vertex momenta
    909           // computed above using centroid_velocities instead of xmomc
    910           // and ymomc as they'll be more representative first order
    911           // values.
    912          
    913           xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    914           ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    915          
    916         }
     886      // Speeds at centroids
     887      if (hc_k > epsilon) {
     888        uc = xmomc[k]/hc_k;
     889        vc = ymomc[k]/hc_k;
     890      } else {
     891        uc = 0.0;
     892        vc = 0.0;
     893      }
     894
     895      // Vertex momenta guaranteed to be consistent with depth guaranteeing
     896      // controlled speed
     897      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     898      xmomv[k3+i] = uc*hv[i];
     899      ymomv[k3+i] = vc*hv[i];
     900     
     901    } else {
     902      // Update momentum as a linear combination of
     903      // xmomc and ymomc (shallow) and momentum
     904      // from extrapolator xmomv and ymomv (deep).
     905      // This assumes that values from xmomv and ymomv have
     906      // been established e.g. by the gradient limiter.
     907
     908      // FIXME (Ole): I think this should be used with vertex momenta
     909      // computed above using centroid_velocities instead of xmomc
     910      // and ymomc as they'll be more representative first order
     911      // values.
     912
     913      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     914      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     915
     916    }
    917917      }
    918918    }
    … …  
    945945      if (hc < minimum_allowed_height) {
    946946       
    947         // Set momentum to zero and ensure h is non negative
    948         xmomc[k] = 0.0;
    949         ymomc[k] = 0.0;
    950         if (hc <= 0.0) wc[k] = zc[k];
     947    // Set momentum to zero and ensure h is non negative
     948    xmomc[k] = 0.0;
     949    ymomc[k] = 0.0;
     950    if (hc <= 0.0) wc[k] = zc[k];
    951951      }
    952952    }
    … …  
    968968        } else {
    969969          //Reduce excessive speeds derived from division by small hc
    970           //FIXME (Ole): This may be unnecessary with new slope limiters
    971           //in effect.
     970      //FIXME (Ole): This may be unnecessary with new slope limiters
     971      //in effect.
    972972         
    973973          u = xmomc[k]/hc;
    974           if (fabs(u) > maximum_allowed_speed) {
    975             reduced_speed = maximum_allowed_speed * u/fabs(u);
    976             //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    977             //   u, reduced_speed);
    978             xmomc[k] = reduced_speed * hc;
    979           }
    980          
     974      if (fabs(u) > maximum_allowed_speed) {
     975        reduced_speed = maximum_allowed_speed * u/fabs(u);
     976        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     977        //   u, reduced_speed);
     978        xmomc[k] = reduced_speed * hc;
     979      }
     980
    981981          v = ymomc[k]/hc;
    982           if (fabs(v) > maximum_allowed_speed) {
    983             reduced_speed = maximum_allowed_speed * v/fabs(v);
    984             //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    985             //   v, reduced_speed);
    986             ymomc[k] = reduced_speed * hc;
    987           }
     982      if (fabs(v) > maximum_allowed_speed) {
     983        reduced_speed = maximum_allowed_speed * v/fabs(v);
     984        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     985        //   v, reduced_speed);
     986        ymomc[k] = reduced_speed * hc;
     987      }
    988988        }
    989989      }
    … …  
    10561056              // But evidence suggests that h0 can be as little as
    10571057              // epsilon!
    1058              
     1058
    10591059  limiting_threshold = 10*H0; // Avoid applying limiter below this
    10601060                              // threshold for performance reasons.
    … …  
    10621062 
    10631063  _flux_function_central((double*) ql -> data,
    1064                          (double*) qr -> data,
    1065                          zl,
    1066                          zr,                         
    1067                         ((double*) normal -> data)[0],
    1068                          ((double*) normal -> data)[1],         
    1069                         epsilon, h0, limiting_threshold,
    1070                         g,
    1071                          (double*) edgeflux -> data,
    1072                         &max_speed);
     1064             (double*) qr -> data,
     1065             zl,
     1066             zr,
     1067            ((double*) normal -> data)[0],
     1068             ((double*) normal -> data)[1],
     1069            epsilon, h0, limiting_threshold,
     1070            g,
     1071             (double*) edgeflux -> data,
     1072            &max_speed);
    10731073 
    10741074  return Py_BuildValue("d", max_speed); 
    … …  
    11561156
    11571157  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
    1158                         &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
     1158            &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
    11591159      report_python_error(AT, "could not parse input arguments");
    11601160      return NULL;
    … …  
    11741174
    11751175  _manning_friction_sloped(g, eps, N,
    1176                         (double*) x -> data,
    1177                         (double*) w -> data,
    1178                         (double*) z -> data,
    1179                         (double*) uh -> data,
    1180                         (double*) vh -> data,
    1181                         (double*) eta -> data,
    1182                         (double*) xmom -> data,
    1183                         (double*) ymom -> data);
     1176            (double*) x -> data,
     1177            (double*) w -> data,
     1178            (double*) z -> data,
     1179            (double*) uh -> data,
     1180            (double*) vh -> data,
     1181            (double*) eta -> data,
     1182            (double*) xmom -> data,
     1183            (double*) ymom -> data);
    11841184
    11851185  return Py_BuildValue("");
    … …  
    11981198
    11991199  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
    1200                         &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
     1200            &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
    12011201      report_python_error(AT, "could not parse input arguments");
    12021202      return NULL;
    … …  
    12161216
    12171217  _chezy_friction(g, eps, N,
    1218                         (double*) x -> data,
    1219                         (double*) w -> data,
    1220                         (double*) z -> data,
    1221                         (double*) uh -> data,
    1222                         (double*) vh -> data,
    1223                         (double*) chezy -> data,
    1224                         (double*) xmom -> data,
    1225                         (double*) ymom -> data);
     1218            (double*) x -> data,
     1219            (double*) w -> data,
     1220            (double*) z -> data,
     1221            (double*) uh -> data,
     1222            (double*) vh -> data,
     1223            (double*) chezy -> data,
     1224            (double*) xmom -> data,
     1225            (double*) ymom -> data);
    12261226
    12271227  return Py_BuildValue("");
    … …  
    12401240
    12411241  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    1242                         &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
     1242            &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
    12431243      report_python_error(AT, "could not parse input arguments");
    12441244      return NULL;
    … …  
    12571257
    12581258  _manning_friction_flat(g, eps, N,
    1259                         (double*) w -> data,
    1260                         (double*) z -> data,
    1261                         (double*) uh -> data,
    1262                         (double*) vh -> data,
    1263                         (double*) eta -> data,
    1264                         (double*) xmom -> data,
    1265                         (double*) ymom -> data);
     1259            (double*) w -> data,
     1260            (double*) z -> data,
     1261            (double*) uh -> data,
     1262            (double*) vh -> data,
     1263            (double*) eta -> data,
     1264            (double*) xmom -> data,
     1265            (double*) ymom -> data);
    12661266
    12671267  return Py_BuildValue("");
    … …  
    14841484      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
    14851485     
    1486           inv_area2 = 1.0/area2;
     1486      inv_area2 = 1.0/area2;
    14871487      // Calculate the gradient of stage on the auxiliary triangle
    14881488      a = dy2*dq1 - dy1*dq2;
    … …  
    15171517      //for (i=0;i<3;i++)
    15181518      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
    1519           stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
    1520           stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
     1519      stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
     1520      stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
    15211521     
    15221522     

  • trunk/anuga_core/source/anuga/structures/structure_operator.py

    r8095 r8129  
    66from anuga.utilities.system_tools import log_to_file
    77from anuga.utilities.numerical_tools import ensure_numeric
    8 
    9 
    10 
    11 class Structure_operator:
     8from anuga.operators.operator import Operator
     9
     10
     11
     12class Structure_operator(anuga.Operator):
    1213    """Structure Operator - transfer water from one rectangular box to another.
    1314    Sets up the geometry of problem
    … …  
    3738                 logging,
    3839                 verbose):
    39        
    40         self.domain = domain
    41         self.domain.set_fractional_step_operator(self)
     40
     41        anuga.Operator.__init__(self,domain)
     42       
    4243        self.end_points = ensure_numeric(end_points)
    4344        self.exchange_lines = ensure_numeric(exchange_lines)
    4445        self.enquiry_points = ensure_numeric(enquiry_points)
    45 
    46 
    47         if domain.numproc > 1:
    48             msg = 'Not implemented to run in parallel'
    49             assert self.__parallel_safe(), msg
    5046
    5147       
    … …  
    259255            self.enquiry_points.append(centre_point1 + gap)
    260256           
    261 
    262     def __parallel_safe(self):
    263 
    264         return False
    265257
    266258    def discharge_routine(self):

  • trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py

    r8073 r8129  
    33
    44import unittest
    5 import os.path
    6 import sys
    7 
    8 from anuga.utilities.system_tools import get_pathname_from_package
     5
     6
    97from anuga.structures.boyd_box_operator import Boyd_box_operator
    108from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    119from anuga.shallow_water.shallow_water_domain import Domain
    12 from anuga.shallow_water.forcing import Rainfall, Inflow
    1310import numpy
    1411

  • trunk/anuga_core/source/anuga/structures/test_inlet_operator.py

    r8094 r8129  
    155155        finaltime = 3.0
    156156        line1 = [[95.0, 10.0], [105.0, 10.0]]
    157         Q1 = file_function('inlet_operator_test1.tms', quantities=['hydrograph'])
     157
     158       
     159        import os
     160        baseDir = os.getcwd()
     161
     162        try:
     163            os.chdir('structures')
     164        except:
     165            pass
     166
     167        Q1 = file_function(filename='inlet_operator_test1.tms', quantities=['hydrograph'])
    158168       
    159169        line2 = [[10.0, 90.0], [20.0, 90.0]]
    160         Q2 = file_function('inlet_operator_test2.tms', quantities=['hydrograph'])
     170        Q2 = file_function(filename='inlet_operator_test2.tms', quantities=['hydrograph'])
     171
     172        os.chdir(baseDir)
    161173       
    162174        Inlet_operator(domain, line1, Q1)

  • trunk/anuga_core/source/anuga/test_all.py

    r8099 r8129  
    2323
    2424# Directories that should not be searched for test files.
    25 exclude_dirs = ['pypar_dist',  'shallow_water_balanced',  'structures', # Special requirements
     25exclude_dirs = ['shallow_water_balanced', 'operators' , # Special requirements
    2626                '.svn',          # subversion
    2727                'props', 'wcprops', 'prop-base', 'text-base', 'tmp']
Note: See TracChangeset for help on using the changeset viewer.