Changeset 1594


Ignore:
Timestamp:
Jul 11, 2005, 12:15:59 PM (20 years ago)
Author:
chris
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r1562 r1594  
    4747}
    4848
    49 int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
    50   //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
    51   //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
    52   //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
    53   //and qc is the centroid
    54   //dq0=q(vertex0)-q(centroid of FV triangle)
    55   //dq1=q(vertex1)-q(vertex0)
    56   //dq2=q(vertex2)-q(vertex0)
    57   if (dq0>=0.0){
    58     if (dq1>=dq2){
    59       if (dq1>=0.0)
    60         *qmax=dq0+dq1;
    61       else
    62         *qmax=dq0;
    63       if ((*qmin=dq0+dq2)<0)
    64         ;//qmin is already set to correct value
    65       else
    66         *qmin=0.0;
    67     }
    68     else{//dq1<dq2
    69       if (dq2>0)
    70         *qmax=dq0+dq2;
    71       else
    72         *qmax=dq0;
    73       if ((*qmin=dq0+dq1)<0)
    74         ;//qmin is the correct value
    75       else
    76         *qmin=0.0;
    77     }
    78   }
    79   else{//dq0<0
    80     if (dq1<=dq2){
    81       if (dq1<0.0)
    82         *qmin=dq0+dq1;
    83       else
    84         *qmin=dq0;
    85       if ((*qmax=dq0+dq2)>0.0)
    86         ;//qmax is already set to the correct value
    87       else
    88         *qmax=0.0;
    89     }
    90     else{//dq1>dq2
    91       if (dq2<0.0)
    92         *qmin=dq0+dq2;
    93       else
    94         *qmin=dq0;
    95       if ((*qmax=dq0+dq1)>0.0)
    96         ;//qmax is already set to the correct value
    97       else
    98         *qmax=0.0;
    99     }
    100   }
    101   return 0;
    102 }
    103 
    104 int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
    105   //given provisional jumps dqv from the FV triangle centroid to its vertices and
    106   //jumps qmin (qmax) between the centroid of the FV triangle and the
    107   //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
    108   //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
    109   int i;
    110   double r=1000.0, r0=1.0, phi=1.0;
    111   static double TINY = 1.0e-100;//to avoid machine accuracy problems.
    112   //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
    113   for (i=0;i<3;i++){
    114     if (dqv[i]<-TINY)
    115       r0=qmin/dqv[i];
    116     if (dqv[i]>TINY)
    117       r0=qmax/dqv[i];
    118     r=min(r0,r);
    119     //
    120   }
    121   phi=min(r*beta_w,1.0);
    122   for (i=0;i<3;i++)
    123     dqv[i]=dqv[i]*phi;
    124   return 0;
    125 }
     49
    12650
    12751// Computational function for flux computation (using stage w=z+h)
    … …  
    246170      if (h >= eps) {
    247171        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    248         //S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    249         S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     172        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    250173        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    251174
    … …  
    503426}
    504427
    505 PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
    506   /*Compute the vertex values based on a linear reconstruction on each triangle
    507     These values are calculated as follows:
    508     1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
    509     formed by the centroids of its three neighbours.
    510     2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
    511     of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
    512     of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
    513     original triangle has gradient (a,b).
    514     3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
    515     find_qmin_and_qmax and limit_gradient
    516 
    517     Python call:
    518     extrapolate_second_order_sw(domain.surrogate_neighbours,
    519                                 domain.number_of_boundaries
    520                                 domain.centroid_coordinates,
    521                                 Stage.centroid_values
    522                                 Xmom.centroid_values
    523                                 Ymom.centroid_values
    524                                 domain.vertex_coordinates,
    525                                 Stage.vertex_values,
    526                                 Xmom.vertex_values,
    527                                 Ymom.vertex_values)
    528 
    529     Post conditions:
    530             The vertices of each triangle have values from a limited linear reconstruction
    531             based on centroid values
    532 
    533   */
    534   PyArrayObject *surrogate_neighbours,
    535     *number_of_boundaries,
    536     *centroid_coordinates,
    537     *stage_centroid_values,
    538     *xmom_centroid_values,
    539     *ymom_centroid_values,
    540     *vertex_coordinates,
    541     *stage_vertex_values,
    542     *xmom_vertex_values,
    543     *ymom_vertex_values;
    544   PyObject *domain, *Tmp;
    545   //Local variables
    546   double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
    547   int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
    548   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
    549   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    550   double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
    551   //by which these jumps are limited
    552   // Convert Python arguments to C
    553   if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
    554                         &domain,
    555                         &surrogate_neighbours,
    556                         &number_of_boundaries,
    557                         &centroid_coordinates,
    558                         &stage_centroid_values,
    559                         &xmom_centroid_values,
    560                         &ymom_centroid_values,
    561                         &vertex_coordinates,
    562                         &stage_vertex_values,
    563                         &xmom_vertex_values,
    564                         &ymom_vertex_values)) {
    565     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    566     return NULL;
    567   }
    568 
    569   //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
    570   Tmp = PyObject_GetAttrString(domain, "beta_w");
    571   if (!Tmp)
    572     return NULL;
    573   beta_w = PyFloat_AsDouble(Tmp);
    574   Py_DECREF(Tmp);
    575   number_of_elements = stage_centroid_values -> dimensions[0];
    576   for (k=0; k<number_of_elements; k++) {
    577     k3=k*3;
    578     k6=k*6;
    579 
    580     if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
    581       ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
    582       ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
    583       ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
    584       ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
    585       ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
    586       ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
    587       ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
    588       ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
    589       ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
    590       continue;
    591     }
    592     else {//we will need centroid coordinates and vertex coordinates of the triangle
    593       //get the vertex coordinates of the FV triangle
    594       xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
    595       xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
    596       xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
    597       //get the centroid coordinates of the FV triangle
    598       coord_index=2*k;
    599       x=((double *)centroid_coordinates->data)[coord_index];
    600       y=((double *)centroid_coordinates->data)[coord_index+1];
    601       //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
    602       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    603       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
    604     }
    605     if (((long *)number_of_boundaries->data)[k]<=1){
    606       //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
    607       //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
    608       k0=((long *)surrogate_neighbours->data)[k3];
    609       k1=((long *)surrogate_neighbours->data)[k3+1];
    610       k2=((long *)surrogate_neighbours->data)[k3+2];
    611       //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
    612       coord_index=2*k0;
    613       x0=((double *)centroid_coordinates->data)[coord_index];
    614       y0=((double *)centroid_coordinates->data)[coord_index+1];
    615       coord_index=2*k1;
    616       x1=((double *)centroid_coordinates->data)[coord_index];
    617       y1=((double *)centroid_coordinates->data)[coord_index+1];
    618       coord_index=2*k2;
    619       x2=((double *)centroid_coordinates->data)[coord_index];
    620       y2=((double *)centroid_coordinates->data)[coord_index+1];
    621       //store x- and y- differentials for the vertices of the auxiliary triangle
    622       dx1=x1-x0; dx2=x2-x0;
    623       dy1=y1-y0; dy2=y2-y0;
    624       //calculate 2*area of the auxiliary triangle
    625       area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
    626       //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
    627       if (area2<=0)
    628         return NULL;
    629 
    630       //### stage ###
    631       //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    632       dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
    633       //calculate differentials between the vertices of the auxiliary triangle
    634       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
    635       dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
    636       //calculate the gradient of stage on the auxiliary triangle
    637       a = dy2*dq1 - dy1*dq2;
    638       a /= area2;
    639       b = dx1*dq2 - dx2*dq1;
    640       b /= area2;
    641       //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    642       dqv[0]=a*dxv0+b*dyv0;
    643       dqv[1]=a*dxv1+b*dyv1;
    644       dqv[2]=a*dxv2+b*dyv2;
    645       //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    646       //and compute jumps from the centroid to the min and max
    647       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    648       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    649       for (i=0;i<3;i++)
    650         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    651 
    652       //### xmom ###
    653       //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    654       dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
    655       //calculate differentials between the vertices of the auxiliary triangle
    656       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
    657       dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
    658       //calculate the gradient of xmom on the auxiliary triangle
    659       a = dy2*dq1 - dy1*dq2;
    660       a /= area2;
    661       b = dx1*dq2 - dx2*dq1;
    662       b /= area2;
    663       //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    664       dqv[0]=a*dxv0+b*dyv0;
    665       dqv[1]=a*dxv1+b*dyv1;
    666       dqv[2]=a*dxv2+b*dyv2;
    667       //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    668       //and compute jumps from the centroid to the min and max
    669       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    670       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    671       for (i=0;i<3;i++)
    672         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    673 
    674       //### ymom ###
    675       //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
    676       dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
    677       //calculate differentials between the vertices of the auxiliary triangle
    678       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
    679       dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
    680       //calculate the gradient of xmom on the auxiliary triangle
    681       a = dy2*dq1 - dy1*dq2;
    682       a /= area2;
    683       b = dx1*dq2 - dx2*dq1;
    684       b /= area2;
    685       //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
    686       dqv[0]=a*dxv0+b*dyv0;
    687       dqv[1]=a*dxv1+b*dyv1;
    688       dqv[2]=a*dxv2+b*dyv2;
    689       //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
    690       //and compute jumps from the centroid to the min and max
    691       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
    692       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    693       for (i=0;i<3;i++)
    694         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    695     }//if (number_of_boundaries[k]<=1)
    696     else{//number_of_boundaries==2
    697       //one internal neighbour and gradient is in direction of the neighbour's centroid
    698       //find the only internal neighbour
    699       for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
    700         if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
    701           break;
    702       }
    703       if ((k2==k3+3))//if we didn't find an internal neighbour
    704         return NULL;//error
    705       k1=((long *)surrogate_neighbours->data)[k2];
    706       //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
    707       coord_index=2*k1;
    708       x1=((double *)centroid_coordinates->data)[coord_index];
    709       y1=((double *)centroid_coordinates->data)[coord_index+1];
    710       //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
    711       dx1=x1-x; dy1=y1-y;
    712       //set area2 as the square of the distance
    713       area2=dx1*dx1+dy1*dy1;
    714       //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
    715       //respectively correspond to the x- and y- gradients of the conserved quantities
    716       dx2=1.0/area2;
    717       dy2=dx2*dy1;
    718       dx2*=dx1;
    719 
    720       //## stage ###
    721       //compute differentials
    722       dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
    723       //calculate the gradient between the centroid of the FV triangle and that of its neighbour
    724       a=dq1*dx2;
    725       b=dq1*dy2;
    726       //calculate provisional vertex jumps, to be limited
    727       dqv[0]=a*dxv0+b*dyv0;
    728       dqv[1]=a*dxv1+b*dyv1;
    729       dqv[2]=a*dxv2+b*dyv2;
    730       //now limit the jumps
    731       if (dq1>=0.0){
    732         qmin=0.0;
    733         qmax=dq1;
    734       }
    735       else{
    736         qmin=dq1;
    737         qmax=0.0;
    738       }
    739       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    740       for (i=0;i<3;i++)
    741         ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
    742 
    743       //## xmom ###
    744       //compute differentials
    745       dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
    746       //calculate the gradient between the centroid of the FV triangle and that of its neighbour
    747       a=dq1*dx2;
    748       b=dq1*dy2;
    749       //calculate provisional vertex jumps, to be limited
    750       dqv[0]=a*dxv0+b*dyv0;
    751       dqv[1]=a*dxv1+b*dyv1;
    752       dqv[2]=a*dxv2+b*dyv2;
    753       //now limit the jumps
    754       if (dq1>=0.0){
    755         qmin=0.0;
    756         qmax=dq1;
    757       }
    758       else{
    759         qmin=dq1;
    760         qmax=0.0;
    761       }
    762       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    763       for (i=0;i<3;i++)
    764         ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
    765 
    766       //## ymom ###
    767       //compute differentials
    768       dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
    769       //calculate the gradient between the centroid of the FV triangle and that of its neighbour
    770       a=dq1*dx2;
    771       b=dq1*dy2;
    772       //calculate provisional vertex jumps, to be limited
    773       dqv[0]=a*dxv0+b*dyv0;
    774       dqv[1]=a*dxv1+b*dyv1;
    775       dqv[2]=a*dxv2+b*dyv2;
    776       //now limit the jumps
    777       if (dq1>=0.0){
    778         qmin=0.0;
    779         qmax=dq1;
    780       }
    781       else{
    782         qmin=dq1;
    783         qmax=0.0;
    784       }
    785       limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
    786       for (i=0;i<3;i++)
    787         ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
    788     }//else [number_of_boudaries==2]
    789   }//for k=0 to number_of_elements-1
    790   return Py_BuildValue("");
    791 }//extrapolate_second-order_sw
    792 
    793428PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
    794429  //
    … …  
    848483}
    849484
     485
    850486PyObject *compute_fluxes(PyObject *self, PyObject *args) {
    851487  /*Compute all fluxes and the timestep suitable for all volumes
    … …  
    866502    domain.timestep = compute_fluxes(timestep,
    867503                                     domain.epsilon,
    868                                      domain.g,
     504                                     domain.g,
    869505                                     domain.neighbours,
    870506                                     domain.neighbour_edges,
    … …  
    882518                                     Stage.explicit_update,
    883519                                     Xmom.explicit_update,
    884                                      Ymom.explicit_update,
    885                                      already_computed_flux)
     520                                     Ymom.explicit_update)
    886521
    887522
    … …  
    905540    *stage_explicit_update,
    906541    *xmom_explicit_update,
    907     *ymom_explicit_update,
    908     *already_computed_flux;//tracks whether the flux across an edge has already been computed
     542    *ymom_explicit_update;
    909543
    910544
    … …  
    912546  double timestep, max_speed, epsilon, g;
    913547  double normal[2], ql[3], qr[3], zl, zr;
    914   double edgeflux[3]; //Work arrays for summing up fluxes
    915 
    916   int number_of_elements, k, i, m, n;
    917   int ki, nm=0, ki2; //Index shorthands
    918   static long call=1;
     548  double flux[3], edgeflux[3]; //Work arrays for summing up fluxes
     549
     550  int number_of_elements, k, i, j, m, n;
     551  int ki, nm, ki2; //Index shorthands
    919552
    920553
    921554  // Convert Python arguments to C
    922   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
     555  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
    923556                        &timestep,
    924557                        &epsilon,
    … …  
    937570                        &stage_explicit_update,
    938571                        &xmom_explicit_update,
    939                         &ymom_explicit_update,
    940                         &already_computed_flux)) {
     572                        &ymom_explicit_update)) {
    941573    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    942574    return NULL;
    943575  }
     576
    944577  number_of_elements = stage_edge_values -> dimensions[0];
    945   call++;//a static local variable to which already_computed_flux is compared
    946   //set explicit_update to zero for all conserved_quantities.
    947   //This assumes compute_fluxes called before forcing terms
     578
     579
    948580  for (k=0; k<number_of_elements; k++) {
    949     ((double *) stage_explicit_update -> data)[k]=0.0;
    950     ((double *) xmom_explicit_update -> data)[k]=0.0;
    951     ((double *) ymom_explicit_update -> data)[k]=0.0;
    952   }
    953   //Loop through neighbours and compute edge flux for each
    954   for (k=0; k<number_of_elements; k++) {
     581
     582    //Reset work array
     583    for (j=0; j<3; j++) flux[j] = 0.0;
     584
     585    //Loop through neighbours and compute edge flux for each
    955586    for (i=0; i<3; i++) {
    956587      ki = k*3+i;
    957       if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
    958         continue;
    959588      ql[0] = ((double *) stage_edge_values -> data)[ki];
    960589      ql[1] = ((double *) xmom_edge_values -> data)[ki];
    … …  
    972601      } else {
    973602        m = ((long *) neighbour_edges -> data)[ki];
     603
    974604        nm = n*3+m;
    975605        qr[0] = ((double *) stage_edge_values -> data)[nm];
    … …  
    978608        zr =    ((double *) bed_edge_values -> data)[nm];
    979609      }
     610
     611      //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
     612      //             ql[0], ql[1], ql[2]);
     613
     614
    980615      // Outward pointing normal vector
    981616      // normal = domain.normals[k, 2*i:2*i+2]
    … …  
    983618      normal[0] = ((double *) normals -> data)[ki2];
    984619      normal[1] = ((double *) normals -> data)[ki2+1];
     620
    985621      //Edge flux computation
    986622      flux_function(ql, qr, zl, zr,
    … …  
    988624                    epsilon, g,
    989625                    edgeflux, &max_speed);
    990       //update triangle k
    991       ((long *) already_computed_flux->data)[ki]=call;
    992       ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
    993       ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
    994       ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
    995       //update the neighbour n
    996       if (n>=0){
    997         ((long *) already_computed_flux->data)[nm]=call;
    998         ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
    999         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
    1000         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
     626
     627
     628      //flux -= edgeflux * edgelengths[k,i]
     629      for (j=0; j<3; j++) {
     630        flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
    1001631      }
    1002       ///for (j=0; j<3; j++) {
    1003         ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
    1004         ///}
    1005         //Update timestep
    1006         //timestep = min(timestep, domain.radii[k]/max_speed)
    1007         //FIXME: SR Add parameter for CFL condition
    1008         if (max_speed > epsilon) {
    1009           timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    1010           //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
    1011           if (n>=0)
    1012             timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    1013         }
     632
     633      //Update timestep
     634      //timestep = min(timestep, domain.radii[k]/max_speed)
     635      //FIXME: SR Add parameter for CFL condition
     636      if (max_speed > epsilon) {
     637        timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     638      }
    1014639    } // end for i
     640
    1015641    //Normalise by area and store for when all conserved
    1016642    //quantities get updated
    1017     ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
    1018     ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
    1019     ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     643    // flux /= areas[k]
     644    for (j=0; j<3; j++) {
     645      flux[j] /= ((double *) areas -> data)[k];
     646    }
     647
     648    ((double *) stage_explicit_update -> data)[k] = flux[0];
     649    ((double *) xmom_explicit_update -> data)[k] = flux[1];
     650    ((double *) ymom_explicit_update -> data)[k] = flux[2];
     651
    1020652  } //end for k
     653
    1021654  return Py_BuildValue("d", timestep);
    1022655}
     656
     657
    1023658
    1024659PyObject *protect(PyObject *self, PyObject *args) {
    … …  
    1178813}
    1179814
    1180 PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
    1181   //a faster version of h_limiter above
    1182   PyObject *domain, *Tmp;
    1183   PyArrayObject
    1184     *hv, *hc, //Depth at vertices and centroids
    1185     *hvbar,   //Limited depth at vertices (return values)
    1186     *neighbours;
    1187 
    1188   int k, i, N, k3,k0,k1,k2;
    1189   int dimensions[2];
    1190   double beta_h; //Safety factor (see config.py)
    1191   double hmin, hmax, dh[3];
    1192   // Convert Python arguments to C
    1193   if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
    1194     return NULL;
    1195   neighbours = get_consecutive_array(domain, "neighbours");
    1196 
    1197   //Get safety factor beta_h
    1198   Tmp = PyObject_GetAttrString(domain, "beta_h");
    1199   if (!Tmp)
    1200     return NULL;
    1201   beta_h = PyFloat_AsDouble(Tmp);
    1202 
    1203   Py_DECREF(Tmp);
    1204 
    1205   N = hc -> dimensions[0];
    1206 
    1207   //Create hvbar
    1208   dimensions[0] = N;
    1209   dimensions[1] = 3;
    1210   hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    1211   for (k=0;k<N;k++){
    1212     k3=k*3;
    1213     //get the ids of the neighbours
    1214     k0 = ((long*) neighbours -> data)[k3];
    1215     k1 = ((long*) neighbours -> data)[k3+1];
    1216     k2 = ((long*) neighbours -> data)[k3+2];
    1217     //set hvbar provisionally
    1218     for (i=0;i<3;i++){
    1219       ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
    1220       dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
    1221     }
    1222     hmin=((double*) hc -> data)[k];
    1223     hmax=hmin;
    1224     if (k0>=0){
    1225       hmin=min(hmin,((double*) hc -> data)[k0]);
    1226       hmax=max(hmax,((double*) hc -> data)[k0]);
    1227     }
    1228     if (k1>=0){
    1229       hmin=min(hmin,((double*) hc -> data)[k1]);
    1230       hmax=max(hmax,((double*) hc -> data)[k1]);
    1231     }
    1232     if (k2>=0){
    1233       hmin=min(hmin,((double*) hc -> data)[k2]);
    1234       hmax=max(hmax,((double*) hc -> data)[k2]);
    1235     }
    1236     hmin-=((double*) hc -> data)[k];
    1237     hmax-=((double*) hc -> data)[k];
    1238     limit_gradient(dh,hmin,hmax,beta_h);
    1239     for (i=0;i<3;i++)
    1240       ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
    1241   }
    1242   return PyArray_Return(hvbar);
    1243 }
     815
     816
    1244817
    1245818PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
    … …  
    1292865
    1293866  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
    1294   {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
    1295867  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
    1296868  {"gravity", gravity, METH_VARARGS, "Print out"},
    … …  
    1299871   METH_VARARGS, "Print out"},
    1300872  {"h_limiter", h_limiter,
    1301    METH_VARARGS, "Print out"},
    1302   {"h_limiter_sw", h_limiter_sw,
    1303873   METH_VARARGS, "Print out"},
    1304874  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
Note: See TracChangeset for help on using the changeset viewer.