TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
b_solver.c
Go to the documentation of this file.
1 
27 #include "turtle.h"
28 #include "t_alloc.h"
29 #include "t_datamanipulation.h"
30 #include "t_utilities.h"
31 #include "rw_maps.h"
32 
33 #include "doublevector_utilities.h"
34 #include "preconditioned_conjugate_gradient.h"
35 #include "keywords_file_b.h"
36 
37 #include "geometry.h"
38 #include "g_raster2plvector.h"
39 #include "bigcells2.h"
40 #include "geometry2.h"
41 #include "b_utilities.h"
42 
43 #include "b_volumes.h"
44 #include "b_sources.h"
45 #include "b_v_advection.h"
46 #include "b_solver.h"
47 
48 #define MIN_T_DIAG_NEGATIVE_INDEX -20
49 
50 #define MIN_T_DIAG 1.0e-10
51 extern STRINGBIN *filenames;
52 extern DOUBLESQUARE_GRID *dsq;
53 extern DOUBLE_GRID *dgrid; // ec 20100413
54 extern DOUBLERASTER_MAP *draster;
55 extern char *wpath;
56 
57 extern PARAM *param;
58 extern FLAG *flag;
59 extern DOUBLEVECTOR *elevation_bottom_fine;
60 extern DOUBLEVECTOR *elevation_bottom_coarse;
61 extern DOUBLEVECTOR *elevation_bottom_flines;
62 extern DOUBLEVECTOR *elevation_bottom_bottom;
63 extern DOUBLEVECTOR *porosity_fine;
64 extern DOUBLEVECTOR *outlet_coefficient;
65 extern DOUBLEVECTOR *water_surface_elevation;
66 extern DOUBLEVECTOR *water_mass_error;
67 extern DOUBLEVECTOR *water_depth;
68 
69 extern DOUBLEVECTOR *surface_water_velocity;
70 extern DOUBLEVECTOR *F1_wet_vert_area;
71 DOUBLEVECTOR *eta_v, *sources, *dirichlet, *t_diagonal,*t_diagonal_no_dirichlet;
72 S_TIMES *s_times, *dirichlet_times;
73 //double Ks=10;
74 
75 double t_st_operator_element(long i,DOUBLEVECTOR *eta){
82  double cond_dirichlet=MAX_ELEVATION_VALUE;
83 
84  return (t_st_operator_element_subs(i,eta,cond_dirichlet)+t_st_advection_operator_element(i,eta,cond_dirichlet)); // ec 20100517 add double cond_dirichlet // corrected by ec 20100321 +1.0e-10*eta->element[i]);// (thirdtrial_2) corrected by ec 20100319
85 }
86 
87 double t_st_operator_element_no_dirichlet(long i,DOUBLEVECTOR *eta){
94  double cond_dirichlet=param->null_dirichlet;
95 
96  return (t_st_operator_element_subs(i,eta,cond_dirichlet)+t_st_advection_operator_element(i,eta,cond_dirichlet)); // ec 20100517 add double cond_dirichlet // corrected by ec 20100321 +1.0e-10*eta->element[i]);// (thirdtrial_2) corrected by ec 20100319
97 }
98 
99 double t_st_operator_element_subs(long i,DOUBLEVECTOR *eta,double cond_dirichlet) {
117  double sum=0.0;
118  double eta_previous;
119  long j,je,line_index,polygon2_index;
120  double eta1,eta2,t_coeff,coeff,dist;
121 
122 
123  for (j=dgrid->coarse->polygons->element[i]->edge_indices->nl;j<=dgrid->coarse->polygons->element[i]->edge_indices->nh;j++) {
124 
125  je=j;
126  if (dgrid->coarse->links->element[i]->connections->element[je]!=dgrid->coarse->boundary_indicator) {
127 
128  line_index=dgrid->coarse->polygons->element[i]->edge_indices->element[je];
129  polygon2_index=dgrid->coarse->links->element[i]->connections->element[je];
130  if (dirichlet->element[polygon2_index]<=cond_dirichlet) { // ec 20100517
131  dist=dgrid->coarse->links->element[i]->d_connections->element[je];
132  eta1=water_surface_elevation->element[i];
133  eta2=water_surface_elevation->element[polygon2_index];
134  eta_previous=water_surface_elevation_mean(eta1,eta2);
135  t_coeff=param->Ks*vertical_area_subs(eta_previous,line_index); // MODIFY HERE param->p_outlet
136  // t_coeff=param->Ks*fmax(vertical_area_subs(eta1,line_index),vertical_area_subs(eta2,line_index));
137  // t_coeff=param->Ks*(vertical_area(eta1,line_index)+vertical_area(eta2,line_index))/2.0; /* modified on 28 July 2008 */
138  coeff=-param->dt*t_coeff*(eta->element[polygon2_index]-eta->element[i])/dist;
139  sum=sum+coeff;
140  }
141  }
142  }
143 
144  return sum;
145 }
146 /*
147 int T_st_operator(DOUBLEVECTOR *y, DOUBLEVECTOR *eta) {
148  *!
149  * \author Emanuele Cordano
150  * \date 19 April 2009 modified on 29 june 2009
151  *
152  * \param y - (DOUBLEVECTOR *) output
153  * \param eta - (DOUBLEVECTOR *) eta
154  *
155  *
156  *
157 
158  long i;
159 // long je; modified by Emanuele Cordano on 24 Jun 2009
160 // double sum,coeff;
161 // long i_cell=(eta->nh+1)/2;
162 // double t_coeff;
163 // double dist,eta1,eta2;
164 // long line_index,polygon2_index;
165 
166  if ((y->nh!=eta->nh) ||(eta->nh!=dgrid->coare->big->polygons->nh)) {
167  printf("Error in T_st_operator function y and etahas different size [%ld,%ld expected: %ld!\n]",y->nh,eta->nh,dsq->big->grid->polygons->nh);
168  return -1;
169  }
170 
171  for (i=eta->nl;i<=eta->nh;i++) {
172  y->element[i]=y->element[i]+t_st_operator_element(i,eta);
173  }
174 
175  return 0;
176 
177 }
178 */
179 
180 /*
181 int wet_area_operator(DOUBLEVECTOR *y, DOUBLEVECTOR *x) {
182  *!
183  *\author Emanuele Cordano
184  *\date 19 April 2009
185  *
186  * * \param y - (DOUBLEVECTOR *) output
187  * \param x - (DOUBLEVECTOR *) input
188  *
189  *
190  long i;
191  double area;
192 
193 
194  if ((y->nh!=x->nh) ||(x->nh!=dsq->big->grid->polygons->nh)) {
195  printf("Error in wet_area_operator function y and etahas different size [%ld,%ld expected: %ld!\n]",y->nh,x->nh,dsq->big->grid->polygons->nh);
196  return -1;
197  }
198 
199  for (i=x->nl;i<=x->nh;i++) {
200 // vol=volume(eta_v->element[i],i);
201  area=wet_area(eta_v->element[i],i,param->deta);
202  y->element[i]=y->element[i]+area*x->element[i];
203  }
204 
205  return 0;
206 }
207 */
208 /*int volume_operator(DOUBLEVECTOR *y, DOUBLEVECTOR *eta) {
209  *!
210  *\author Emanuele Cordano
211  *\date 19 April 2009
212  *
213  * * \param y - (DOUBLEVECTOR *) output
214  * \param eta - (DOUBLEVECTOR *) eta
215  *
216  *
217  long i;
218  double vol;
219 
220 
221  if ((y->nh!=eta->nh) ||(eta->nh!=dsq->big->grid->polygons->nh)) {
222  printf("Error in volume_operator function y and eta have different size [%ld,%ld expected: %ld!\n]",y->nh,eta->nh,dsq->big->grid->polygons->nh);
223  return -1;
224  }
225 
226  for (i=eta->nl;i<=eta->nh;i++) {
227  vol=0.0;
228 
229  vol=volume(eta->element[i],i);
230 // printf("\nDebug y[%ld]=%le vol=%le area=%lf \n",i,y->element[i],vol,area);
231  y->element[i]=y->element[i]+vol;
232 // printf("\nDebug 2 y[%ld]=%le vol=%le \n",i,y->element[i],vol);
233 
234  }
235 
236  return 0;
237 }
238 */
239 /*
240 int volume_operator_minus(DOUBLEVECTOR *y, DOUBLEVECTOR *eta) {
241  *!
242  *\author Emanuele Cordano
243  *\date 19 April 2009 modified on 1 June 2009
244  *
245  * \param y - (DOUBLEVECTOR *) output
246  * \param eta - (DOUBLEVECTOR *) eta
247  *
248  *
249  *
250  long i;
251  double vol,area;
252 
253 
254  if ((y->nh!=eta->nh) ||(eta->nh!=dsq->big->grid->polygons->nh)) {
255  printf("Error in volume_operator_minus function y and eta have different size [%ld,%ld expected: %ld!\n]",y->nh,eta->nh,dsq->big->grid->polygons->nh);
256  return -1;
257  }
258  double t_inv=param->t-param->dt/2.0;
259  //printf("t_inv= %lf \n",t_inv);
260  get_sources(t_inv,sources);
261  for (i=eta->nl;i<=eta->nh;i++) {
262  vol=volume(eta->element[i],i)+sources->element[i]*dsq->big->grid->polygons->element[i]->area2D*param->dt;
263 // vol=0.0;
264 // printf("\nDebug y[%ld]=%le vol=%le \n",i,y->element[i],vol);
265  y->element[i]=y->element[i]-vol;
266 // printf("\nDebug 2 y[%ld]=%le vol=%le \n",i,y->element[i],vol);
267 
268  }
269 
270  return 0;
271 }
272 */
273 
274 double b_smatrix_element (long i,DOUBLEVECTOR *x){
275 
291  if ((t_diagonal->element[i]!=MIN_T_DIAG_NEGATIVE_INDEX) && (dirichlet->element[i]<=param->null_dirichlet)) { // ec 10100323
292  return (t_st_operator_element_no_dirichlet(i,x)+wet_area(eta_v->element[i],i,param->deta)*x->element[i]-t_diagonal_no_dirichlet->element[i]*x->element[i]+t_diagonal->element[i]*x->element[i]); //corrected by ec 20100323 //corrected by ec 20100321 for thridtrial_3
293  } else if (dirichlet->element[i]>param->null_dirichlet) {
294 
295  return (wet_area(eta_v->element[i],i,param->deta)*x->element[i]);
296  } else {
297  return (x->element[i]);
298  }
299 
300 }
301 
302 /*
303 
304 int b_smatrix(DOUBLEVECTOR *y, DOUBLEVECTOR *x) {
305  *!
306  * \author Emanuele Cordano
307  * \date 29 june 2009
308  *
309  * \param y - (DOUBLEVECTOR *) output
310  * \param x - (DOUBLEVECTOR *) input
311  *
312  * \details It fills the vector y using the function b_smatrix_element() for each element of y.
313  *
314 
315  long i;
316 
317  if ((y->nh!=x->nh) ||(x->nh!=dsq->big->grid->polygons->nh)) {
318  printf("Error in b_smatrix function y and x has different size [%ld,%ld expected: %ld!\n]",y->nh,x->nh,dsq->big->grid->polygons->nh);
319  return -1;
320  }
321 
322  for (i=x->nl;i<=x->nh;i++) {
323  y->element[i]=b_smatrix_element(i,x);
324  }
325 
326  return 0;
327 
328 }
329 */
330 //int b_knownterm_0(DOUBLEVECTOR *be0)
331 /*int b_knownterm(DOUBLEVECTOR *be) {
332  *!
333  *
334  *\author Emanuele Cordano
335  *\date April 2009
336  *
337  *\param be : (DOUBLEVECTOR *)
338  *
339  * /DA METTERE A POSTO!!!
340  *
341 
342 
343  int l,s,m,p;
344 
345  l=zeros(be);
346 // print_doublevector_elements(eta_v,100);
347 // print_doublevector_elements(water_surface_elevation,100);
348 
349  m=volume_operator(be,eta_v);
350 
351  s=volume_operator_minus(be,water_surface_elevation);
352 
353  p=T_st_operator(be,eta_v);
354 
355 
356  if ( (p==0) && (m==0) && (s==0)) {
357  // printf("b_knownterm exit 0 \n");
358  return 0;
359  } else {
360  printf("b_knownterm exit -1 \n");
361  }
362 
363  return -1;
364 
365 }
366 */
367 long Newton_convergence(DOUBLEVECTOR *x_temp,DOUBLEVECTOR *be, DOUBLEVECTOR *be0) { //doublee ?????
387 /*
388  * The introduced vector is $\mathbf{Ws}( ^{m}\eta^{n+1})$ is the vector of the wet area for each cell, in fact from (\ref{def_Vw}), it is:
389 \begin{equation} \label{def_Vw}
390  Ws_{i}(\psi_i) =\frac {\mathrm{d} V_{i}(\eta_i)}{\mathrm{d} \eta_i}
391 \end{equation}
392 The solution of (\ref{def_solver}) is unique and can be implemented \citep{Casulli2008} .
393  */
394  char *function_name="Newton_convergence";
395  double a=0;
396  double max_b,max_x;
397  long i;
398  double x_temp_max,x_temp_adm=param->x_temp_adm;
399  long kkm,kkmm=0;
400  long time_0=0;
401  long time_1=0;
402  double elapsed;
403 #ifdef WRITE_ITERATION_NUMBER
404  char *LOG_FILE_NAME=join_strings(wpath,"boussinesq_iterations.log");
405  long kcnt=0;
406  FILE *FILE_LOG;
407  double volume_term=0.0,be_term=0.0;
408  FILE_LOG=fopen(LOG_FILE_NAME,"a");
409 #endif
410  double initial_residual;
411  long newton_cnt=0;
412  /* initialization of known term */
413 
414  for(i=eta_v->nl;i<=eta_v->nh;i++) {
415  if (dirichlet->element[i]<=param->null_dirichlet) {
416  volume_term=volume(eta_v->element[i],i);
417  be_term=be0->element[i];
418  be->element[i]=volume_term+t_st_operator_element(i,eta_v)-be_term;
419  if ((t_diagonal->element[i]<=MIN_T_DIAG) && (volume_term<=MIN_T_DIAG) && (be_term<=MIN_T_DIAG)) { //ec 20100323
420  t_diagonal->element[i]=MIN_T_DIAG_NEGATIVE_INDEX;
421  be->element[i]=eta_v->element[i]-water_surface_elevation->element[i];
422  } //end ec 20100323 PUT AN ELSE CONDITION HERE!!!
423  } else {
424  be->element[i]=volume(eta_v->element[i],i)-volume(dirichlet->element[i],i); //eta_v->element[i]-dirichlet->element[i];
425  // printf ("be[%ld]=%lf\n",i,be->element[i]);
426  // stop_execution();Initial residual
427  }
428 
429  }
430  initial_residual=max_doublevector(be);
431  //printf ("Initial residual: %le",initial_residual);
432  //stop_execution();
433  /* END initialization of known term */
434  double newton_toll=param->max_errors;
435 
436  do {
437  for (i=eta_v->nl;i<=eta_v->nh;i++) {
438  // if (dirichlet->element[i]>param->null_dirichlet) printf("eta[%ld]=%lf \n",i,eta_v->element[i]);
439  //eta_v->element[i]=dirichlet->element[i];
440  } // added_by_ec_20100517
441  // zeros(x_temp);//fvector
442  for(i=eta_v->nl;i<=eta_v->nh;i++) {
443  if (dirichlet->element[i]<=param->null_dirichlet) {
444  be->element[i]=volume(eta_v->element[i],i)+t_st_operator_element(i,eta_v)-be0->element[i];//here_is_the errorr ec 20100518
445  } else {
446  be->element[i]=volume(eta_v->element[i],i)-volume(dirichlet->element[i],i); //eta_v->element[i]-dirichlet->element[i];
447  // printf ("be[%ld]=%lf,eta_v->element[%ld]=%lf \n",i,be->element[i],i,eta_v->element[i]);
448  // stop_execution();
449  }
450 
451  }
452 
453  // b_knownterm(be);
454 
455  kkm=0;
456  // kkmm=kkmm+1;
457  newton_cnt++;
458  //if (max_doublevector(be)>param->max_errors) {
459  if (max_doublevector(be)>newton_toll) {
460  time_0=clock();
461  // kkm=conjugate_gradient_search(kkm,param->max_errors,x_temp,be,b_smatrix); /* please uncomment for the unconditioned CG solver */
462  //kkm=jacobi_preconditioned_conjugate_gradient_search(kkm,param->max_errors, x_temp,be,(*b_smatrix_element));
463  kkm=jacobi_preconditioned_conjugate_gradient_search(kkm,1e-6, x_temp,be,(*b_smatrix_element));
464  time_1=clock();
465  elapsed=(double)((time_1-time_0)/CLOCKS_PER_SEC);
466 // max_residual printf ("Number of reiterations %ld (execution time : %lf )\n",kkm,elapsed);
467 #ifdef WRITE_ITERATION_NUMBER
468  fprintf(FILE_LOG,"%ld,",kkm);
469 #endif
470  if (newton_cnt>1) {
471  for (i=x_temp->nl;i<=x_temp->nh;i++) {
472  if (x_temp->element[i]<0.0) {
473  // printf("WARNING in Newton_covergence: x_temp[%ld] is negative (%le) (be=%le) (eta0=%lf, eta1=%lf) at %ld iteration !! \n",i,x_temp->element[i],be->element[i],water_surface_elevation->element[i],eta_v->element[i],kkmm);
474  //stop_execution();
475  // controllare errato!!!
476  x_temp->element[i]=0.0;
477  }
478  }
479 
480  }
481  for(i=x_temp->nl;i<=x_temp->nh;i++) {
482 
483  eta_v->element[i]=eta_v->element[i]-x_temp->element[i];
484  //if (dirichlet->element[i]>param->null_dirichlet) eta_v->element[i]=dirichlet->element[i];
485 
486  }
487  //printf ("xtemo %le",max_doublevector(x_temp));
488  DOUBLEVECTOR* residual=new_doublevector(eta_v->nh);
489 // double sum_tmp=0.0;
490  for(i=eta_v->nl;i<=eta_v->nh;i++) {
491  // residual->element[i]=be->element[i]-b_smatrix_element(i,x_temp);
492  if (dirichlet->element[i]<=param->null_dirichlet) {
493  residual->element[i]=volume(eta_v->element[i],i)-be0->element[i]+t_st_operator_element(i,eta_v);
494  // residual->element[i]=volume(eta_v->element[i],i)-be0->element[i]+t_st_operator_element_no_dirichlet(i,eta_v);
495  // if ((be0->element[i]<=0.0) && (t_diagonal->element[i]<1.0e-10)) {
496  // be->element[i]=0.0;
497  // t_diagonal->element[i]=;
498  // }
499  } else {
500  residual->element[i]=volume(eta_v->element[i],i)-volume(dirichlet->element[i],i);
501  // x_temp->element[i];
502  //double residual2=eta_v->element[i]-dirichlet->element[i];
503  // printf("residual[%ld]=%le residual2=%le \n",i,residual->element[i],residual2);
504  // stop_execution();
505  }
506  // sum_tmp=sum_tmp+residual->element[i]+be0->element[i];
507  // if (sum_tmp<=0.0) {
508  // printf("res[%ld]=%le be[%ld]=%le sum=%le\n",i,residual->element[i],i,be0->element[i],sum_tmp);
509 
510  // stop_execution();
511  // }
512  }
513  // if (sum_tmp<=0.0) {
514  // printf("sum=%le \n",sum_tmp);
515  //if (sum_tmp<=0.0) {printf("sum_tmp=%f",sum_tmp);stop_execution();}
516  // }
517  // free_doublevector(residual);
518  max_b=max_doublevector(residual);
519  free_doublevector(residual);
520  // ec 20100324 time 00:48
521  max_x=max_doublevector(x_temp);
522  if (max_x==0.0) {
523 
524  max_b=0.0;
525  }
526  // end ec 20100324 time 00:48
527  } else {
528  // kkm=-1;
529  // printf ("Mass balance verified with error=%le (max %le)\n",max_doublevector(be),param->max_errors);
530  max_b=0.0;
531  }
532 
533  // x_temp_max=max_doublevector(x_temp);
534 
535 
536  // printf("max_residual=%le adm=%le \n",max_b,newton_toll);
537 
538  }while (max_b>newton_toll); //(max_b>1.0e-6);
539 
540  //while (max_b>param->max_errors);//((kkm>0) || (x_temp_max>x_temp_adm));
541 #ifdef WRITE_ITERATION_NUMBER
542  fprintf(FILE_LOG," %ld \n",newton_cnt);
543  fclose(FILE_LOG);
544 #endif
545  return newton_cnt; // kkmm;
546 }
547 
548 int time_loop(short print,int (*write_output)(void *v1, void *v2)){
567  double t=param->t;
568  double t_start=param->t_start;
569  double t_end=param->t_end;
570  double dt_print=param->dt_print;
571  double dt=param->dt;
572  double t_inv;
573  OUTPUT_FILENAMES *fn;
574  long i;
575  int s;
576  FILE *fd;
577 
578  char SSSS[ ]={"SSSS"}; /* string containing th number of the printed map */
579  long kks,lt;
580  DOUBLEVECTOR *x,*be,*be0;
581  x=new_doublevector(water_surface_elevation->nh);
582  for (i=x->nl;i<=x->nh;i++) { // ec 20110419
583  x->element[i]=0.0;
584  } // end ec 20110419
585  be=new_doublevector(water_surface_elevation->nh);
586  be0=new_doublevector(water_surface_elevation->nh);
587  fn=new_output_filenames(print);
588 
589  eta_v=new_doublevector(water_surface_elevation->nh);
590  sources=new_doublevector(water_surface_elevation->nh);
591  dirichlet=new_doublevector(water_surface_elevation->nh);
592  t_diagonal=new_doublevector(water_surface_elevation->nh); // ec 20100322
593  t_diagonal_no_dirichlet=new_doublevector(water_surface_elevation->nh); // ec 20100518
594  for (i=t_diagonal->nl;i<=t_diagonal->nh;i++) {// ec 20100322
595  t_diagonal->element[i]=0.0;
596  t_diagonal_no_dirichlet->element[i]=0.0;
597  }// end ec 20100322
598 
599  fd=fopen(filenames->element[I_TIMES_FT_FILE]+1,"r");
600  if (fd!=NULL) {
601  fclose(fd);
602  s_times=get_s_times(filenames->element[I_TIMES_FT_FILE]+1,print);
603  t=param->t_start;
604  get_sources(t,sources);
605  } else {
606  for (i=sources->nl;i<=sources->nh;i++) {
607  sources->element[i]=0.0;
608  }
609  }
610  fd=fopen(filenames->element[I_DIRICHLETTIMES_FT_FILE]+1,"r");
611 
612  if (fd!=NULL) {
613  fclose(fd);
614  printf("cxxx\n");
615  dirichlet_times=get_s_times(filenames->element[I_DIRICHLETTIMES_FT_FILE]+1,print);
616  t=param->t_start;
617  get_dirichletsnode(t,dirichlet);
618  } else {
619  printf("%s\n",filenames->element[I_DIRICHLETTIMES_FT_FILE]+1);
620  for (i=dirichlet->nl;i<=dirichlet->nh;i++) {
621  dirichlet->element[i]=param->null_dirichlet;
622  }
623  }
624 
625 
626  printf("water_table_surface_elevation cells:%ld\n",water_surface_elevation->nh);
627 
628 
629  while ((t>=t_start) && (t<=t_end)) {
630  t=t+dt;
631  param->t=t;
632  t_inv=t-dt/2.0;
633 
634  if (s_times!=NULL) get_sources(t_inv,sources);
635  if (dirichlet_times!=NULL) get_dirichletsnode(t,dirichlet);
636 
637  update_F1_wet_vert_area();
638 
639  for (i=water_surface_elevation->nl;i<=water_surface_elevation->nh;i++){
640  eta_v->element[i]=water_surface_elevation->element[i]; /* modified by Emanuele Cordano on may 20, 2009 */
641  // printf(" %le",sources->element[i]);
642  // stop_execution();
643  if (sources->element[i]>0.0) eta_v->element[i]=elevation_bottom_bottom->element[i]+sources->element[i]*dt; /* modified by Emanuele Cordano on may 11, 2010 */
644 
645  if (dirichlet->element[i]>param->null_dirichlet) eta_v->element[i]=dirichlet->element[i];
646 
647 
648  if (water_depth==NULL) {
649  be0->element[i]=volume(water_surface_elevation->element[i],i)+sources->element[i]*dgrid->coarse->polygons->element[i]->area2D*dt-q_discharge_from_outlet_cell(water_surface_elevation->element[i],i)*dt;
650 
651  }else {
652  be0->element[i]=(water_depth->element[i]+sources->element[i]*dt)*dgrid->coarse->polygons->element[i]->area2D-q_discharge_from_outlet_cell(water_surface_elevation->element[i],i)*dt;
653 
654  }
655  /* known term for advection in surface flow */
656 
657  be0->element[i]=be0->element[i]+b_advection(i); // modified by Emanuele Cordano on 29 Oct 2009
658 
659  /* added by Emanuele Cordano */
660  }
661  //printf ("stop here11");
662  //stop_execution();
663  if (water_depth!=NULL) {
664  printf("Warning: volumes at the nodes were calculated according to water depth value. Water depth values are now deleting \n");
665  free_doublevector(water_depth);
666  water_depth=NULL;
667 
668  }
669  get_diagonal(t_diagonal,t_st_operator_element); //get digonal of t_st_element ec-20100323
670  get_diagonal(t_diagonal_no_dirichlet,t_st_operator_element_no_dirichlet); //get digonal of t_st_element_no_dirichlet ec-20100518
671  kks=Newton_convergence(x,be,be0);
672 
673  /* update of wet_area_surface_velocity */
674 
675  update_velocity(eta_v); // added by Emanuele Cordano on 29 Oct 2009
676  for (i=water_surface_elevation->nl;i<=water_surface_elevation->nh;i++){
677 
678  water_surface_elevation->element[i]=eta_v->element[i];
679  water_mass_error->element[i]=be->element[i];
680  }
681  if (print==1) printf("Time_loop function: time %le [%le,%le] executed with %ld iterations \n",t,t_start,t_end,kks);
682  if ((long)t%(long)dt_print==0) {
683  //if (print==1) printf("Time_loop function: time %le [%le,%le] executed with %ld iterations \n",t,t_start,t_end,kks);
684 
685  lt=(long)(t/dt_print);
686  // SSSS="0000";
687 
688  // printf(" t= %ld %le \n ",lt,t);
689 // stop_execution();
690  write_suffix( SSSS, (long)t/(long)dt_print,0);
691  s=(*write_output)((void *)fn,(void *)SSSS);
692  if (s!=0) printf("Error in time_loop: results at %s was not correctly written (exit %d)",SSSS,s);
693  }
694 
695 
696  }
697 
698  if (dirichlet_times!=NULL) free_s_times(dirichlet_times);
699  if (s_times!=NULL) free_s_times(s_times);
700  free_doublevector(eta_v);
701  free_doublevector(be);
702  free_doublevector(be0);
703  free_doublevector(x);
704  free_doublevector(dirichlet);
705  free_doublevector(sources);
706  free_doublevector(t_diagonal);
707  free_OUTPUT_FILENAMES(fn);
708  if (print==1) printf ("Function time_loop was successfully executed!! ");
709  //free(SSSS);
710  return 0;
711 }
712 
713 
714 double water_surface_elevation_mean(double eta1,double eta2) {
722  double val=eta1;
723 
724  if (flag->arithmetic_mean0==1) {
725  val=(eta1+eta2)/2.0;
726  } else {
727  val=fmax(eta1,eta2);
728  }
729 
730  return val;
731 }
732 
733 
734 //int write_map_results(void *file_results, void *file_error,void *argument3)
735 int write_map_results(void *output, void *SSSS) {
736  /*
737  *\author Emanuele Cordano
738  *\date 21 April 2009
739  *
740  *\param (void *) file_results - name of files where results are plotted
741  *\param (void *) file_error - name of files where errors are plotted
742  *\param (void *) argument3 - pointer which points to information about the time step
743  *
744  *
745  *
746  */
747  //char *filename_results=(char*)file_results;
748  //char *filename_error=(char*)file_error;
749  //char *SSSS;
750 
751 
752  OUTPUT_FILENAMES *fu=(OUTPUT_FILENAMES *)output;
753  int l,s,g;
754  char *results;
755 
756  char *error;
757  char *velocity_name;
758 
759  results=join_strings(fu->file_result,(char *)SSSS);
760  error=join_strings(fu->file_error,(char *)SSSS);
761  velocity_name=join_strings(results,"_velocity.ft");
762 
763  g=write_doublevector_in_a_ascii_ft_format_single_file(velocity_name,surface_water_velocity);
764 
765  l=write_raster_from_doublevector(results,water_surface_elevation,draster->coarse->UV,dsq->big->indices_pixel,draster->coarse->layer[DTM_MASK]);
766  s=write_raster_from_doublevector(error,water_mass_error,draster->coarse->UV,dsq->big->indices_pixel,draster->coarse->layer[DTM_MASK]);
767  if ((l==0) && (s==0)) {
768  free(results);
769  free(error);
770  return 0;
771  } else {
772  printf ("Error in write_map_results l=%d s=%d \n ",l,s);
773  }
774  free(results);
775  free(error);
776  free(velocity_name);
777  // free(SSSS);
778  return -1;
779 }
780 /* output are seved and then written in ascii format files through the following functions. */
781 
782 OUTPUT_FILENAMES *new_output_filenames(short print){
783 /*
784  * \author Emanuele Cordano
785  * \data April 2006
786  *
787  */
788  OUTPUT_FILENAMES *fn;
789 
790  fn=(OUTPUT_FILENAMES *)malloc((size_t)(sizeof(OUTPUT_FILENAMES)));
791  if (!fn) t_error("Fn (OUTPUT_FILENAMES) was not allocated");
792 
793  fn->file_result=copy_string(filenames->element[O_COARSEMAP_WATERSURFACE_ELEVATION]+1);
794  fn->file_error=copy_string(filenames->element[O_COARSEMAP_WATERMASS_ERROR]+1);
795 
796  if (print==1) printf("OUTPUT_FILENAMES was successfully initialized file_results: %s file_error: %s !! \n",fn->file_result,fn->file_error);
797 
798  return fn;
799 }
800 
801 void free_OUTPUT_FILENAMES(OUTPUT_FILENAMES *fn){
802 /*
803  *
804  * \author Emanuele Cordano
805  * \date April 2006
806  *
807  */
808 
809 //if (fn->SSSS!=NULL) free(fn->SSSS);
810 if (fn->file_error!=NULL) free(fn->file_error);
811 if (fn->file_result!=NULL) free(fn->file_result);
812 
813 free(fn);
814 }
815 
816