Changeset 1556 for inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
- Timestamp:
- Jun 29, 2005, 4:49:54 PM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
r1509 r1556 48 48 49 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 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 52 //four values (at the centroid of the FV triangle and the auxiliary triangle vertices), 53 53 //and qc is the centroid … … 71 71 else 72 72 *qmax=dq0; 73 if ((*qmin=dq0+dq1)<0) 73 if ((*qmin=dq0+dq1)<0) 74 74 ;//qmin is the correct value 75 75 else … … 93 93 else 94 94 *qmin=dq0; 95 if ((*qmax=dq0+dq1)>0.0) 95 if ((*qmax=dq0+dq1)>0.0) 96 96 ;//qmax is already set to the correct value 97 else 97 else 98 98 *qmax=0.0; 99 99 } … … 109 109 int i; 110 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. 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 113 for (i=0;i<3;i++){ 114 114 if (dqv[i]<-TINY) … … 508 508 1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle 509 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 510 2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product 511 511 of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average 512 512 of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the … … 525 525 Stage.vertex_values, 526 526 Xmom.vertex_values, 527 Ymom.vertex_values) 527 Ymom.vertex_values) 528 528 529 529 Post conditions: 530 The vertices of each triangle have values from a limited linear reconstruction 530 The vertices of each triangle have values from a limited linear reconstruction 531 531 based on centroid values 532 532 533 533 */ 534 PyArrayObject *surrogate_neighbours, 534 PyArrayObject *surrogate_neighbours, 535 535 *number_of_boundaries, 536 536 *centroid_coordinates, … … 592 592 else {//we will need centroid coordinates and vertex coordinates of the triangle 593 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]; 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 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]; 598 coord_index=2*k; 599 x=((double *)centroid_coordinates->data)[coord_index]; 600 y=((double *)centroid_coordinates->data)[coord_index+1]; 601 601 //store x- and y- differentials for the vertices of the FV triangle relative to the centroid 602 602 dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x; … … 604 604 } 605 605 if (((long *)number_of_boundaries->data)[k]<=1){ 606 //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours 606 //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours 607 607 //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours 608 608 k0=((long *)surrogate_neighbours->data)[k3]; … … 668 668 //and compute jumps from the centroid to the min and max 669 669 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 670 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 670 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 671 671 for (i=0;i<3;i++) 672 672 ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i]; … … 690 690 //and compute jumps from the centroid to the min and max 691 691 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 692 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 692 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 693 693 for (i=0;i<3;i++) 694 694 ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i]; … … 706 706 //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1) 707 707 coord_index=2*k1; 708 x1=((double *)centroid_coordinates->data)[coord_index]; 708 x1=((double *)centroid_coordinates->data)[coord_index]; 709 709 y1=((double *)centroid_coordinates->data)[coord_index+1]; 710 710 //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour … … 717 717 dy2=dx2*dy1; 718 718 dx2*=dx1; 719 719 720 720 //## stage ### 721 721 //compute differentials … … 737 737 qmax=0.0; 738 738 } 739 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 739 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 740 740 for (i=0;i<3;i++) 741 741 ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i]; … … 760 760 qmax=0.0; 761 761 } 762 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 762 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 763 763 for (i=0;i<3;i++) 764 764 ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i]; … … 783 783 qmax=0.0; 784 784 } 785 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 785 limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited 786 786 for (i=0;i<3;i++) 787 787 ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i]; … … 866 866 domain.timestep = compute_fluxes(timestep, 867 867 domain.epsilon, 868 868 domain.g, 869 869 domain.neighbours, 870 870 domain.neighbour_edges, … … 883 883 Xmom.explicit_update, 884 884 Ymom.explicit_update, 885 885 already_computed_flux) 886 886 887 887 … … 944 944 number_of_elements = stage_edge_values -> dimensions[0]; 945 945 call++;//a static local variable to which already_computed_flux is compared 946 //set explicit_update to zero for all conserved_quantities. 946 //set explicit_update to zero for all conserved_quantities. 947 947 //This assumes compute_fluxes called before forcing terms 948 948 for (k=0; k<number_of_elements; k++) { … … 961 961 ql[2] = ((double *) ymom_edge_values -> data)[ki]; 962 962 zl = ((double *) bed_edge_values -> data)[ki]; 963 963 964 964 //Quantities at neighbour on nearest face 965 965 n = ((long *) neighbours -> data)[ki]; … … 1017 1017 ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1018 1018 ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1019 ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1019 ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k]; 1020 1020 } //end for k 1021 1021 return Py_BuildValue("d", timestep); … … 1194 1194 return NULL; 1195 1195 neighbours = get_consecutive_array(domain, "neighbours"); 1196 1196 1197 1197 //Get safety factor beta_h 1198 1198 Tmp = PyObject_GetAttrString(domain, "beta_h");
Note: See TracChangeset
for help on using the changeset viewer.
