TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
b_sources.c
Go to the documentation of this file.
1 
27 #include "turtle.h"
28 
29 #include "t_alloc.h"
30 #include "t_io.h"
31 #include "t_datamanipulation.h"
32 #include "t_utilities.h"
33 #include "rw_maps.h"
34 #include "geometry.h"
35 #include "g_raster2plvector.h"
36 #include "bigcells2.h"
37 #include "geometry2.h"
38 #include "b_utilities.h"
39 #include "b_solver.h"
40 #include "b_volumes.h"
41 #include "doublevector_utilities.h"
42 #include "preconditioned_conjugate_gradient.h"
43 #include "keywords_file_b.h"
44 #include "b_sources.h"
45 
46 extern STRINGBIN *filenames;
47 extern DOUBLESQUARE_GRID *dsq;
48 extern DOUBLE_GRID *dgrid; // ec 20100413
49 extern DOUBLERASTER_MAP *draster;
50 
51 
52 extern PARAM *param;
53 extern FLAG *flag;
54 extern DOUBLEVECTOR *elevation_bottom_fine;
55 extern DOUBLEVECTOR *elevation_bottom_coarse;
56 extern DOUBLEVECTOR *elevation_bottom_flines;
57 extern DOUBLEVECTOR *elevation_bottom_bottom;
58 extern DOUBLEVECTOR *porosity_fine;
59 extern DOUBLEVECTOR *water_surface_elevation;
60 extern DOUBLEVECTOR *water_mass_error;
62 extern S_TIMES *s_times, *dirichlet_times;
63 //fra una funzione di lettura tempi suffissi;
64 
65 
66 
67 int get_sources(double t,DOUBLEVECTOR *sources){
80  static double t0;
81  static double t1;
82  static long k;
83  static DOUBLEVECTOR *s0;
84  static DOUBLEVECTOR *s1;
85  char *map0=NULL;
86  char *map1=NULL; /* modified by Emanuele Cordano on October 2009 */
87  long j;
88 
89 
90 
91  if ((t<=s_times->times->element[s_times->times->nl]) || (s1==NULL) || (s0==NULL)) {
92 
93  k=s_times->times->nl;
94  t0=s_times->times->element[k];
95  t1=s_times->times->element[k+1];
96  map0=join_strings(filenames->element[I_SOURCEMAPSERIES_COARSE]+1,s_times->s_suffixes->element[k]+1);
97  if (s0==NULL) s0=read_doublevector_from_raster(A_FLAG,map0,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
98  map1=join_strings(filenames->element[I_SOURCEMAPSERIES_COARSE]+1,s_times->s_suffixes->element[k+1]+1);
99  if (s1==NULL) s1=read_doublevector_from_raster(A_FLAG,map1,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
100 
101  } else while ((t>t1) && (k<s_times->times->nh)) {
102 
103 
104  copy_doublevector(s1,s0);
105  t0=t1;
106  free_doublevector(s1);
107  k++;
108  free(map1);
109  map1=join_strings(filenames->element[I_SOURCEMAPSERIES_COARSE]+1,s_times->s_suffixes->element[k+1]+1);
110  s1=read_doublevector_from_raster(A_FLAG,map1,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
111  t1=s_times->times->element[k+1];
112  }
113 
114  if ((s1->nh!=s0->nh) || (sources->nh!=s0->nh)) printf("Error in get_sources: vector s0 [%ld]and s1 [%ld] and s[%ld] have different size!! \n",s1->nh,s0->nh,sources->nh);
115  for(j=sources->nl;j<=sources->nh;j++) {
116  sources->element[j]=(t-t0)/(t1-t0)*s1->element[j]+(t1-t)/(t1-t0)*s0->element[j];
117 
118  }
119 
120  if (map0!=NULL) free(map0);
121  if (map1!=NULL) free(map1);
122 
123  return 0;
124 }
125 
126 
127 S_TIMES *get_s_times(char *filename,short print) {
136  S_TIMES *s_t;
137  FILE *f;
138 // leggere s_times
139  s_t=(S_TIMES *)malloc(sizeof(S_TIMES));
140  if (!s_t) t_error("In get_s_times: s_t was not allocated! ");
141 
142  f=t_fopen(filename,"r");
143  int ifile=read_index(f,print);
144  s_t->times=read_doublearray(f,print);
145  s_t->s_suffixes=read_stringarray(f,print);
146 
147  t_fclose(f);
148 
149  /* verifies*/
150 
151  if (s_t->times->nh!=s_t->s_suffixes->index->nh) {
152  printf("Error in get_sorces: times and s_suffixes has different size %ld ad %ld respectively!! \n",s_t->times->nh,s_t->s_suffixes->index->nh);
153  }
154  if (s_t->times->nh<2) {
155  printf("Error in get_sorces: times %ld are less than two values!! \n",s_t->times->nh);
156  }
157 
158  long k;
159  for (k=s_t->times->nl;k<s_t->times->nh;k++) {
160  if (s_t->times->element[k+1]<s_t->times->element[k]) printf("Error in get_sources: times[%ld]=%le is lower than times[%ld]=%le !! \n",k+1,s_t->times->element[k+1],k,s_t->times->element[k]);
161  }
162  return s_t;
163 
164 }
165 
166 void free_s_times(S_TIMES* s_t) {
173  if (s_t->times!=NULL) free_doublevector(s_t->times);
174  if (s_t->s_suffixes!=NULL) free_stringbin(s_t->s_suffixes);
175 
176  if (s_t!=NULL) free(s_t);
177 
178 }
179 
180 
181 
182 int get_dirichletsnode(double t,DOUBLEVECTOR *dirichlet){
195  static double t0,tp0;
196  static double t1,tp1;
197  static long k;
198  static DOUBLEVECTOR *d0;
199  static DOUBLEVECTOR *d1;
200  char *map0,*map1;
201  double tp;
202  long j;
203 
204  if ((t<=dirichlet_times->times->element[dirichlet_times->times->nl]) || (d1==NULL) || (d0==NULL)) {
205 // printf("h %ld t=%lf",dirichlet_times->times->nl,t);
206 
207  k=dirichlet_times->times->nl;
208  t0=dirichlet_times->times->element[k];
209  t1=dirichlet_times->times->element[k+1];
210  map0=join_strings(filenames->element[I_DIRICHLETMAPSERIES_COARSE]+1,dirichlet_times->s_suffixes->element[k]+1);
211  if (d0==NULL) d0=read_doublevector_from_raster(A_FLAG,map0,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
212  map1=join_strings(filenames->element[I_DIRICHLETMAPSERIES_COARSE]+1,dirichlet_times->s_suffixes->element[k+1]+1);
213  if (d1==NULL) d1=read_doublevector_from_raster(A_FLAG,map1,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
214  tp0=pow(t0,param->exp_dirichlet);
215  tp1=pow(t1,param->exp_dirichlet);
216  } else while ((t>t1) && (k<dirichlet_times->times->nh)) {
217 
218 
219  copy_doublevector(d1,d0);
220  t0=t1;
221  free_doublevector(d1);
222  k++;
223  map1=join_strings(filenames->element[I_DIRICHLETMAPSERIES_COARSE]+1,dirichlet_times->s_suffixes->element[k+1]+1);
224  d1=read_doublevector_from_raster(A_FLAG,map1,draster->coarse->layer[DTM_MASK],draster->coarse->UV,dsq->big->indices_pixel);
225  t1=dirichlet_times->times->element[k+1];
226  tp0=pow(t0,param->exp_dirichlet);
227  tp1=pow(t1,param->exp_dirichlet);
228  //tp0=t0; //ec_20100511
229  //tp1=t1; //ec_20100511
230 
231  }
232 
233  if ((d1->nh!=d0->nh) || (dirichlet->nh!=d0->nh)) printf("Error in get_sources: vector d0 [%ld]and d1 [%ld] and dirichlet[%ld] have different size!! \n",d1->nh,d0->nh,dirichlet->nh);
234  for(j=dirichlet->nl;j<=dirichlet->nh;j++) {
235  if ((d0->element[j]<=param->null_dirichlet) || (d1->element[j]<=param->null_dirichlet)) {
236 
237  dirichlet->element[j]=param->null_dirichlet;
238  } else {
239 
240  tp=pow(t,param->exp_dirichlet);
241  // tp=t; //ec_20100511
242  dirichlet->element[j]=(tp-tp0)/(tp1-tp0)*d1->element[j]+(tp1-tp)/(tp1-tp0)*d0->element[j];
243  // printf("tp0=%le tp1=%le t1=%le tp=%le val0=%le val1=%le val=%le \n",tp0,tp1,t1,tp,d0->element[j],d1->element[j],dirichlet->element[j]);
244  // stop_execution();
245  }
246 
247 // printf("s[%ld]=%lf (system %lf) \n",j,sources->element[j],t);
248  }
249 
250  return 0;
251 }
252 
253