Ignore:
Timestamp:
Apr 1, 2009, 8:41:37 PM (14 years ago)
Author:
sudi
Message:

Revised codes for discontinuous bed.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

    r5832 r6694  
    33#include "math.h"
    44#include <stdio.h>
     5#include "util_ext.h"
    56const double pi = 3.14159265358979;
    67
    7 
    8 // Shared code snippets
    9 #include "util_ext.h"
    10 
    11 
    12 //double max(double a, double b) {
    13 //      double z;
    14 //      z=(a>b)?a:b;
    15 //      return z;}
    16        
    17 //double min(double a, double b) {
    18 //      double z;
    19 //      z=(a<b)?a:b;
    20 //      return z;}
    21 
    22 
    23 
    24 
    25 
    268//Innermost flux function (using w=z+h)
    27 int _flux_function_wellbalanced(double *q_left, double *q_right,
    28                 double normals, double g, double epsilon,
    29                 double *edgeflux, double *max_speed) {
    30                
     9int _flux_function_wellbalanced(double *q_left,
     10                                                                double *q_right,
     11                                                                double normals,
     12                                                                double g,
     13                                                                double epsilon,
     14                                                                double *edgeflux,
     15                                                                double *max_speed) {           
    3116                int i;
    3217        double z_left, z_right;
    3318                double flux_left[2], flux_right[2];
    3419                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
    35                 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
     20                double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;
    3621                double s_max, s_min, denom;
    3722               
    38                 w_left = q_left[0];
     23                w_left  = q_left[0];
    3924                uh_left = q_left[1];
    4025                uh_left = uh_left*normals;
    41         z_left = q_left[2];
    42        
    43                 w_right = q_right[0];
     26        z_left  = q_left[2];
     27                h_left  = q_left[3];
     28                u_left  = q_left[4];
     29       
     30                w_right  = q_right[0];
    4431                uh_right = q_right[1];
    4532                uh_right = uh_right*normals;
    46         z_right = q_right[2];
     33        z_right  = q_right[2];
     34                h_right  = q_right[3];
     35                u_right  = q_right[4];
    4736         
    4837                if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
    4938       
    50                 //z = (z_left+z_right)/2.0;
    5139                z_star = max(z_left, z_right);                                                                                                          //equation(7) 
    5240                                   
    53                 // Compute speeds in x-direction
    54                 h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1].
    55                 h_left_star = max(0, w_left-z_star);                                                                                                            //(8)
    56    
     41                // Compute speeds in x-direction.
     42                h_left_star = max(0.0, w_left-z_star);                                                                                                          //equation(8)
     43               
    5744                if (h_left_star < epsilon) {
    5845                        u_left = 0.0; h_left_star = 0.0;
    … …  
    6148                        u_left = uh_left/h_left_star;
    6249                }
    63        
    64                 h_right = w_right-z_right;                                                                                                                                                                                                     
    65                 h_right_star = max(0, w_right-z_star);                                                                                                          //(9)
     50                                                                                                                                                                                                   
     51                h_right_star = max(0.0, w_right-z_star);                                                                                                        //equation(9)
    6652
    6753                if (h_right_star < epsilon) {
    … …  
    9278                // Flux computation
    9379                denom = s_max-s_min;
    94                 if (denom < epsilon) {
     80                if (denom < epsilon && z_right > z_left) {
    9581                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
     82                        // Maximal wavespeed                   
     83                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
     84                        edgeflux[1] *= normals;
    9685                        *max_speed = 0.0;
    97                 } else {
     86                } else if (denom < epsilon){
     87                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
     88                        // Maximal wavespeed
     89                        *max_speed = 0.0;               
     90                }else {
    9891                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
    9992                        edgeflux[0] += s_max*s_min*(h_right_star-h_left_star);                   //s_max*s_min*(qr[0]-ql[0]);
    … …  
    10295                        edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
    10396                        edgeflux[1] /= denom;
     97                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
    10498                        edgeflux[1] *= normals;
    105    
    106                 //edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star;
    107                 //edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star;
    108                 edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star;           //what about the h_right and h_right_star??????????????????????
    109                 // Maximal wavespeed
    110         *max_speed = max(fabs(s_max), fabs(s_min));
     99                        // Maximal wavespeed
     100                        *max_speed = max(fabs(s_max), fabs(s_min));
    111101                } 
    112                 //printf("max_speed = %g \n",*max_speed);
     102               
    113103    return 0;           
    114104        }
    … …  
    121111                double epsilon,
    122112                double g,
     113               
    123114                long* neighbours,
    124115                long* neighbour_vertices,
    125116                double* normals,
    126117                double* areas,
     118               
    127119                double* stage_edge_values,
    128120                double* xmom_edge_values,
    129121                double* bed_edge_values,
     122                double* height_edge_values,
     123                double* vel_edge_values,
     124               
    130125                double* stage_boundary_values,
    131126                double* xmom_boundary_values,
     127                double* bed_boundary_values,
     128                double* height_boundary_values,
     129                double* vel_boundary_values,
     130               
    132131                double* stage_explicit_update,
    133132                double* xmom_explicit_update,
     133                double* bed_explicit_update,
     134                double* height_explicit_update,
     135                double* vel_explicit_update,
     136               
    134137                int number_of_elements,
    135138                double* max_speed_array) {
    136139               
    137                 double flux[2], ql[3], qr[3], edgeflux[2];
     140                double flux[2], ql[5], qr[5], edgeflux[2];
    138141                double max_speed, normal;
    139142                int k, i, ki, n, m, nm=0;
    … …  
    149152                                ql[1] = xmom_edge_values[ki];
    150153                                ql[2] = bed_edge_values[ki];
    151                                 //ql[3] = velocity_edge_values[ki];
    152                 //ql[4] = height_edge_values[ki];
     154                ql[3] = height_edge_values[ki];
     155                                ql[4] = vel_edge_values[ki];
    153156               
    154157                                n = neighbours[ki];
     158                                printf("neighbours[ki] = %li \n",neighbours[ki]);
    155159                                if (n<0) {
    156160                                        m = -n-1;
    157161                                        qr[0] = stage_boundary_values[m];
    158162                                        qr[1] = xmom_boundary_values[m];
    159                                         qr[2] = ql[2];
     163                                        qr[2] = ql[2];                                                                                                                                                 
     164                                        qr[3] = height_edge_values[m];                                                                                                         
     165                                        qr[4] = vel_edge_values[m];                                                                                                                             
    160166                   
    161167                                } else {
    … …  
    164170                                        qr[0] = stage_edge_values[nm];
    165171                                        qr[1] = xmom_edge_values[nm];
    166                                         qr[2] = bed_edge_values[nm];                           
     172                                        qr[2] = bed_edge_values[nm];
     173                                        qr[3] = height_edge_values[nm];                                                                                                                 
     174                                        qr[4] = vel_edge_values[nm];                                                                                                                   
    167175                                }
    168176                               
    … …  
    175183                                if (max_speed > epsilon) {
    176184                        // Original CFL calculation
    177                                        
    178                                         timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
     185                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
    179186                                        if (n>=0) {
    180                                                 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
     187                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed.
    181188                                        }
    182189                                }
    … …  
    190197                        max_speed_array[k]=max_speed;
    191198                }
    192                 //printf("timestep = %f \n",timestep);
     199               
    193200                return timestep;       
    194201        }
    … …  
    209216        *normals,
    210217        *areas,
     218       
    211219        *stage_edge_values,
    212220        *xmom_edge_values,
    213221        *bed_edge_values,
     222        *height_edge_values,
     223        *vel_edge_values,
     224       
    214225        *stage_boundary_values,
    215226        *xmom_boundary_values,
     227        *bed_boundary_values,
     228        *height_boundary_values,
     229        *vel_boundary_values,
     230       
    216231        *stage_explicit_update,
    217232        *xmom_explicit_update,
     233        *bed_explicit_update,
     234        *height_explicit_update,
     235        *vel_explicit_update,
     236       
    218237        *max_speed_array;
    219238   
    … …  
    222241   
    223242  // Convert Python arguments to C
    224   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
     243  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
    225244                        &timestep,
    226245                        &epsilon,
    227246                        &g,
     247                       
    228248                        &neighbours,
    229249                        &neighbour_vertices,
    230250                        &normals,
    231251                        &areas,
     252                       
    232253                        &stage_edge_values,
    233254                        &xmom_edge_values,
    234255                        &bed_edge_values,
     256                        &height_edge_values,
     257                        &vel_edge_values,
     258                       
    235259                        &stage_boundary_values,
    236260                        &xmom_boundary_values,
     261                        &bed_boundary_values,
     262                        &height_boundary_values,
     263                        &vel_boundary_values,
     264                       
    237265                        &stage_explicit_update,
    238266                        &xmom_explicit_update,
     267                        &bed_explicit_update,
     268                        &height_explicit_update,
     269                        &vel_explicit_update,
     270                       
    239271                        &number_of_elements,
    240272                        &max_speed_array)) {
    241     PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
     273    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input");
    242274    return NULL;
    243275  }
    … …  
    249281                                 epsilon,
    250282                                 g,
     283                                 
    251284                                 (long*) neighbours -> data,
    252285                                 (long*) neighbour_vertices -> data,
    253286                                 (double*) normals -> data,
    254287                                 (double*) areas -> data,
     288                                 
    255289                                 (double*) stage_edge_values -> data,
    256290                                 (double*) xmom_edge_values -> data,
    257291                                 (double*) bed_edge_values -> data,
     292                                 (double*) height_edge_values -> data,
     293                                 (double*) vel_edge_values -> data,
     294                                 
    258295                                 (double*) stage_boundary_values -> data,
    259296                                 (double*) xmom_boundary_values -> data,
     297                                 (double*) bed_boundary_values -> data,
     298                                 (double*) height_boundary_values -> data,
     299                                 (double*) vel_boundary_values -> data,
     300                                 
    260301                                 (double*) stage_explicit_update -> data,
    261302                                 (double*) xmom_explicit_update -> data,
     303                                 (double*) bed_explicit_update -> data,
     304                                 (double*) height_explicit_update -> data,
     305                                 (double*) vel_explicit_update -> data,
     306                                 
    262307                                 number_of_elements,
    263308                                 (double*) max_speed_array -> data);
    … …  
    266311  return Py_BuildValue("d", timestep);
    267312}
    268 
    269 
    270313
    271314
    … …  
    285328  import_array();
    286329}
    287        
    288        
Note: See TracChangeset for help on using the changeset viewer.