TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
b_volumes.c
Go to the documentation of this file.
1 
27 #include "turtle.h"
28 #include "t_alloc.h"
29 #include "t_utilities.h"
30 #include "rw_maps.h"
31 #include "geometry.h"
32 #include "g_raster2plvector.h"
33 #include "bigcells2.h"
34 #include "geometry2.h"
35 #include "b_utilities.h"
36 #include "b_solver.h"
37 #include "keywords_file_b.h"
38 #include "b_volumes.h"
39 
40 extern STRINGBIN *filenames;
41 
42 
43 extern DOUBLERASTER_MAP *draster;
44 // extern DOUBLESQUARE_GRID *dsq;
45 extern DOUBLE_GRID *dgrid; // ec 20100413
46 extern PARAM *param;
47 extern FLAG *flag;
48 
49 extern DOUBLEVECTOR *elevation_bottom_fine;
50 extern DOUBLEVECTOR *elevation_bottom_coarse;
51 extern DOUBLEVECTOR *elevation_bottom_flines;
52 extern DOUBLEVECTOR *elevation_bottom_bottom;
53 extern DOUBLEVECTOR *elevation_terrain_fine;
54 extern DOUBLEVECTOR *elevation_terrain_flines;
55 extern DOUBLEVECTOR *outlet_coefficient;
56 extern DOUBLEVECTOR *outlet_coefficient_surf;
57 extern DOUBLEVECTOR *porosity_fine;
59 extern DOUBLEVECTOR *water_surface_elevation;
60 extern DOUBLEVECTOR *water_mass_error;
62 //double Ks=10;
63 double volume (double eta, long i) {
73  double vol=volume_subs(eta,i)+volume_surf(eta,i);
74 
75  return vol;
76 
77 }
78 
79 
80 double volume_surf(double eta, long i) {
96  double vol=0.0;
97  long r,c,kp;
98 // long novalue=-99.0; //dgrid->fine->novalue;
99  double hval,aval;
100  char *function_name="volume_surf";
101 
102  if ((i<dgrid->small_polygon_content->index->nl) || (i>dgrid->small_polygon_content->index->nh)) printf("Error in volume function: volume %ld does not exist [%ld,%ld]. \n ",i,dgrid->small_polygon_content->index->nl,dgrid->small_polygon_content->index->nh);
103 
104  for(r=NL;r<=dgrid->small_polygon_content->index->element[i];r++) {
105  // for(r=dsq->small_content_polygon->element[i]->nrl;r<=dsq->small_content_polygon->element[i]->nrh;r++){
106  // for(c=dsq->small_content_polygon->element[i]->ncl;c<=dsq->small_content_polygon->element[i]->nch;c++){
107  // kp=dsq->small_content_polygon->element[i]->element[r][c];
108  kp=dgrid->small_polygon_content->element[i][r];
109  //if (kp!=novalue){
110 
111  if ((kp<elevation_bottom_fine->nl) || (kp>elevation_bottom_fine->nh)) {
112  aval=0.0;
113  hval=0.0;
114  // printf ("Error in volume 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! novalue=%ld \n",kp,elevation_bottom_fine->nh,i,r,c,novalue);
115  printf ("Error in function %s volume 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! \n",function_name,kp,elevation_bottom_fine->nh,i,r,c);
116  stop_execution();
117  }else {
118  // aval=dsq->fine->grid->polygons->element[kp]->area2D;
119  aval=dgrid->fine->polygons->element[kp]->area2D;
120  hval=(eta-elevation_terrain_fine->element[kp])*aval;
121  }
122 
123  if (hval>0.0){
124  vol+=hval;
125 
126  }
127  // }
128  // }
129  }
130 
131 
132 
133  return vol;
134 }
135 
136 double volume_subs(double eta, long i) {
152  double vol=0.0;
153  long r,c,kp;
154 // long novalue=dsq->fine->novalue;
155  double hval,aval;
156  char *function_name="volume_subs";
157  if ((i<dgrid->small_polygon_content->index->nl) || (i>dgrid->small_polygon_content->index->nh)) printf("Error in volume function: volume %ld does not exist [%ld,%ld]. \n ",i,dgrid->small_polygon_content->index->nl,dgrid->small_polygon_content->index->nh);
158 
159  for(r=NL;r<=dgrid->small_polygon_content->index->element[i];r++){
160 // for(r=dsq->small_content_polygon->element[i]->nrl;r<=dsq->small_content_polygon->element[i]->nrh;r++){
161  // for(c=dsq->small_content_polygon->element[i]->ncl;c<=dsq->small_content_polygon->element[i]->nch;c++){
162  // kp=dsq->small_content_polygon->element[i]->element[r][c];
163  kp=dgrid->small_polygon_content->element[i][r];
164  // if (kp!=novalue){
165 
166  if ((kp<elevation_bottom_fine->nl) || (kp>elevation_bottom_fine->nh)) {
167  aval=0.0;
168  hval=0.0;
169  // printf ("Error in volume 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! novalue=%ld \n",kp,elevation_bottom_fine->nh,i,r,c,novalue);
170  printf ("Error in %s volume 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! \n",function_name,kp,elevation_bottom_fine->nh,i,r,c);
171  stop_execution();
172  }else {
173  // aval=dsq->fine->grid->polygons->element[kp]->area2D*porosity_fine->element[kp];
174  aval=dgrid->fine->polygons->element[kp]->area2D*porosity_fine->element[kp];
175  hval=(fmin(eta,elevation_terrain_fine->element[kp])-elevation_bottom_fine->element[kp])*aval;
176  }
177 
178  if (hval>0.0){
179  vol+=hval;
180 
181  }
182  // }
183  //}
184  }
185 
186 
187 
188  return vol;
189 }
190 
191 //double vertical_area()
192 
193 double vertical_area(double eta, long j) {
204  double area_vert=vertical_area_subs(eta,j);
205 
206  return area_vert;
207 
208 }
209 
210 
211 double vertical_area_subs(double eta, long j) {
224  long c,kp;
225  double area_vert=0.0;
226  long novalue=dgrid->novalue;
227  long nl_l=dgrid->fine->lines->nl; // dsq->fine->grid->lines->nl;
228  long nh_l=dgrid->fine->lines->nh; //dsq->fine->grid->lines->nh;
229 
230  if ((j>dgrid->small_line_content->index->nh) || (j<dgrid->small_line_content->index->nl)) printf("Error in vertical_area_subs function: Line %ld does not exist [%ld,%ld] . \n",j,dgrid->small_line_content->index->nl,dgrid->small_line_content->index->nh);
231 
232  for(c=NL;c<=dgrid->small_line_content->index->element[j];c++) {
233 
234  kp=dgrid->small_line_content->element[j][c];
235 
236  if (((kp<nl_l) || (kp>nh_l)) && (kp!=novalue) ) { // && (kp!=novalue)
237  printf ("Error in vertical area lines %ld (of %ld) cannot belong to the big line %ld (c=%ld)!! \n",kp,nh_l,j,c);
238  stop_execution();
239  } else if (kp!=novalue) {
240 // } else {
241  // area_vert+=fmax(fmin(eta,elevation_terrain_flines->element[kp])-elevation_bottom_flines->element[kp],0.0)*dsq->fine->grid->lines->element[kp]->length2d;
242  area_vert+=fmax(fmin(eta,elevation_terrain_flines->element[kp])-elevation_bottom_flines->element[kp],0.0)*dgrid->fine->lines->element[kp]->length2d; //dsq->fine->grid->lines->element[kp]->length2d;
243 
244  }
245  }
246 
247  return area_vert;
248 }
249 
250 
251 double vertical_area_surf(double eta, long j) {
264  long c,kp;
265  double area_vert=0.0;
266  long novalue=dgrid->novalue;
267  long nl_l=dgrid->fine->lines->nl;
268  long nh_l=dgrid->fine->lines->nh;
269 
270  // if ((j>dgrid->coarse->lines->nh) || (j<dgrid->coarse->lines->nl))
271  if ((j>dgrid->small_line_content->index->nh) || (j<dgrid->small_line_content->index->nl)) printf("Error in vertical_area_surf function: Line %ld does not exist [%ld,%ld] . \n",j,dgrid->fine->lines->nl,dgrid->fine->lines->nh);
272 
273  for(c=NL;c<=dgrid->small_line_content->index->element[j];c++) {
274 // for(c=1;c<=dsq->small_content_line->index->element[j];c++) {
275 // kp=dsq->small_content_line->element[j][c];
276  kp=dgrid->small_line_content->element[j][c];
277  if (((kp<nl_l) || (kp>nh_l)) && (kp!=novalue) ) { //
278  printf ("Error in vertical area lines %ld (of %ld) cannot belong to the big line %ld (c=%ld)!! ",kp,nh_l,j,c);
279  stop_execution();
280  } else if (kp!=novalue) {
281  // } else {
282  area_vert+=fmax(eta-elevation_terrain_flines->element[kp],0.0)*dgrid->fine->lines->element[kp]->length2d; //dsq->fine->grid->lines->element[kp]->length2d;
283 
284  }
285  }
286 
287  return area_vert;
288 }
289 
290 
291 
292 
293 
294 double wet_area(double eta,long i, double deta) {
311 // double vol1,vol2;
312 
313 // vol1=volume(eta-deta/2.0,i);
314 
315 // vol2=volume(eta+deta/2.0,i);
316 
317 // return (vol2-vol1)/deta;
319  long r,c,kp;
320  double aval=0.0;
321  double area=0.0;
322  long novalue=dgrid->novalue; // dsq->fine->novalue;
323 
324  if ((i<dgrid->small_polygon_content->index->nl) || (i>dgrid->small_polygon_content->index->nh)) printf("Error in wet_area function: cell %ld does not exist [%ld,%ld]. \n ",i,dgrid->small_polygon_content->index->nl,dgrid->small_polygon_content->index->nh);
325 
326  for(r=NL;r<=dgrid->small_polygon_content->index->element[i];r++){
327 // for(r=dsq->small_content_polygon->element[i]->nrl;r<=dsq->small_content_polygon->element[i]->nrh;r++){
328  // for(c=dsq->small_content_polygon->element[i]->ncl;c<=dsq->small_content_polygon->element[i]->nch;c++){
329  kp=dgrid->small_polygon_content->element[i][r]; //dsq->small_content_polygon->element[i]->element[r][c];
330  if (kp!=novalue){
331  if ((kp<elevation_bottom_fine->nl) || (kp>elevation_bottom_fine->nh)) {
332  aval=0.0;
333 
334  printf ("Error in wet_aera 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! novalue=%ld \n",kp,elevation_bottom_fine->nh,i,r,c,novalue);
335  stop_execution();
336  }else if (eta>=elevation_terrain_fine->element[kp]){
337  aval=dgrid->fine->polygons->element[kp]->area2D; //dsq->fine->grid->polygons->element[kp]->area2D;
338  } else if (eta>=elevation_bottom_fine->element[kp]) {
339  aval=dgrid->fine->polygons->element[kp]->area2D*porosity_fine->element[kp];//dsq->fine->grid->polygons->element[kp]->area2D*porosity_fine->element[kp];
340  } else {
341  aval=0.0;
342  }
343 
344  if (aval>0.0){
345  area+=aval;
346 
347  }
348  }
349  // }
350  }
351 
352  return area;
353 }
354 
355 
356 
357 double min_elevation(long i) {
369  double min_elevation=1.0e+5;
370  long r,c,kp;
371  long novalue=dgrid->novalue; // dsq->fine->novalue;
372 
373 
374  if ((i<dgrid->small_polygon_content->index->nl) || (i>dgrid->small_polygon_content->index->nh)) printf("Error in min_evation: cell %ld does not exist [%ld,%ld]. \n ",i,dgrid->small_polygon_content->index->nl,dgrid->small_polygon_content->index->nh);
375 
376  for(r=NL;r<=dgrid->small_polygon_content->index->element[i];r++){
377 // for(r=dsq->small_content_polygon->element[i]->nrl;r<=dsq->small_content_polygon->element[i]->nrh;r++){
378  // for(c=dsq->small_content_polygon->element[i]->ncl;c<=dsq->small_content_polygon->element[i]->nch;c++){
379  // kp=dsq->small_content_polygon->element[i]->element[r][c];
380  kp=dgrid->small_polygon_content->element[i][r]; //dsq->small_content_polygon->element[i]->element[r][c];
381  if (kp!=novalue){
382 
383  if ((kp<elevation_bottom_fine->nl) || (kp>elevation_bottom_fine->nh)) {
384 
385  printf ("Error in min_elevation 'small' polygon %ld (of %ld) cannot belong to the big polygon %ld (r=%ld,c=%ld)!! novalue=%ld \n",kp,elevation_bottom_fine->nh,i,r,c,novalue);
386  stop_execution();
387  }else {
388  min_elevation=fmin(elevation_bottom_fine->element[kp],min_elevation);
389  }
390 
391 
392  }
393  //}
394  }
395 
396  return min_elevation;
397 
398 }
399 
400 double q_discharge_from_outlet_subs_line(double eta, long j) {
413  long c,kp;
414  double q_outlet=0.0,q_outlet_surf=0.0;
415  double p=param->p_outlet;
416  double p_surf=param->p_outlet_surf;
417  long novalue=dgrid->novalue; // dsq->fine->novalue;
418  long nl_l=dgrid->fine->lines->nl;
419  long nh_l=dgrid->fine->lines->nh;
420 
421  if ((j>dgrid->small_line_content->index->nh) || (j<dgrid->small_line_content->index->nl)) printf("Error in q_discharge_from_outlet: Line %ld does not exist [%ld,%ld] . \n",j,dgrid->small_line_content->index->nl,dgrid->small_line_content->index->nh);
422 
423  for(c=NL;c<=dgrid->small_line_content->index->element[j];c++) {
424  // for(c=1;c<=dsq->small_content_line->index->element[j];c++) {
425  // kp=dsq->small_content_line->element[j][c];
426  kp=dgrid->small_line_content->element[j][c];
427  if (((kp<nl_l) || (kp>nh_l)) && (kp!=novalue) ) {
428  printf ("Error in q_discharge_from_outlet lines %ld (of %ld) cannot belong to the big line %ld (c=%ld)!! novalue=%ld \n",kp,nh_l,j,c,novalue);
429  stop_execution();
430  } else if (kp!=novalue && (outlet_coefficient->element[kp]!=0.0) && (outlet_coefficient!=NULL) && (outlet_coefficient_surf!=NULL)) {
431  // calcolo della portata
432  // area_vert+=fmax(eta-elevation_bottom_flines->element[kp],0.0)*dsq->fine->grid->lines->element[kp]->length2d;
433  q_outlet+=pow(fmax(fmin(eta,elevation_terrain_flines->element[kp])-elevation_bottom_flines->element[kp],0.0),p)*dgrid->fine->lines->element[kp]->length2d*outlet_coefficient->element[kp]; // dsq->fine->grid->lines->element[kp]->length2d
434  q_outlet_surf+=pow(fmax(eta-elevation_terrain_flines->element[kp],0.0),p_surf)*dgrid->fine->lines->element[kp]->length2d*outlet_coefficient_surf->element[kp]; // dsq->fine->grid->lines->element[kp]->length2d
435  }
436  }
437 
438  return q_outlet;
439 }
440 
441 double q_discharge_from_outlet_cell(double eta, long i) {
455  double q_out=0;
456  long j;
457  long nh_edges=dgrid->coarse->polygons->element[i]->edge_indices->nh; // dsq->big->grid->polygons->element[i]->edge_indices->nh;
458  long nl_edges=dgrid->coarse->polygons->element[i]->edge_indices->nl; //dsq->big->grid->polygons->element[i]->edge_indices->nl;
459  long kp,kl;
460  long kboundary=dgrid->coarse->boundary_indicator; //dsq->big->grid->boundary_indicator;
461 
462  for(j=nl_edges;j<=nh_edges;j++) {
463  kp=dgrid->coarse->links->element[i]->connections->element[j]; //dsq->big->grid->links->element[i]->connections->element[j];
464  kl=dgrid->coarse->polygons->element[i]->edge_indices->element[j]; //dsq->big->grid->polygons->element[i]->edge_indices->element[j];
465  if ((kp==kboundary) && (outlet_coefficient!=NULL)) q_out+=q_discharge_from_outlet_subs_line(eta,kl);
466  }
467 
468  return q_out;
469 }
470 
471