Ignore:
Timestamp:
Jun 3, 2009, 1:16:10 PM (15 years ago)
Author:
ole
Message:

Added another memset in quantity_ext.c and fixed some indentation issues to
comply with the Emacs c-style.

Location:
anuga_core/source/anuga
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r6654 r7154  
    3030                        double* b){
    3131
    32         int i, k, k0, k1, k2, index3;
    33         double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
    34 
    35 
    36         for (k=0; k<N; k++) {
    37                 index3 = 3*k;
    38 
    39                 if (number_of_boundaries[k] < 2) {
    40                         // Two or three true neighbours
    41 
    42                         // Get indices of neighbours (or self when used as surrogate)
    43                         // k0, k1, k2 = surrogate_neighbours[k,:]
    44 
    45                         k0 = surrogate_neighbours[index3 + 0];
    46                         k1 = surrogate_neighbours[index3 + 1];
    47                         k2 = surrogate_neighbours[index3 + 2];
    48 
    49 
    50                         if (k0 == k1 || k1 == k2) return -1;
    51 
    52                         // Get data
    53                         q0 = centroid_values[k0];
    54                         q1 = centroid_values[k1];
    55                         q2 = centroid_values[k2];
    56 
    57                         x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    58                         x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    59                         x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
    60 
    61                         // Gradient
    62                         _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
    63 
    64                 } else if (number_of_boundaries[k] == 2) {
    65                         // One true neighbour
    66 
    67                         // Get index of the one neighbour
    68                         i=0; k0 = k;
    69                         while (i<3 && k0==k) {
    70                                 k0 = surrogate_neighbours[index3 + i];
    71                                 i++;
    72                         }
    73                         if (k0 == k) return -1;
    74 
    75                         k1 = k; //self
    76 
    77                         // Get data
    78                         q0 = centroid_values[k0];
    79                         q1 = centroid_values[k1];
    80 
    81                         x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
    82                         x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
    83 
    84                         // Two point gradient
    85                         _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
    86 
    87                 }
    88                 //    else:
    89                 //        #No true neighbours -
    90                 //        #Fall back to first order scheme
    91         }
    92         return 0;
     32  int i, k, k0, k1, k2, index3;
     33  double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
     34
     35
     36  for (k=0; k<N; k++) {
     37    index3 = 3*k;
     38   
     39    if (number_of_boundaries[k] < 2) {
     40      // Two or three true neighbours
     41     
     42      // Get indices of neighbours (or self when used as surrogate)
     43      // k0, k1, k2 = surrogate_neighbours[k,:]
     44     
     45      k0 = surrogate_neighbours[index3 + 0];
     46      k1 = surrogate_neighbours[index3 + 1];
     47      k2 = surrogate_neighbours[index3 + 2];
     48
     49     
     50      if (k0 == k1 || k1 == k2) return -1;
     51     
     52      // Get data
     53      q0 = centroid_values[k0];
     54      q1 = centroid_values[k1];
     55      q2 = centroid_values[k2];
     56
     57      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     58      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     59      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
     60
     61      // Gradient
     62      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
     63
     64    } else if (number_of_boundaries[k] == 2) {
     65      // One true neighbour
     66
     67      // Get index of the one neighbour
     68      i=0; k0 = k;
     69      while (i<3 && k0==k) {
     70        k0 = surrogate_neighbours[index3 + i];
     71        i++;
     72      }
     73      if (k0 == k) return -1;
     74     
     75      k1 = k; //self
     76
     77      // Get data
     78      q0 = centroid_values[k0];
     79      q1 = centroid_values[k1];
     80     
     81      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
     82      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
     83     
     84      // Two point gradient
     85      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
     86     
     87    }
     88    //    else:
     89    //        #No true neighbours -
     90    //        #Fall back to first order scheme
     91  }
     92  return 0;
    9393}
    9494
    9595
    9696int _extrapolate_from_gradient(int N,
    97                 double* centroids,
    98                 double* centroid_values,
    99                 double* vertex_coordinates,
    100                 double* vertex_values,
    101                 double* edge_values,
    102                 double* a,
    103                 double* b) {
    104 
    105         int k, k2, k3, k6;
    106         double x, y, x0, y0, x1, y1, x2, y2;
    107 
    108         for (k=0; k<N; k++){
    109                 k6 = 6*k;
    110                 k3 = 3*k;
    111                 k2 = 2*k;
    112 
    113                 // Centroid coordinates
    114                 x = centroids[k2]; y = centroids[k2+1];
    115 
    116                 // vertex coordinates
    117                 // x0, y0, x1, y1, x2, y2 = X[k,:]
    118                 x0 = vertex_coordinates[k6 + 0];
    119                 y0 = vertex_coordinates[k6 + 1];
    120                 x1 = vertex_coordinates[k6 + 2];
    121                 y1 = vertex_coordinates[k6 + 3];
    122                 x2 = vertex_coordinates[k6 + 4];
    123                 y2 = vertex_coordinates[k6 + 5];
    124 
    125                 // Extrapolate to Vertices
    126                 vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
    127                 vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    128                 vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
    129 
    130                 // Extrapolate to Edges (midpoints)
    131                 edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
    132                 edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
    133                 edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    134 
    135         }
    136         return 0;
     97                              double* centroids,
     98                              double* centroid_values,
     99                              double* vertex_coordinates,
     100                              double* vertex_values,
     101                              double* edge_values,
     102                              double* a,
     103                              double* b) {
     104
     105  int k, k2, k3, k6;
     106  double x, y, x0, y0, x1, y1, x2, y2;
     107 
     108  for (k=0; k<N; k++){
     109    k6 = 6*k;
     110    k3 = 3*k;
     111    k2 = 2*k;
     112   
     113    // Centroid coordinates
     114    x = centroids[k2]; y = centroids[k2+1];
     115   
     116    // vertex coordinates
     117    // x0, y0, x1, y1, x2, y2 = X[k,:]
     118    x0 = vertex_coordinates[k6 + 0];
     119    y0 = vertex_coordinates[k6 + 1];
     120    x1 = vertex_coordinates[k6 + 2];
     121    y1 = vertex_coordinates[k6 + 3];
     122    x2 = vertex_coordinates[k6 + 4];
     123    y2 = vertex_coordinates[k6 + 5];
     124   
     125    // Extrapolate to Vertices
     126    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
     127    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
     128    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     129   
     130    // Extrapolate to Edges (midpoints)
     131    edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
     132    edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
     133    edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
     134   
     135  }
     136  return 0;
    137137}
    138138
    … …  
    149149                                         double* y_gradient) {
    150150
    151         int i, k, k2, k3, k6;
    152         double x, y, x0, y0, x1, y1, x2, y2;
    153         long n;
    154         double qmin, qmax, qc;
    155         double qn[3];
    156         double dq, dqa[3], r;
    157 
    158         for (k=0; k<N; k++){
    159                 k6 = 6*k;
    160                 k3 = 3*k;
    161                 k2 = 2*k;
    162 
    163                 // Centroid coordinates
    164                 x = centroids[k2+0];
    165                 y = centroids[k2+1];
    166 
    167                 // vertex coordinates
    168                 // x0, y0, x1, y1, x2, y2 = X[k,:]
    169                 x0 = vertex_coordinates[k6 + 0];
    170                 y0 = vertex_coordinates[k6 + 1];
    171                 x1 = vertex_coordinates[k6 + 2];
    172                 y1 = vertex_coordinates[k6 + 3];
    173                 x2 = vertex_coordinates[k6 + 4];
    174                 y2 = vertex_coordinates[k6 + 5];
    175 
    176                 // Extrapolate to Vertices
    177                 vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
    178                 vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
    179                 vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
    180 
    181                 // Extrapolate to Edges (midpoints)
    182                 edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
    183                 edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
    184                 edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    185         }
    186 
    187 
    188 
    189         for (k=0; k<N; k++){
    190                 k6 = 6*k;
    191                 k3 = 3*k;
    192                 k2 = 2*k;
    193 
    194 
    195                 qc = centroid_values[k];
    196                
    197                 qmin = qc;
    198                 qmax = qc;
    199                
    200                 for (i=0; i<3; i++) {
    201                     n = neighbours[k3+i];
    202                     if (n < 0) {
    203                         qn[i] = qc;
    204                     } else {
    205                         qn[i] = centroid_values[n];
    206                     }
    207 
    208                     qmin = min(qmin, qn[i]);
    209                     qmax = max(qmax, qn[i]);
    210                 }
    211 
    212                 //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
    213                 //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
    214 
    215 /*              for (i=0; i<3; i++) { */
    216 /*                  n = neighbours[k3+i]; */
    217 /*                  if (n < 0) { */
    218 /*                      qn[i] = qc; */
    219 /*                      qmin[i] = qtmin; */
    220 /*                      qmax[i] = qtmax; */
    221 /*                  }  */
    222 /*              } */
    223                        
    224                 phi[k] = 1.0;
    225 
    226                 for (i=0; i<3; i++) {
    227                     dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
    228                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    229 
    230                     r = 1.0;
    231      
    232                     if (dq > 0.0) r = (qmax - qc)/dq;
    233                     if (dq < 0.0) r = (qmin - qc)/dq;
    234                    
    235                     phi[k] = min( min(r*beta, 1.0), phi[k]);
    236                    
    237                 }
    238 
    239 
    240 
    241                 //Update gradient, edge and vertex values using phi limiter
    242                 x_gradient[k] = x_gradient[k]*phi[k];
    243                 y_gradient[k] = y_gradient[k]*phi[k];
    244 
    245                 edge_values[k3+0] = qc + phi[k]*dqa[0];
    246                 edge_values[k3+1] = qc + phi[k]*dqa[1];
    247                 edge_values[k3+2] = qc + phi[k]*dqa[2];
    248 
    249 
    250                 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
    251                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    252                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
    253                
    254 
    255         }
    256 
    257         return 0;
    258 
     151  int i, k, k2, k3, k6;
     152  double x, y, x0, y0, x1, y1, x2, y2;
     153  long n;
     154  double qmin, qmax, qc;
     155  double qn[3];
     156  double dq, dqa[3], r;
     157
     158  for (k=0; k<N; k++){
     159    k6 = 6*k;
     160    k3 = 3*k;
     161    k2 = 2*k;
     162
     163    // Centroid coordinates
     164    x = centroids[k2+0];
     165    y = centroids[k2+1];
     166
     167    // vertex coordinates
     168    // x0, y0, x1, y1, x2, y2 = X[k,:]
     169    x0 = vertex_coordinates[k6 + 0];
     170    y0 = vertex_coordinates[k6 + 1];
     171    x1 = vertex_coordinates[k6 + 2];
     172    y1 = vertex_coordinates[k6 + 3];
     173    x2 = vertex_coordinates[k6 + 4];
     174    y2 = vertex_coordinates[k6 + 5];
     175   
     176    // Extrapolate to Vertices
     177    vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
     178    vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
     179    vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
     180   
     181    // Extrapolate to Edges (midpoints)
     182    edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
     183    edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
     184    edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
     185  }
     186 
     187 
     188 
     189  for (k=0; k<N; k++){
     190    k6 = 6*k;
     191    k3 = 3*k;
     192    k2 = 2*k;
     193   
     194   
     195    qc = centroid_values[k];
     196   
     197    qmin = qc;
     198    qmax = qc;
     199   
     200    for (i=0; i<3; i++) {
     201      n = neighbours[k3+i];
     202      if (n < 0) {
     203        qn[i] = qc;
     204      } else {
     205        qn[i] = centroid_values[n];
     206      }
     207     
     208      qmin = min(qmin, qn[i]);
     209      qmax = max(qmax, qn[i]);
     210    }
     211   
     212    //qtmin = min(min(min(qn[0],qn[1]),qn[2]),qc);
     213    //qtmax = max(max(max(qn[0],qn[1]),qn[2]),qc);
     214   
     215    /*          for (i=0; i<3; i++) { */
     216    /*              n = neighbours[k3+i]; */
     217    /*              if (n < 0) { */
     218    /*                  qn[i] = qc; */
     219    /*                  qmin[i] = qtmin; */
     220    /*                  qmax[i] = qtmax; */
     221    /*              }  */
     222    /*          } */
     223   
     224    phi[k] = 1.0;
     225   
     226    for (i=0; i<3; i++) {
     227      dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
     228      dqa[i] = dq;                      //Save dq for use in updating vertex values
     229     
     230      r = 1.0;
     231     
     232      if (dq > 0.0) r = (qmax - qc)/dq;
     233      if (dq < 0.0) r = (qmin - qc)/dq;
     234     
     235      phi[k] = min( min(r*beta, 1.0), phi[k]);
     236     
     237    }
     238   
     239   
     240   
     241    //Update gradient, edge and vertex values using phi limiter
     242    x_gradient[k] = x_gradient[k]*phi[k];
     243    y_gradient[k] = y_gradient[k]*phi[k];
     244   
     245    edge_values[k3+0] = qc + phi[k]*dqa[0];
     246    edge_values[k3+1] = qc + phi[k]*dqa[1];
     247    edge_values[k3+2] = qc + phi[k]*dqa[2];
     248   
     249   
     250    vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     251    vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     252    vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     253   
     254   
     255  }
     256 
     257  return 0;
     258 
    259259}
    260260
    … …  
    263263
    264264int _limit_vertices_by_all_neighbours(int N, double beta,
    265                      double* centroid_values,
    266                      double* vertex_values,
    267                      double* edge_values,
    268                      long*   neighbours,
    269                      double* x_gradient,
    270                      double* y_gradient) {
    271 
    272 
    273         int i, k, k2, k3, k6;
    274         long n;
    275         double qmin, qmax, qn, qc;
    276         double dq, dqa[3], phi, r;
    277 
    278         for (k=0; k<N; k++){
    279                 k6 = 6*k;
    280                 k3 = 3*k;
    281                 k2 = 2*k;
    282 
    283                 qc = centroid_values[k];
    284                 qmin = qc;
    285                 qmax = qc;
    286 
    287                 for (i=0; i<3; i++) {
    288                     n = neighbours[k3+i];
    289                     if (n >= 0) {
    290                         qn = centroid_values[n]; //Neighbour's centroid value
    291 
    292                         qmin = min(qmin, qn);
    293                         qmax = max(qmax, qn);
    294                     }
    295                 }
    296 
    297                 phi = 1.0;
    298                 for (i=0; i<3; i++) {   
    299                     r = 1.0;
    300      
    301                     dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
    302                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    303      
    304                     if (dq > 0.0) r = (qmax - qc)/dq;
    305                     if (dq < 0.0) r = (qmin - qc)/dq;     
     265                                      double* centroid_values,
     266                                      double* vertex_values,
     267                                      double* edge_values,
     268                                      long*   neighbours,
     269                                      double* x_gradient,
     270                                      double* y_gradient) {
    306271 
    307                    
    308                     phi = min( min(r*beta, 1.0), phi);   
    309                 }
    310 
    311                 //Update gradient, vertex and edge values using phi limiter
    312                 x_gradient[k] = x_gradient[k]*phi;
    313                 y_gradient[k] = y_gradient[k]*phi;
    314    
    315                 vertex_values[k3+0] = qc + phi*dqa[0];
    316                 vertex_values[k3+1] = qc + phi*dqa[1];
    317                 vertex_values[k3+2] = qc + phi*dqa[2];
    318                
    319                 edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
    320                 edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
    321                 edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
    322 
    323         }
    324 
    325         return 0;
     272 
     273  int i, k, k2, k3, k6;
     274  long n;
     275  double qmin, qmax, qn, qc;
     276  double dq, dqa[3], phi, r;
     277 
     278  for (k=0; k<N; k++){
     279    k6 = 6*k;
     280    k3 = 3*k;
     281    k2 = 2*k;
     282   
     283    qc = centroid_values[k];
     284    qmin = qc;
     285    qmax = qc;
     286   
     287    for (i=0; i<3; i++) {
     288      n = neighbours[k3+i];
     289      if (n >= 0) {
     290        qn = centroid_values[n]; //Neighbour's centroid value
     291       
     292        qmin = min(qmin, qn);
     293        qmax = max(qmax, qn);
     294      }
     295    }
     296   
     297    phi = 1.0;
     298    for (i=0; i<3; i++) {   
     299      r = 1.0;
     300     
     301      dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
     302      dqa[i] = dq;                      //Save dq for use in updating vertex values
     303     
     304      if (dq > 0.0) r = (qmax - qc)/dq;
     305      if (dq < 0.0) r = (qmin - qc)/dq;     
     306     
     307     
     308      phi = min( min(r*beta, 1.0), phi);   
     309    }
     310   
     311    //Update gradient, vertex and edge values using phi limiter
     312    x_gradient[k] = x_gradient[k]*phi;
     313    y_gradient[k] = y_gradient[k]*phi;
     314   
     315    vertex_values[k3+0] = qc + phi*dqa[0];
     316    vertex_values[k3+1] = qc + phi*dqa[1];
     317    vertex_values[k3+2] = qc + phi*dqa[2];
     318   
     319    edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]);
     320    edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]);
     321    edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]);
     322   
     323  }
     324 
     325  return 0;
    326326}
    327327
    … …  
    337337                                   double* y_gradient) {
    338338
    339         int i, k, k2, k3, k6;
    340         long n;
    341         double qmin, qmax, qn, qc;
    342         double dq, dqa[3], phi, r;
    343 
    344         for (k=0; k<N; k++){
    345                 k6 = 6*k;
    346                 k3 = 3*k;
    347                 k2 = 2*k;
    348 
    349                 qc = centroid_values[k];
    350                 qmin = qc;
    351                 qmax = qc;
    352 
    353                 for (i=0; i<3; i++) {
    354                     n = neighbours[k3+i];
    355                     if (n >= 0) {
    356                         qn = centroid_values[n]; //Neighbour's centroid value
    357 
    358                         qmin = min(qmin, qn);
    359                         qmax = max(qmax, qn);
    360                     }
    361                 }
    362 
    363                 phi = 1.0;
    364                 for (i=0; i<3; i++) {   
    365                     r = 1.0;
    366      
    367                     dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
    368                     dqa[i] = dq;                      //Save dq for use in updating vertex values
    369      
    370                     if (dq > 0.0) r = (qmax - qc)/dq;
    371                     if (dq < 0.0) r = (qmin - qc)/dq;     
     339  int i, k, k2, k3, k6;
     340  long n;
     341  double qmin, qmax, qn, qc;
     342  double dq, dqa[3], phi, r;
    372343 
    373                    
    374                     phi = min( min(r*beta, 1.0), phi);   
    375                 }
    376    
    377                 //Update gradient, vertex and edge values using phi limiter
    378                 x_gradient[k] = x_gradient[k]*phi;
    379                 y_gradient[k] = y_gradient[k]*phi;
    380 
    381                 edge_values[k3+0] = qc + phi*dqa[0];
    382                 edge_values[k3+1] = qc + phi*dqa[1];
    383                 edge_values[k3+2] = qc + phi*dqa[2];
    384                
    385                 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
    386                 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
    387                 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
    388 
    389         }
    390 
    391         return 0;
     344  for (k=0; k<N; k++){
     345    k6 = 6*k;
     346    k3 = 3*k;
     347    k2 = 2*k;
     348   
     349    qc = centroid_values[k];
     350    qmin = qc;
     351    qmax = qc;
     352   
     353    for (i=0; i<3; i++) {
     354      n = neighbours[k3+i];
     355      if (n >= 0) {
     356        qn = centroid_values[n]; //Neighbour's centroid value
     357       
     358        qmin = min(qmin, qn);
     359        qmax = max(qmax, qn);
     360      }
     361    }
     362   
     363    phi = 1.0;
     364    for (i=0; i<3; i++) {   
     365      r = 1.0;
     366     
     367      dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
     368      dqa[i] = dq;                      //Save dq for use in updating vertex values
     369     
     370      if (dq > 0.0) r = (qmax - qc)/dq;
     371      if (dq < 0.0) r = (qmin - qc)/dq;     
     372     
     373     
     374      phi = min( min(r*beta, 1.0), phi);   
     375    }
     376   
     377    //Update gradient, vertex and edge values using phi limiter
     378    x_gradient[k] = x_gradient[k]*phi;
     379    y_gradient[k] = y_gradient[k]*phi;
     380   
     381    edge_values[k3+0] = qc + phi*dqa[0];
     382    edge_values[k3+1] = qc + phi*dqa[1];
     383    edge_values[k3+2] = qc + phi*dqa[2];
     384   
     385    vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     386    vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     387    vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     388   
     389  }
     390 
     391  return 0;
    392392}
    393393
    … …  
    730730
    731731
    732         // MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
    733         for (k=0;k<N;k++){
    734                 semi_implicit_update[k]=0.0;
    735         }
     732        // Reset semi_implicit_update here ready for next time step
     733        memset(semi_implicit_update, 0, N*sizeof(double));
    736734
    737735        return 0;

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

    r7143 r7154  
    6868    if (dq1>=dq2){
    6969      if (dq1>=0.0)
    70     *qmax=dq0+dq1;
     70        *qmax=dq0+dq1;
    7171      else
    72     *qmax=dq0;
    73    
     72        *qmax=dq0;
     73     
    7474      *qmin=dq0+dq2;
    7575      if (*qmin>=0.0) *qmin = 0.0;
    … …  
    7777    else{// dq1<dq2
    7878      if (dq2>0)
    79     *qmax=dq0+dq2;
     79        *qmax=dq0+dq2;
    8080      else
    81     *qmax=dq0;
    82    
     81        *qmax=dq0;
     82     
    8383      *qmin=dq0+dq1;   
    8484      if (*qmin>=0.0) *qmin=0.0;
    … …  
    8888    if (dq1<=dq2){
    8989      if (dq1<0.0)
    90     *qmin=dq0+dq1;
     90        *qmin=dq0+dq1;
    9191      else
    92     *qmin=dq0;
    93    
     92        *qmin=dq0;
     93     
    9494      *qmax=dq0+dq2;   
    9595      if (*qmax<=0.0) *qmax=0.0;
    … …  
    9797    else{// dq1>dq2
    9898      if (dq2<0.0)
    99     *qmin=dq0+dq2;
     99        *qmin=dq0+dq2;
    100100      else
    101     *qmin=dq0;
    102    
     101        *qmin=dq0;
     102     
    103103      *qmax=dq0+dq1;
    104104      if (*qmax<=0.0) *qmax=0.0;
    … …  
    683683      // FIXME: Try with this one precomputed
    684684      for (i=0; i<3; i++) {
    685     dz = max(dz, fabs(zv[k3+i]-zc[k]));
     685        dz = max(dz, fabs(zv[k3+i]-zc[k]));
    686686      }
    687687    }
    … …  
    711711     
    712712      if (dz > 0.0) {
    713     alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     713        alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
    714714      } else {
    715     alpha = 1.0;  // Flat bed
     715        alpha = 1.0;  // Flat bed
    716716      }
    717717      //printf("Using old style limiter\n");
    … …  
    726726   
    727727      if (hmin < H0) {
    728     alpha = 1.0;
    729     for (i=0; i<3; i++) {
    730      
    731       h_diff = hc_k - hv[i];     
    732       if (h_diff <= 0) {
    733         // Deep water triangle is further away from bed than
    734         // shallow water (hbar < h). Any alpha will do
    735      
    736       } else { 
    737         // Denominator is positive which means that we need some of the
    738         // h-limited stage.
    739        
    740         alpha = min(alpha, (hc_k - H0)/h_diff);     
     728        alpha = 1.0;
     729        for (i=0; i<3; i++) {
     730         
     731          h_diff = hc_k - hv[i];     
     732          if (h_diff <= 0) {
     733            // Deep water triangle is further away from bed than
     734            // shallow water (hbar < h). Any alpha will do
     735     
     736          } else { 
     737            // Denominator is positive which means that we need some of the
     738            // h-limited stage.
     739           
     740            alpha = min(alpha, (hc_k - H0)/h_diff);     
     741          }
     742        }
     743
     744        // Ensure alpha in [0,1]
     745        if (alpha>1.0) alpha=1.0;
     746        if (alpha<0.0) alpha=0.0;
     747       
     748      } else {
     749        // Use w-limited stage exclusively in deeper water.
     750        alpha = 1.0;       
    741751      }
    742752    }
    743 
    744     // Ensure alpha in [0,1]
    745     if (alpha>1.0) alpha=1.0;
    746     if (alpha<0.0) alpha=0.0;
    747    
    748       } else {
    749     // Use w-limited stage exclusively in deeper water.
    750     alpha = 1.0;       
    751       }
    752     }
    753 
    754 
     753   
     754   
    755755    //  Let
    756756    //
    … …  
    770770    //   Momentum is balanced between constant and limited
    771771
    772 
     772   
    773773    if (alpha < 1) {     
    774774      for (i=0; i<3; i++) {
    775      
    776     wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
    777 
    778     // Update momentum at vertices
    779     if (use_centroid_velocities == 1) {
    780       // This is a simple, efficient and robust option
    781       // It uses first order approximation of velocities, but retains
    782       // the order used by stage.
    783    
    784       // Speeds at centroids
    785       if (hc_k > epsilon) {
    786         uc = xmomc[k]/hc_k;
    787         vc = ymomc[k]/hc_k;
    788       } else {
    789         uc = 0.0;
    790         vc = 0.0;
    791       }
    792      
    793       // Vertex momenta guaranteed to be consistent with depth guaranteeing
    794       // controlled speed
    795       hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    796       xmomv[k3+i] = uc*hv[i];
    797       ymomv[k3+i] = vc*hv[i];
    798      
    799     } else {
    800       // Update momentum as a linear combination of
    801       // xmomc and ymomc (shallow) and momentum
    802       // from extrapolator xmomv and ymomv (deep).
    803       // This assumes that values from xmomv and ymomv have
    804       // been established e.g. by the gradient limiter.
    805 
    806       // FIXME (Ole): I think this should be used with vertex momenta
    807       // computed above using centroid_velocities instead of xmomc
    808       // and ymomc as they'll be more representative first order
    809       // values.
    810      
    811       xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    812       ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    813    
    814     }
     775       
     776        wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     777
     778        // Update momentum at vertices
     779        if (use_centroid_velocities == 1) {
     780          // This is a simple, efficient and robust option
     781          // It uses first order approximation of velocities, but retains
     782          // the order used by stage.
     783   
     784          // Speeds at centroids
     785          if (hc_k > epsilon) {
     786            uc = xmomc[k]/hc_k;
     787            vc = ymomc[k]/hc_k;
     788          } else {
     789            uc = 0.0;
     790            vc = 0.0;
     791          }
     792         
     793          // Vertex momenta guaranteed to be consistent with depth guaranteeing
     794          // controlled speed
     795          hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     796          xmomv[k3+i] = uc*hv[i];
     797          ymomv[k3+i] = vc*hv[i];
     798     
     799        } else {
     800          // Update momentum as a linear combination of
     801          // xmomc and ymomc (shallow) and momentum
     802          // from extrapolator xmomv and ymomv (deep).
     803          // This assumes that values from xmomv and ymomv have
     804          // been established e.g. by the gradient limiter.
     805         
     806          // FIXME (Ole): I think this should be used with vertex momenta
     807          // computed above using centroid_velocities instead of xmomc
     808          // and ymomc as they'll be more representative first order
     809          // values.
     810         
     811          xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     812          ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     813         
     814        }
    815815      }
    816816    }
    … …  
    843843      if (hc < minimum_allowed_height) {
    844844       
    845     // Set momentum to zero and ensure h is non negative
    846     xmomc[k] = 0.0;
    847     ymomc[k] = 0.0;
    848     if (hc <= 0.0) wc[k] = zc[k];
     845        // Set momentum to zero and ensure h is non negative
     846        xmomc[k] = 0.0;
     847        ymomc[k] = 0.0;
     848        if (hc <= 0.0) wc[k] = zc[k];
    849849      }
    850850    }
    … …  
    854854    for (k=0; k<N; k++) {
    855855      hc = wc[k] - zc[k];
    856 
     856     
    857857      if (hc < minimum_allowed_height) {
    858858
    … …  
    866866        } else {
    867867          //Reduce excessive speeds derived from division by small hc
    868         //FIXME (Ole): This may be unnecessary with new slope limiters
    869         //in effect.
     868          //FIXME (Ole): This may be unnecessary with new slope limiters
     869          //in effect.
    870870         
    871871          u = xmomc[k]/hc;
    872       if (fabs(u) > maximum_allowed_speed) {
    873         reduced_speed = maximum_allowed_speed * u/fabs(u);
    874         //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    875         //   u, reduced_speed);
    876         xmomc[k] = reduced_speed * hc;
    877       }
    878 
     872          if (fabs(u) > maximum_allowed_speed) {
     873            reduced_speed = maximum_allowed_speed * u/fabs(u);
     874            //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     875            //   u, reduced_speed);
     876            xmomc[k] = reduced_speed * hc;
     877          }
     878         
    879879          v = ymomc[k]/hc;
    880       if (fabs(v) > maximum_allowed_speed) {
    881         reduced_speed = maximum_allowed_speed * v/fabs(v);
    882         //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    883         //   v, reduced_speed);
    884         ymomc[k] = reduced_speed * hc;
    885       }
     880          if (fabs(v) > maximum_allowed_speed) {
     881            reduced_speed = maximum_allowed_speed * v/fabs(v);
     882            //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     883            //   v, reduced_speed);
     884            ymomc[k] = reduced_speed * hc;
     885          }
    886886        }
    887887      }
Note: See TracChangeset for help on using the changeset viewer.