Changeset 6140


Ignore:
Timestamp:
Jan 13, 2009, 8:58:24 AM (16 years ago)
Author:
wilson
Message:

version 1.0

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/channel_domain_ext.c

    r6038 r6140  
    1414/*      z=(a>b)?a:b; */
    1515/*      return z;} */
    16        
     16
    1717/* double min(double a, double b) { */
    1818/*      double z; */
    … …  
    2424// This is used by flux functions
    2525// Input parameters uh and h may be modified by this function.
    26 double _compute_speed(double *uh, 
    27                       double *h, 
    28                       double epsilon, 
     26double _compute_speed(double *uh,
     27                      double *h,
     28                      double epsilon,
    2929                      double h0) {
    30  
     30
    3131  double u;
    32  
     32
    3333  if (*h < epsilon) {
    3434    *h = 0.0;  //Could have been negative
    3535     u = 0.0;
    3636  } else {
    37     u = *uh/(*h + h0/ *h);   
     37    u = *uh/(*h + h0/ *h);
    3838  }
    39  
     39
    4040
    4141  // Adjust momentum to be consistent with speed
    4242  *uh = u * *h;
    43  
     43
    4444  return u;
    4545}
    4646
    47 //-------------------------------------------------------------
    48 // New vel based code
    49 //-------------------------------------------------------------
     47
     48
    5049//Innermost flux function (using w=z+h)
    51 int _flux_function_vel(double *q_left, double *q_right,
    52                        double normals, double g, double epsilon, double h0,
    53                        double *edgeflux, double *max_speed) {
    54    
     50int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm,
     51                           double *q_rightp, double g, double epsilon, double h0,                           double *edgeflux, double *max_speed) {
     52
    5553    int i;
    5654    double flux_left[2], flux_right[2];
    57     double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
    58     double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
    59     double z, s_max, s_min, denom;
    60    
    61 
    62     w_left  = q_left[0];
    63     uh_left = q_left[1]*normals;
    64     z_left  = q_left[2];
    65     h_left  = q_left[3];
    66     u_left  = q_left[4]*normals;
    67 
    68 /*     printf("w_left  = %f \n",w_left); */
    69 /*     printf("uh_left = %f \n",uh_left); */
    70 /*     printf("z_left  = %f \n",z_left); */
    71 /*     printf("h_left  = %f \n",h_left); */
    72 /*     printf("u_left  = %f \n",u_left); */
    73    
    74     w_right  = q_right[0];
    75     uh_right = q_right[1]*normals;
    76     z_right  = q_right[2];
    77     h_right  = q_right[3];
    78     u_right  = q_right[4]*normals;             
    79    
    80     z = (z_left+z_right)/2.0;             
    81    
    82     soundspeed_left = sqrt(g*h_left);
    83     soundspeed_right = sqrt(g*h_right);
    84    
    85     s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
    86     if (s_max < 0.0) s_max = 0.0;
    87    
    88     s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    89     if (s_min > 0.0) s_min = 0.0;
    90    
    91    
    92     // Flux formulas
    93     flux_left[0] = u_left*h_left;
    94     flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
    95    
    96     flux_right[0] = u_right*h_right;
    97     flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
    98    
    99     // Flux computation
    100     denom = s_max-s_min;
     55    double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm;
     56     double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp;
     57 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm;
     58 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp;
     59    double s_maxl, s_minl,s_maxr,s_minr, denom;
     60    double zphalf,zmhalf,hleftstar,hrightstar;
     61    double fluxtemp1,fluxtemp0,speedtemp;
     62
     63    zmhalf = max(q_leftm[2],q_leftp[2]);
     64    zphalf = max(q_rightm[2],q_rightp[2]);   
     65
     66    a_leftm  = q_leftm[0];
     67    d_leftm = q_leftm[1];
     68    z_leftm  = q_leftm[2];
     69    h_leftm  = max(0,q_leftm[3]+q_leftm[2]-zmhalf);
     70    u_leftm  = q_leftm[4];
     71    b_leftm  = q_leftm[5];
     72    w_leftm  = h_leftm+z_leftm;
     73
     74    a_leftp  = q_leftp[0];
     75    d_leftp = q_leftp[1];
     76    z_leftp  = q_leftp[2];
     77    h_leftp  = max(0,q_leftp[3]+q_leftp[2]-zmhalf);
     78    u_leftp  = q_leftp[4];
     79    b_leftp  = q_leftp[5];
     80    w_leftp  = h_leftp+z_leftp;
     81
     82    a_rightm  = q_rightm[0];
     83    d_rightm = q_rightm[1];
     84    z_rightm  = q_rightm[2];
     85    h_rightm  = max(0,q_rightm[3]+q_rightm[2]-zphalf);
     86    u_rightm  = q_rightm[4];
     87    b_rightm  = q_rightm[5];
     88    w_rightm  = h_rightm+z_rightm;
     89
     90    a_rightp  = q_rightp[0];
     91    d_rightp = q_rightp[1];
     92    z_rightp  = q_rightp[2];
     93    h_rightp  = max(0,q_rightp[3]+q_rightp[2]-zphalf);
     94    u_rightp  = q_rightp[4];
     95    b_rightp  = q_rightp[5];
     96    w_rightp  = h_rightp+z_rightp;
     97
     98    hleftstar = q_leftp[3];
     99    hrightstar = q_rightm[3];   
     100
     101
     102
     103
     104    soundspeed_leftp = sqrt(g*h_leftp);
     105    soundspeed_leftm = sqrt(g*h_leftm);
     106    soundspeed_rightp = sqrt(g*h_rightp);
     107    soundspeed_rightm = sqrt(g*h_rightm);
     108
     109    s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp);
     110    if (s_maxl < 0.0) s_maxl = 0.0;
     111
     112    s_minl = min(u_leftm-soundspeed_leftm, u_leftp-soundspeed_leftp);
     113    if (s_minl > 0.0) s_minl = 0.0;
     114
     115    s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp);
     116    if (s_maxr < 0.0) s_maxr = 0.0;
     117
     118    s_minr = min(u_rightm-soundspeed_rightm, u_rightp-soundspeed_rightp);
     119    if (s_minr > 0.0) s_minr = 0.0;
     120
     121    // Flux formulas for left hand side
     122    flux_left[0] = d_leftm;
     123    flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm;
     124
     125    flux_right[0] = d_leftp;
     126    flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp;
     127
     128   
     129    // Flux computation for left hand side
     130    denom = s_maxl-s_minl;
    101131    if (denom < epsilon) {
    102132        for (i=0; i<2; i++) edgeflux[i] = 0.0;
    103         *max_speed = 0.0;
    104133    } else {
    105         edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
    106         edgeflux[0] += s_max*s_min*(w_right-w_left);
     134        edgeflux[0] = s_maxl*flux_left[0] - s_minl*flux_right[0];
     135//      edgeflux[0] += s_maxl*s_minl*(a_leftp-a_leftm);
    107136        edgeflux[0] /= denom;
    108         edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
    109         edgeflux[1] += s_max*s_min*(uh_right-uh_left);
     137        edgeflux[1] = s_maxl*flux_left[1] - s_minl*flux_right[1];
     138        edgeflux[1] += s_maxl*s_minl*(d_leftp-d_leftm);
    110139        edgeflux[1] /= denom;
    111         edgeflux[1] *= normals;
    112        
    113         // Maximal wavespeed
    114         *max_speed = max(fabs(s_max), fabs(s_min));
    115     } 
    116     return 0;           
     140   
     141
     142    }
     143   
     144        fluxtemp0 = edgeflux[0];
     145        fluxtemp1 = edgeflux[1];
     146   
     147
     148    // Flux formulas for right hand side
     149    flux_left[0] = d_rightm;
     150    flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm;
     151
     152    flux_right[0] = d_rightp;
     153    flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp;
     154   
     155   
     156   
     157    // Flux computation for right hand side
     158    denom = s_maxr-s_minr;
     159    if (denom < epsilon) {
     160        for (i=0; i<2; i++) edgeflux[i] = 0.0;
     161    } else {
     162        edgeflux[0] = s_maxr*flux_left[0] - s_minr*flux_right[0];
     163//              edgeflux[0] += s_maxr*s_minr*(a_rightp-a_rightm);
     164        edgeflux[0] /= denom;
     165        edgeflux[1] = s_maxr*flux_left[1] - s_minr*flux_right[1];
     166        edgeflux[1] += s_maxr*s_minr*(d_rightp-d_rightm);
     167        edgeflux[1] /= denom;
     168   
     169
     170    }
     171   
     172   
     173    edgeflux[0]=edgeflux[0]-fluxtemp0;
     174    edgeflux[1]=edgeflux[1]-fluxtemp1;
     175   
     176   
     177     
     178    edgeflux[1]-=0.5*g*h_rightm*h_rightm-0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar-0.5*g*h_leftp*h_leftp;
     179
     180   
     181
     182     
     183
     184     //edgeflux[1]-=0.5*g*b_rightm*h_rightm*h_rightm-0.5*g*b_leftp*h_leftp*h_leftp;
     185                // Maximal wavespeed
     186        if ( (s_maxl-s_minl)<epsilon && (s_maxr-s_minr)<epsilon ){
     187            *max_speed = 0.0;
     188        }else{
     189        speedtemp = max(fabs(s_maxl),fabs(s_minl));
     190        speedtemp = max(speedtemp,fabs(s_maxr));
     191        speedtemp = max(speedtemp,fabs(s_minr));
     192        *max_speed = speedtemp;
     193        }
     194
     195//printf("%f\n",h_right);
     196    return 0;
    117197}
    118198
     199
     200
     201
     202
    119203// Computational function for flux computation
    120 double _compute_fluxes_vel_ext(double cfl,
     204double _compute_fluxes_channel_ext(double cfl,
    121205                               double timestep,
    122206                               double epsilon,
    … …  
    127211                               double* normals,
    128212                               double* areas,
    129                                double* stage_edge_values,
    130                                double* xmom_edge_values,
     213                               double* area_edge_values,
     214                               double* discharge_edge_values,
    131215                               double* bed_edge_values,
    132216                               double* height_edge_values,
    133217                               double* velocity_edge_values,
    134                                double* stage_boundary_values,
    135                                double* xmom_boundary_values,
     218                               double* width_edge_values,
     219                               double* area_boundary_values,
     220                               double* discharge_boundary_values,
    136221                               double* bed_boundary_values,
    137222                               double* height_boundary_values,
    138223                               double* velocity_boundary_values,
    139                                double* stage_explicit_update,
    140                                double* xmom_explicit_update,
     224                               double* width_boundary_values,
     225                               double* area_explicit_update,
     226                               double* discharge_explicit_update,
    141227                               int number_of_elements,
    142228                               double* max_speed_array) {
    143                
    144     double flux[2], ql[5], qr[5], edgeflux[2];
     229
     230    double flux[2], qlm[6],qlp[6], qrm[6],qrp[6], edgeflux[2];
    145231    double max_speed, normal;
    146232    int k, i, ki, n, m, nm=0;
    147    
    148    
     233    double zstar;
    149234    for (k=0; k<number_of_elements; k++) {
    150235        flux[0] = 0.0;
    151236        flux[1] = 0.0;
    152        
    153         for (i=0; i<2; i++) {
    154             ki = k*2+i;
    155            
    156             ql[0] = stage_edge_values[ki];
    157             ql[1] = xmom_edge_values[ki];
    158             ql[2] = bed_edge_values[ki];
    159             ql[3] = height_edge_values[ki];
    160             ql[4] = velocity_edge_values[ki];
    161            
    162             n = neighbours[ki];
     237
     238       
     239            ki = k*2;
     240           
     241
     242         n = neighbours[ki];
    163243            if (n<0) {
    164244                m = -n-1;
    165                 qr[0] = stage_boundary_values[m];
    166                 qr[1] = xmom_boundary_values[m];
    167                 qr[2] = bed_boundary_values[m];
    168                 qr[3] = height_boundary_values[m];
    169                 qr[4] = velocity_boundary_values[m];
    170             } else {
    171                 m = neighbour_vertices[ki];
     245
     246            qlm[0] = area_boundary_values[m];
     247            qlm[1] = discharge_boundary_values[m];
     248            qlm[2] = bed_boundary_values[m];
     249            qlm[3] = height_boundary_values[m];
     250            qlm[4] = velocity_boundary_values[m];
     251            qlm[5] = width_boundary_values[m];
     252
     253            }else{
     254        m = neighbour_vertices[ki];
    172255                nm = n*2+m;
    173                 qr[0] = stage_edge_values[nm];
    174                 qr[1] = xmom_edge_values[nm];
    175                 qr[2] = bed_edge_values[nm];
    176                 qr[3] = height_edge_values[nm];
    177                 qr[4] = velocity_edge_values[nm];
    178             }
    179256           
    180             normal = normals[ki];
    181             _flux_function_vel(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
     257
     258            qlm[0] = area_edge_values[nm];
     259            qlm[1] = discharge_edge_values[nm];
     260            qlm[2] = bed_edge_values[nm];
     261            qlm[3] = height_edge_values[nm];
     262            qlm[4] = velocity_edge_values[nm];
     263            qlm[5] = width_edge_values[nm];
     264    }
     265            qlp[0] = area_edge_values[ki];
     266            qlp[1] = discharge_edge_values[ki];
     267            qlp[2] = bed_edge_values[ki];
     268            qlp[3] = height_edge_values[ki];
     269            qlp[4] = velocity_edge_values[ki];
     270            qlp[5] = width_edge_values[ki];
     271
     272         ki = k*2+1;
     273           
     274
     275         n = neighbours[ki];
     276            if (n<0) {
     277                m = -n-1;
     278            qrp[0] = area_boundary_values[m];
     279            qrp[1] = discharge_boundary_values[m];
     280            qrp[2] = bed_boundary_values[m];
     281            qrp[3] = height_boundary_values[m];
     282            qrp[4] = velocity_boundary_values[m];
     283            qrp[5] = width_boundary_values[m];
     284         
     285
     286
     287            }else{
     288        m = neighbour_vertices[ki];
     289                nm = n*2+m;
     290           
     291           
     292            qrp[0] = area_edge_values[nm];
     293            qrp[1] = discharge_edge_values[nm];
     294            qrp[2] = bed_edge_values[nm];
     295            qrp[3] = height_edge_values[nm];
     296            qrp[4] = velocity_edge_values[nm];
     297            qrp[5] = width_edge_values[nm];
     298            }
     299            qrm[0] = area_edge_values[ki];
     300            qrm[1] = discharge_edge_values[ki];
     301            qrm[2] = bed_edge_values[ki];
     302            qrm[3] = height_edge_values[ki];
     303            qrm[4] = velocity_edge_values[ki];
     304            qrm[5] = width_edge_values[ki];
     305
     306             _flux_function_channel(qlm,qlp,qrm,qrp,g,epsilon,h0,edgeflux,&max_speed);
    182307            flux[0] -= edgeflux[0];
    183308            flux[1] -= edgeflux[1];
    184            
     309          
    185310            // Update timestep based on edge i and possibly neighbour n
    186311            if (max_speed > epsilon) {
    187312                // Original CFL calculation
    188                
    189                 timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); 
     313
     314                timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
    190315                if (n>=0) {
    191                     timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); 
     316                    timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
    192317                }
    193318            }
    194         } // End edge i (and neighbour n)
     319        // End edge i (and neighbour n)
    195320        flux[0] /= areas[k];
    196         stage_explicit_update[k] = flux[0];
     321        area_explicit_update[k] = flux[0];
    197322        flux[1] /= areas[k];
    198         xmom_explicit_update[k] = flux[1];
    199        
     323        discharge_explicit_update[k] = flux[1];
    200324        //Keep track of maximal speeds
    201325        max_speed_array[k]=max_speed;
    202326    }
    203     return timestep;   
     327    return timestep;
     328
    204329}
     330
    205331
    206332//-------------------------------------------------------------
    … …  
    208334//------------------------------------------------------------
    209335//Innermost flux function (using w=z+h)
    210 int _flux_function(double *q_left, double *q_right,
    211         double z_left, double z_right,
    212                    double normals, double g, double epsilon, double h0,
    213                 double *edgeflux, double *max_speed) {
    214                
    215                 int i;
    216                 double ql[2], qr[2], flux_left[2], flux_right[2];
    217                 double z, w_left, h_left, uh_left, soundspeed_left, u_left;
    218                 double w_right, h_right, uh_right, soundspeed_right, u_right;
    219                 double s_max, s_min, denom;
    220                
    221                 //printf("h0 = %f \n",h0);
    222                 ql[0] = q_left[0];
    223                 ql[1] = q_left[1];
    224                 ql[1] = ql[1]*normals;
    225        
    226                 qr[0] = q_right[0];
    227                 qr[1] = q_right[1];
    228                 qr[1] = qr[1]*normals;
    229          
    230                 z = (z_left+z_right)/2.0;                 
    231                                    
    232                 //w_left = ql[0];
    233                 //h_left = w_left-z;
    234                 //uh_left = ql[1];
    235                
    236 
    237 
    238                 // Compute speeds in x-direction
    239                 w_left = ql[0];         
    240                 h_left = w_left-z;
    241                 uh_left = ql[1];
    242 
    243                 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
    244 
    245                 w_right = qr[0];
    246                 h_right = w_right-z;
    247                 uh_right = qr[1];
    248 
    249                 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
    250  
    251                 soundspeed_left = sqrt(g*h_left);
    252                 soundspeed_right = sqrt(g*h_right);
    253                
    254                 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
    255                 if (s_max < 0.0) s_max = 0.0;
    256        
    257                 s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    258                 if (s_min > 0.0) s_min = 0.0;
    259                
    260                
    261                 // Flux formulas
    262                 flux_left[0] = u_left*h_left;
    263                 flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
    264 
    265                 flux_right[0] = u_right*h_right;
    266                 flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
    267 
    268                 // Flux computation
    269                 denom = s_max-s_min;
    270                 if (denom < epsilon) {
    271                         for (i=0; i<2; i++) edgeflux[i] = 0.0;
    272                         *max_speed = 0.0;
    273                 } else {
    274                         edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
    275                         edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
    276                         edgeflux[0] /= denom;
    277                         edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
    278                         edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
    279                         edgeflux[1] /= denom;
    280                         edgeflux[1] *= normals;
    281    
    282                 // Maximal wavespeed
    283         *max_speed = max(fabs(s_max), fabs(s_min));
    284                 } 
    285     return 0;           
    286         }
    287                
    288                
    289        
    290        
     336
     337
     338
     339
     340
    291341// Computational function for flux computation
    292 double _compute_fluxes_ext(
    293                            double cfl,
    294                            double timestep,
    295                            double epsilon,
    296                            double g,
    297                            double h0,
    298                            long* neighbours,
    299                            long* neighbour_vertices,
    300                            double* normals,
    301                            double* areas,
    302                            double* stage_edge_values,
    303                            double* xmom_edge_values,
    304                            double* bed_edge_values,
    305                            double* stage_boundary_values,
    306                            double* xmom_boundary_values,
    307                            double* stage_explicit_update,
    308                            double* xmom_explicit_update,
    309                            int number_of_elements,
    310                            double* max_speed_array) {
    311                
    312                 double flux[2], ql[2], qr[2], edgeflux[2];
    313                 double zl, zr, max_speed, normal;
    314                 int k, i, ki, n, m, nm=0;
    315                
    316                
    317                 for (k=0; k<number_of_elements; k++) {
    318                         flux[0] = 0.0;
    319                         flux[1] = 0.0;
    320                        
    321                         for (i=0; i<2; i++) {
    322                                 ki = k*2+i;
    323                                
    324                                 ql[0] = stage_edge_values[ki];
    325                                 ql[1] = xmom_edge_values[ki];
    326                                 zl = bed_edge_values[ki];
    327                                
    328                                 n = neighbours[ki];
    329                                 if (n<0) {
    330                                         m = -n-1;
    331                                         qr[0] = stage_boundary_values[m];
    332                                         qr[1] = xmom_boundary_values[m];
    333                                         zr = zl;
    334                                 } else {
    335                                         m = neighbour_vertices[ki];
    336                                         nm = n*2+m;
    337                                         qr[0] = stage_edge_values[nm];
    338                                         qr[1] = xmom_edge_values[nm];
    339                                         zr = bed_edge_values[nm];                               
    340                                 }
    341                                
    342                                 normal = normals[ki];
    343                                 _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
    344                                 flux[0] -= edgeflux[0];
    345                                 flux[1] -= edgeflux[1];
    346                                
    347                                 // Update timestep based on edge i and possibly neighbour n
    348                                 if (max_speed > epsilon) {
    349                                     // Original CFL calculation
    350                                    
    351                                     timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
    352                                     if (n>=0) {
    353                                         timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
    354                                     }
    355                                 }
    356             } // End edge i (and neighbour n)
    357                         flux[0] /= areas[k];
    358                         stage_explicit_update[k] = flux[0];
    359                         flux[1] /= areas[k];
    360                         xmom_explicit_update[k] = flux[1];
    361                        
    362                         //Keep track of maximal speeds
    363                         max_speed_array[k]=max_speed;
    364                 }
    365                 return timestep;       
    366         }
     342
    367343
    368344//=========================================================================
    369345// Python Glue
    370346//=========================================================================
    371 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
    372  
    373    PyObject
     347
     348
     349
     350//------------------------------------------------
     351// New velocity based compute fluxes
     352//------------------------------------------------
     353
     354PyObject *compute_fluxes_channel_ext(PyObject *self, PyObject *args) {
     355
     356    PyObject
    374357        *domain,
    375         *stage,
    376         *xmom,
    377         *bed;
    378 
    379     PyArrayObject
    380         *neighbours,
     358        *area,
     359        *discharge,
     360        *bed,
     361        *height,
     362        *velocity,
     363        *width;
     364
     365    PyArrayObject
     366        *neighbours,
    381367        *neighbour_vertices,
    382         *normals, 
     368        *normals,
    383369        *areas,
    384         *stage_vertex_values,
    385         *xmom_vertex_values,
     370        *area_vertex_values,
     371        *discharge_vertex_values,
    386372        *bed_vertex_values,
    387         *stage_boundary_values,
    388         *xmom_boundary_values,
    389         *stage_explicit_update,
    390         *xmom_explicit_update,
     373        *height_vertex_values,
     374        *velocity_vertex_values,
     375        *width_vertex_values,
     376        *area_boundary_values,
     377        *discharge_boundary_values,
     378        *bed_boundary_values,
     379        *height_boundary_values,
     380        *velocity_boundary_values,
     381        *width_boundary_values,
     382        *area_explicit_update,
     383        *discharge_explicit_update,
    391384        *max_speed_array;
    392    
     385
    393386  double timestep, epsilon, g, h0, cfl;
    394387  int number_of_elements;
    395388
    396    
     389
    397390  // Convert Python arguments to C
    398   if (!PyArg_ParseTuple(args, "dOOOO",
     391  if (!PyArg_ParseTuple(args, "dOOOOOOO",
    399392                        &timestep,
    400393                        &domain,
    401                         &stage,
    402                         &xmom,
    403                         &bed)) {
    404       PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_ext could not parse input");
     394                        &area,
     395                        &discharge,
     396                        &bed,
     397                        &height,
     398                        &velocity,
     399                        &width)) {
     400      PyErr_SetString(PyExc_RuntimeError, "comp_flux_channel_ext.c: compute_fluxes_channel_ext could not parse input");
    405401      return NULL;
    406402  }
    … …  
    411407    h0                = get_python_double(domain,"h0");
    412408    cfl               = get_python_double(domain,"CFL");
    413  
    414    
     409
     410
    415411    neighbours        = get_consecutive_array(domain, "neighbours");
    416     neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
     412    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
    417413    normals           = get_consecutive_array(domain, "normals");
    418     areas             = get_consecutive_array(domain, "areas");   
     414    areas             = get_consecutive_array(domain, "areas");
    419415    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
    420    
    421     stage_vertex_values = get_consecutive_array(stage, "vertex_values");   
    422     xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");   
    423     bed_vertex_values   = get_consecutive_array(bed, "vertex_values");   
    424 
    425     stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
    426     xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
    427 
    428 
    429     stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
    430     xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
    431 
    432 
    433 
    434     number_of_elements = stage_vertex_values -> dimensions[0];
    435 
    436 
    437  
    438     // Call underlying flux computation routine and update
    439     // the explicit update arrays
    440     timestep = _compute_fluxes_ext(
    441         cfl,
    442         timestep,
    443         epsilon,
    444         g,
    445         h0,
    446         (long*) neighbours -> data,
    447         (long*) neighbour_vertices -> data,
    448         (double*) normals -> data,
    449         (double*) areas -> data,
    450         (double*) stage_vertex_values -> data,
    451         (double*) xmom_vertex_values -> data,
    452         (double*) bed_vertex_values -> data,
    453         (double*) stage_boundary_values -> data,
    454         (double*) xmom_boundary_values -> data,
    455         (double*) stage_explicit_update -> data,
    456         (double*) xmom_explicit_update -> data,
    457         number_of_elements,
    458         (double*) max_speed_array -> data);
    459 
    460 
    461   Py_DECREF(neighbours);
    462   Py_DECREF(neighbour_vertices);
    463   Py_DECREF(normals);
    464   Py_DECREF(areas);
    465   Py_DECREF(stage_vertex_values);
    466   Py_DECREF(xmom_vertex_values);
    467   Py_DECREF(bed_vertex_values);
    468   Py_DECREF(stage_boundary_values);
    469   Py_DECREF(xmom_boundary_values);
    470   Py_DECREF(stage_explicit_update);
    471   Py_DECREF(xmom_explicit_update);
    472   Py_DECREF(max_speed_array);
    473 
    474 
    475 
    476 
    477   // Return updated flux timestep
    478   return Py_BuildValue("d", timestep);
    479 }
    480 
    481 
    482 //------------------------------------------------
    483 // New velocity based compute fluxes
    484 //------------------------------------------------
    485 PyObject *compute_fluxes_vel_ext(PyObject *self, PyObject *args) {
    486  
    487     PyObject
    488         *domain,
    489         *stage,
    490         *xmom,
    491         *bed,
    492         *height,
    493         *velocity;
    494 
    495     PyArrayObject
    496         *neighbours,
    497         *neighbour_vertices,
    498         *normals,
    499         *areas,
    500         *stage_vertex_values,
    501         *xmom_vertex_values,
    502         *bed_vertex_values,
    503         *height_vertex_values,
    504         *velocity_vertex_values,
    505         *stage_boundary_values,
    506         *xmom_boundary_values,
    507         *bed_boundary_values,
    508         *height_boundary_values,
    509         *velocity_boundary_values,
    510         *stage_explicit_update,
    511         *xmom_explicit_update,
    512         *max_speed_array;
    513    
    514   double timestep, epsilon, g, h0, cfl;
    515   int number_of_elements;
    516 
    517    
    518   // Convert Python arguments to C
    519   if (!PyArg_ParseTuple(args, "dOOOOOO",
    520                         &timestep,
    521                         &domain,
    522                         &stage,
    523                         &xmom,
    524                         &bed,
    525                         &height,
    526                         &velocity)) {
    527       PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
    528       return NULL;
    529   }
    530 
    531 
    532     epsilon           = get_python_double(domain,"epsilon");
    533     g                 = get_python_double(domain,"g");
    534     h0                = get_python_double(domain,"h0");
    535     cfl               = get_python_double(domain,"CFL");
    536  
    537    
    538     neighbours        = get_consecutive_array(domain, "neighbours");
    539     neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
    540     normals           = get_consecutive_array(domain, "normals");
    541     areas             = get_consecutive_array(domain, "areas");   
    542     max_speed_array   = get_consecutive_array(domain, "max_speed_array");
    543    
    544     stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");   
    545     xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");   
    546     bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");   
    547     height_vertex_values     = get_consecutive_array(height,   "vertex_values");   
    548     velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");   
    549 
    550     stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");   
    551     xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");   
    552     bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");   
    553     height_boundary_values    = get_consecutive_array(height,    "boundary_values");   
    554     velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");   
    555 
    556 
    557     stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
    558     xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");   
    559 
    560     number_of_elements = stage_vertex_values -> dimensions[0];
    561  
    562     // Call underlying flux computation routine and update
    563     // the explicit update arrays
    564     timestep = _compute_fluxes_vel_ext(cfl,
     416
     417    area_vertex_values      = get_consecutive_array(area,    "vertex_values");
     418    discharge_vertex_values       = get_consecutive_array(discharge,     "vertex_values");
     419    bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
     420    height_vertex_values     = get_consecutive_array(height,   "vertex_values");
     421    velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
     422    width_vertex_values      = get_consecutive_array(width, "vertex_values");
     423
     424    area_boundary_values     = get_consecutive_array(area,     "boundary_values");
     425    discharge_boundary_values      = get_consecutive_array(discharge,      "boundary_values");
     426    bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
     427    height_boundary_values    = get_consecutive_array(height,    "boundary_values");
     428    velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
     429    width_boundary_values       = get_consecutive_array(width,     "boundary_values");
     430
     431
     432    area_explicit_update = get_consecutive_array(area, "explicit_update");
     433    discharge_explicit_update  = get_consecutive_array(discharge,  "explicit_update");
     434
     435    number_of_elements = area_vertex_values -> dimensions[0];
     436
     437    // Call underlying flux computation routine and update
     438    // the explicit update arrays
     439    timestep = _compute_fluxes_channel_ext(cfl,
    565440                                       timestep,
    566441                                       epsilon,
    … …  
    571446                                       (double*) normals -> data,
    572447                                       (double*) areas -> data,
    573                                        (double*) stage_vertex_values -> data,
    574                                        (double*) xmom_vertex_values -> data,
     448                                       (double*) area_vertex_values -> data,
     449                                       (double*) discharge_vertex_values -> data,
    575450                                       (double*) bed_vertex_values -> data,
    576451                                       (double*) height_vertex_values -> data,
    577452                                       (double*) velocity_vertex_values -> data,
    578                                        (double*) stage_boundary_values -> data,
    579                                        (double*) xmom_boundary_values -> data,
     453                                       (double*) width_vertex_values -> data,
     454                                       (double*) area_boundary_values -> data,
     455                                       (double*) discharge_boundary_values -> data,
    580456                                       (double*) bed_boundary_values -> data,
    581457                                       (double*) height_boundary_values -> data,
    582458                                       (double*) velocity_boundary_values -> data,
    583                                        (double*) stage_explicit_update -> data,
    584                                        (double*) xmom_explicit_update -> data,
     459                                       (double*) width_boundary_values -> data,
     460                                       (double*) area_explicit_update -> data,
     461                                       (double*) discharge_explicit_update -> data,
    585462                                       number_of_elements,
    586463                                       (double*) max_speed_array -> data);
    587    
    588    
     464
     465
    589466    Py_DECREF(neighbours);
    590467    Py_DECREF(neighbour_vertices);
    591468    Py_DECREF(normals);
    592469    Py_DECREF(areas);
    593     Py_DECREF(stage_vertex_values);
    594     Py_DECREF(xmom_vertex_values);
     470    Py_DECREF(area_vertex_values);
     471    Py_DECREF(discharge_vertex_values);
    595472    Py_DECREF(bed_vertex_values);
    596473    Py_DECREF(height_vertex_values);
    597474    Py_DECREF(velocity_vertex_values);
    598     Py_DECREF(stage_boundary_values);
    599     Py_DECREF(xmom_boundary_values);
     475    Py_DECREF(width_vertex_values);
     476    Py_DECREF(area_boundary_values);
     477    Py_DECREF(discharge_boundary_values);
    600478    Py_DECREF(bed_boundary_values);
    601479    Py_DECREF(height_boundary_values);
    602480    Py_DECREF(velocity_boundary_values);
    603     Py_DECREF(stage_explicit_update);
    604     Py_DECREF(xmom_explicit_update);
     481    Py_DECREF(width_boundary_values);
     482    Py_DECREF(area_explicit_update);
     483    Py_DECREF(discharge_explicit_update);
    605484    Py_DECREF(max_speed_array);
    606    
    607    
     485
     486
    608487    // Return updated flux timestep
    609488    return Py_BuildValue("d", timestep);
    … …  
    611490
    612491
    613 
    614492//-------------------------------
    615493// Method table for python module
    … …  
    617495
    618496static struct PyMethodDef MethodTable[] = {
    619   {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
    620   {"compute_fluxes_vel_ext", compute_fluxes_vel_ext, METH_VARARGS, "Print out"},
    621   {NULL, NULL}
     497  {"compute_fluxes_channel_ext", compute_fluxes_channel_ext, METH_VARARGS, "Print out"},
     498  {NULL}
    622499};
    623500
    624 // Module initialisation
    625 void initcomp_flux_vel_ext(void){
    626   Py_InitModule("comp_flux_vel_ext", MethodTable);
     501/* // Module initialisation */
     502/* void initcomp_flux_vel_ext(void){ */
     503/*   Py_InitModule("comp_flux_vel_ext", MethodTable); */
     504/*   import_array(); */
     505/* } */
     506
     507void initchannel_domain_ext(void){
     508  Py_InitModule("channel_domain_ext", MethodTable);
    627509  import_array();
    628510}
Note: See TracChangeset for help on using the changeset viewer.