source: issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp@ 18473

Last change on this file since 18473 was 18473, checked in by abuzzi, 11 years ago

MODIF: vfsa that works
File size: 22.2 KB
Line 
1#include <iostream>
2#include <cmath>
3#include <fstream>
4#include <string>
5#include <cstdio>
6#include <cstdlib>
7#include <time.h>
8#include <cassert>
9#include <gsl/gsl_multifit.h>
10#include <pthread.h>
11
12using namespace std;
13
14class Matrix{/*{{{*/
15 private:
16 int M; /*Number of lines */
17 int N; /*Number if Columns */
18 double *values;
19 public:
20 Matrix(int m_in,int n_in){/*{{{*/
21 this->M = m_in;
22 this->N = n_in;
23 this->values = new double[M*N];
24 }/*}}}*/
25 ~Matrix(){/*{{{*/
26 delete [] this->values;
27 }/*}}}*/
28 void Echo(void){/*{{{*/
29 for(int i=0;i<M;i++){
30 for(int j=0;j<N;j++){
31 cout << " " << this->values[i*N+j];
32 }
33 cout << endl;
34 }
35 }/*}}}*/
36 void SetValue(int i,int j,double value){/*{{{*/
37 this->values[i*N+j] = value;
38 }/*}}}*/
39 double GetValue(int i,int j){/*{{{*/
40 return this->values[i*N+j];
41 }/*}}}*/
42 void GetSize(int* pM,int* pN){/*{{{*/
43 *pM = this->M;
44 *pN = this->N;
45 }/*}}}*/
46 double* GetPointer(void){/*{{{*/
47 return this->values;
48 }/*}}}*/
49 void MatrixSum(Matrix* A,Matrix* B){/*{{{*/
50 /*Check that sizes are compatible*/
51 int M_B,N_B,M_A,N_A;
52 B->GetSize(&M_B,&N_B);
53 A->GetSize(&M_A,&N_A);
54 assert(this->M==M_B && this->N==N_B);
55 assert(this->M==M_A && this->N==N_A);
56 for(int i=0;i<this->M;i++){
57 for(int j=0;j<this->N;j++){
58 this->SetValue(i,j,A->GetValue(i,j) + B->GetValue(i,j));
59 }
60 }
61 }/*}}}*/
62 void MatrixDiff(Matrix* A,Matrix* B){/*{{{*/
63 /*Check that sizes are compatible*/
64 int M_B,N_B,M_A,N_A;
65 B->GetSize(&M_B,&N_B);
66 A->GetSize(&M_A,&N_A);
67 assert(this->M==M_B && this->N==N_B);
68 assert(this->M==M_A && this->N==N_A);
69 for(int i=0;i<this->M;i++){
70 for(int j=0;j<this->N;j++){
71 this->SetValue(i,j,A->GetValue(i,j) - B->GetValue(i,j));
72 }
73 }
74 }/*}}}*/
75 void MatrixAbs(Matrix* A){/*{{{*/
76 for(int i=0;i<this->M;i++){
77 for(int j=0;j<this->N;j++){
78 this->SetValue(i,j,fabs(A->GetValue(i,j)));
79 }
80 }
81 }/*}}}*/
82 void MatrixEqual(Matrix* A){/*{{{*/
83 for(int i=0;i<this->M;i++){
84 for(int j=0;j<this->N;j++){
85 this->SetValue(i,j,A->GetValue(i,j));
86 }
87 }
88 }/*}}}*/
89 double MatrixInternSum(){/*{{{*/
90 double sum=0;
91 for(int i=0;i<this->M;i++){
92 for(int j=0;j<this->N;j++){
93 sum+=this->GetValue(i,j);
94 }
95 }
96 return sum;
97 }/*}}}*/
98 void ExtractLine(Matrix* A,int i){/*{{{*/
99 /* Check that the size of A is compatible */
100 int M_A,N_A;
101 A->GetSize(&M_A,&N_A);
102 assert(M_A==1 && this->N==N_A);
103 for(int j=0;j<this->N;j++){
104 A->SetValue(0,j,this->GetValue(i,j));
105 }
106 }/*}}}*/
107 void ExtractColumn(Matrix* A,int j){/*{{{*/
108 /* Check that the size of A is compatible */
109 int M_A,N_A;
110 A->GetSize(&M_A,&N_A);
111 assert(N_A==1 && this->M==M_A);
112 for(int i=0;i<this->M;i++){
113 A->SetValue(i,0,this->GetValue(i,j));
114 }
115 }/*}}}*/
116 void AddNumber(double a){/*{{{*/
117 for(int i=0;i<this->M;i++){
118 for(int j=0;j<this->N;j++){
119 this->SetValue(i,j,this->GetValue(i,j) + a);
120 }
121 }
122 }/*}}}*/
123};/*}}}*/
124
125/*Local prototypes{{{*/
126void makep(Matrix *Pobs,int nx,int ny, int dx, int dy);
127void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny);
128void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel);
129void* plouffT(void* vpthread_handle);
130void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny);
131void reshape(Matrix* V,Matrix* V1,int nx,int ny);
132double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,double dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my);
133void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B);
134double signe(double a);
135void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
136void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd);
137double coolshed(double T0,double k,double c,double D);
138void LaunchThread(void* function(void*), void* usr,int num_threads);
139/*}}}*/
140
141/*Multithreading structures {{{*/
142typedef struct{
143 void* usr;
144 int my_thread;
145 int num_threads;
146} pthread_handle;
147
148typedef struct{
149 Matrix *g;
150 Matrix *Pobs;
151 Matrix *Pp;
152 Matrix *mesh;
153 Matrix *rho;
154 int dx;
155 int dy;
156 int dn;
157 int m;
158 int n;
159 int l;
160} AppStruct;
161/*}}}*/
162
163int main(){/*{{{*/
164
165 /* Seed the random number generator {{{*/
166 srand (time(NULL)); /*}}}*/
167 /* Define the variables {{{*/
168
169 int dx = 500; /* prism dimension in x-direction */
170 int dy = 500; /* prism dimension in y-direction */
171 int mx = 39; /* number of prisms in x-direction */
172 int my = 60; /* number of prisms in y-direction */
173 int nx = 39; /* number of data points in x-direction */
174 int ny = 60; /* number of data points in y-direction */
175 int dn = 15000; /* distance for neighbouting prisms for gravity calculation */
176 double ptval = 100.; /* max. amount to perturb model */
177 double ptval2 = 100.;
178
179 Matrix *Pobs=new Matrix(nx*ny,2); /* data positions */
180 makep(Pobs,nx,ny,dx,dy);
181// Pobs->Echo();
182
183
184 Matrix *Pp=new Matrix(mx*my,2); /* data positions */
185 makep(Pp,mx,my,dx,dy);
186// Pp->Echo();
187
188 double rhoi = 890; /* ice density */
189 double rhow = 1030; /* water density */
190// double rhos = 2013; /* sediment density */
191 double rhoc = 2670; /* bedrock density */
192
193 Matrix *Rho = new Matrix(1,2);
194 Rho->SetValue(0,0,rhoi);
195 Rho->SetValue(0,1,rhow);
196 Matrix *rho1 = new Matrix(1,3);
197 rho1->SetValue(0,0,rhoi);
198 rho1->SetValue(0,1,rhow);
199 rho1->SetValue(0,2,rhoc);
200 Matrix *rho2 = new Matrix(1,2);
201 rho2->SetValue(0,0,rhoi-rhoc);
202 rho2->SetValue(0,1,rhow-rhoc);
203
204 double ctr=1; /* parameter for filtering */
205 double sd=0.1;
206
207 Matrix *xobs= new Matrix(ny,nx);
208 Matrix *yobs= new Matrix(ny,nx);
209
210 vec2grid(Pobs,xobs,yobs,nx,ny);
211// xobs->Echo();
212// yobs->Echo();
213
214 double dlevel= 860; /* mean level of data acquisition */
215
216 double mmax = 1000; /* max value for layer interfaces */
217 double mmax2 = 1000;
218 double mmax3 = 1000;
219
220 /* control parameter for temperature schedule */
221
222 double ca=0.9; /* for acceptance */
223 double cm=0.5; /* for model perturbation */
224
225 double T0a = 0.1; /* initial temperature for acceptance */
226 double T0m = 0.9; /* initial temperature for model perturbation */
227 double D = 2; /* dimension of the model */
228 int maxconsecrej = 1000; /* max consecutive rejection */
229 int maxsuccess = 100; /* max number of success within one temperature */
230 double T_min = 1e-10; /* stopping temp */
231 double Tred = 1;
232 double E_min = -1000000;
233 double E_exp = 0.0291; /* expected misfit */
234 int maxiter = 10000;
235 int maxtotaliter = 1000000;
236 double Tol = 1e-10; /* tolerance on misfit */
237 int sfreq = 100;
238
239 /*}}}*/
240 /* load the data {{{*/
241
242 /* Observed gravity anomaly */
243
244 ifstream file1("store_fa500_36s.txt");
245 Matrix * gobs= new Matrix(nx*ny,1);
246 double inputnumber;
247 for(int i=0;i<ny*nx; i++){
248 file1 >> inputnumber;
249 gobs->SetValue(i,0, inputnumber*1e-5);
250 }
251 file1.close();
252// gobs->Echo();
253
254 /* load data about the ice thickness */
255
256 ifstream file2("store_flag_icethick500.txt");
257 Matrix * icethick= new Matrix(mx*my,1);
258 for(int s=0;s<mx*my; s++){
259 file2 >> inputnumber;
260 icethick->SetValue(s,0,inputnumber);
261 }
262 file2.close();
263// icethick->Echo();
264
265 /* load the batimethry data */
266
267 ifstream file3("store_flag_bathy500.txt");
268 Matrix * bathy= new Matrix(mx*my,1);
269 for(int s=0;s<mx*my; s++){
270 file3 >> inputnumber;
271 bathy->SetValue(s,0,inputnumber);
272 }
273 file3.close();
274// bathy->Echo();
275
276 /* id of grid to evaluate misfit */
277
278
279 ifstream file4("store_flag_eval500.txt");
280 Matrix * evalid= new Matrix(mx*my,1);
281 for(int s=0;s<mx*my; s++){
282 file4 >> inputnumber;
283 evalid->SetValue(s,0,inputnumber);
284 }
285 file4.close();
286// evalid->Echo();
287
288 /* initial guess of the model */
289
290 ifstream file5("m0_140114b.txt");
291 Matrix * mesh_ini= new Matrix(mx*my,3);
292 for(int s=0;s<mx*my; s++){
293 for(int j=0;j<3;j++){
294 file5 >> inputnumber;
295 mesh_ini->SetValue(s,j,inputnumber);
296 }
297 }
298 file5.close();
299// mesh_ini->Echo();
300 /*}}}*/
301 /* VFSA {{{ */
302
303 /* name of the files to save results */
304 std::ofstream savefile1 ("r_140114b.txt");
305 std::ofstream savefile2("m_140114b.txt");
306
307 /* counters initialization */
308 int success = 0;
309 int finished = 0;
310 int consec = 0;
311 double Ta = T0a;
312 double Tm = T0m;
313 int iterT = 0; /* iteration within a T */
314 int total = 0; /* total number of iteration */
315 int totaliter = 0;
316 int msave = 0;
317 double E_new;
318 double E_final;
319 double dE;
320 double P;
321 double rn;
322 Matrix* m_old = new Matrix(mx *my,3);
323 Matrix* m_min = new Matrix(mx *my,3);
324 Matrix* m_new = new Matrix(mx *my,3);
325 m_old->MatrixEqual(mesh_ini);
326
327 /* calculate initial misfit */
328 double E_old=misfit(m_old,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my);
329
330 /* record initial settings */
331 savefile1 << "P "<< "Ta "<< "Tm "<< "Eold "<< "totaliter "<< "Tred "<< endl;
332 savefile1 << "nan "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
333 savefile2 << totaliter<< endl;
334 for(int i=0;i<mx*my;i++){
335 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
336 }
337 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
338
339 /* beginning of the loop */
340
341 while(finished==0){
342
343 iterT++;
344 totaliter++;
345
346 /* stop or reduce T */
347 if(iterT>=maxiter || success>maxsuccess){
348 if(Ta<T_min || total>maxtotaliter || fabs(E_old)<=Tol){
349 finished=1;
350 total+=iterT;
351 break;
352 }
353 else{ /* reduce T */
354 Ta=coolshed(T0a,Tred,ca,D);
355 Tm=coolshed(T0m,Tred,cm,D);
356 total+=iterT;
357 iterT=0;
358 success=1;
359 Tred++;
360 consec=0;
361 }
362 }
363
364 /* update model and calculate energy */
365
366 newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd); /* new model */
367 E_new=misfit(m_new,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my); /* new energy */
368 dE=E_new-E_old; /* energy difference */
369
370 /* acceptance probability */
371
372 P=exp(-dE/Ta);
373
374 /* stop if energy is lower than specified minimum */
375 if (E_new<E_min){
376 m_old->MatrixEqual(m_new);
377 E_old=E_new;
378 break;
379 }
380
381 rn=rand()/double (RAND_MAX);
382
383 /* accept new model or not */
384 if(dE<=0){
385 m_old->MatrixEqual(m_new);
386 E_old=E_new;
387 success++;
388 consec=0;
389 savefile1 << P<<" "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
390 if(Ta<1e-3){
391 savefile2 << totaliter<< endl;
392 for(int i=0;i<mx*my;i++){
393 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
394 }
395 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
396 }
397 }
398 else{
399 if(P>rn){
400 m_old->MatrixEqual(m_new);
401 E_old=E_new;
402 success++;
403 consec=0;
404 savefile1 << P<<" "<< Ta<<" "<< Tm<<" "<< E_old<<" "<< totaliter<<" "<< Tred <<" "<< endl;
405 if(Ta<1e-3){
406 savefile2 << totaliter<< endl;
407 for(int i=0;i<mx*my;i++){
408 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
409 }
410 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
411 }
412 }
413 else{
414 consec++;
415 }
416 }
417 cout<<totaliter<<endl;
418 }
419
420 m_min->MatrixEqual(m_old);
421 E_final=E_old;
422 savefile1 << "nan"<<" "<< "nan"<<" "<< "nan"<<" "<< E_final<<" "<< "nan"<<" "<< "nan" <<" "<< endl;
423 savefile2 << " Mesh final"<< endl;
424 for(int i=0;i<mx*my;i++){
425 savefile2 << m_min->GetValue(i,0)<<" "<< m_min->GetValue(i,1)<<" "<< m_min->GetValue(i,2)<<endl;
426 }
427 savefile1.close();
428 savefile2.close();
429
430 delete m_old;
431 delete m_min;
432 delete m_new;
433 delete Pobs;
434 delete Pp;
435 delete Rho;
436 delete rho1;
437 delete rho2;
438 delete xobs;
439 delete yobs;
440 delete mesh_ini;
441 delete bathy;
442 delete icethick;
443 delete evalid;
444
445 /*}}}*/
446 return 0;
447}/*}}}*/
448
449void GSLsquarefit(Matrix** pX,Matrix* A,Matrix* B){/*{{{*/
450
451 /*GSL Matrices and vectors: */
452 int M,N;
453 double chisq;
454
455 /*Get system size*/
456 A->GetSize(&M,&N);
457
458 /*Initialize gsl matrices and vectors: */
459 gsl_matrix* a = gsl_matrix_alloc(M,N);
460 for(int i=0;i<M;i++){
461 for(int j=0;j<N;j++){
462 gsl_matrix_set (a,i,j,A->GetValue(i,j));
463 }
464 }
465 gsl_vector* b = gsl_vector_alloc(M);
466 for(int i=0;i<M;i++){
467 gsl_vector_set(b,i,B->GetValue(i,0));
468 }
469
470 gsl_vector* x = gsl_vector_alloc(N);
471 gsl_matrix* cov = gsl_matrix_alloc(N,N);
472
473 /*Least square fit: */
474 gsl_multifit_linear_workspace* work = gsl_multifit_linear_alloc(M,N);
475 gsl_multifit_linear (a, b, x, cov, &chisq, work);
476 gsl_multifit_linear_free (work);
477
478 /*Clean up and assign output pointer*/
479 Matrix* X = new Matrix(N,1);
480 for(int j=0;j<N;j++){
481 X->SetValue(j,0,gsl_vector_get(x,j));
482 }
483 *pX = X;
484
485 gsl_matrix_free(a);
486 gsl_vector_free(b);
487 gsl_vector_free(x);
488 gsl_matrix_free(cov);
489
490}/*}}}*/
491void makep(Matrix *Pobs,int nx,int ny, int dx, int dy){/*{{{*/
492 for(int i=0;i<ny;i++){
493 for(int j=0;j<nx;j++){
494 Pobs->SetValue(j+nx*i,0,j*dx);
495 Pobs->SetValue(j+nx*i,1,i*dy);
496 }
497 }
498}/*}}}*/
499void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny){/*{{{*/
500 for(int i=0;i<ny;i++){
501 for (int j=0;j<nx;j++){
502 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
503 V2->SetValue(i,j, V->GetValue(j+nx*i,1));
504 }
505 }
506}/*}}}*/
507void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel){/*{{{*/
508 int sizem1,sizem2;
509 m->GetSize(&sizem1,&sizem2);
510 for(int i=0;i<sizem1;i++){
511 for(int j=0;j<sizem2+1;j++){
512 if(j<sizem2){
513 m1->SetValue(i,j,1e-10*(sizem2+1-j));
514 m2->SetValue(i,j,m->GetValue(i,j));
515 if(m->GetValue(i,j)<0){
516 m1->SetValue(i,j,m->GetValue(i,j));
517 m2->SetValue(i,j,i*1e-10);
518 }
519 }
520 else{
521 m1->SetValue(i,j,1e-10);
522 }
523 }
524 }
525 m1->AddNumber(dlevel);
526 m2->AddNumber(dlevel);
527}/*}}}*/
528void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny){/*{{{*/
529 for(int i=0;i<ny;i++){
530 for (int j=0;j<nx;j++){
531 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
532 }
533 }
534}/*}}}*/
535void reshape(Matrix* V,Matrix* V1,int nx,int ny){/*{{{*/
536 for (int i=0;i<ny;i++){
537 for(int j=0;j<nx;j++){
538 V1->SetValue(j+nx*i,0,V->GetValue(i,j));
539 }
540 }
541}/*}}}*/
542double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,double dlevel,Matrix* Pobs,Matrix* xobs,Matrix* yobs,Matrix* Pp,Matrix* rho1, Matrix* rho2,int dx,int dy,int dn,int nx,int ny, int mx,int my){/*{{{*/
543 Matrix* m1=new Matrix(mx*my,4);
544 Matrix* m2=new Matrix(mx*my,3);
545 Matrix* g1=new Matrix(nx*ny,1);
546 Matrix* g2=new Matrix(nx*ny,1);
547 Matrix* g=new Matrix(nx*ny,1);
548 Matrix* gcalgr=new Matrix(ny,nx);
549 Matrix* gcalvec=new Matrix(mx*my,1);
550 Matrix* df=new Matrix(mx*my,1);
551 Matrix* G=new Matrix(mx*my,3);
552 double a=0;
553 double b=0;
554 double e=0;
555 msplit(m0,m1,m2,dlevel);
556
557 /*Multithreaded core*/
558 int num_threads = 8;
559 AppStruct usr;
560 usr.g = g1;
561 usr.Pobs = Pobs;
562 usr.Pp=Pp;
563 usr.mesh= m1;
564 usr.rho= rho1;
565 usr.dx=dx;
566 usr.dy=dy;
567 usr.dn=dn;
568 usr.m=mx*my;
569 usr.n=nx*ny;
570 usr.l=4;
571 LaunchThread(plouffT,(void*)&usr,num_threads);
572 usr.g = g2;
573 usr.mesh= m2;
574 usr.rho= rho2;
575 usr.l=3;
576 LaunchThread(plouffT,(void*)&usr,num_threads);
577 g->MatrixSum(g1,g2);
578 vec2gridsimple(g,gcalgr,nx,ny);
579 reshape(gcalgr,gcalvec,mx,my);
580 for (int i=0;i<mx*my;i++){
581 df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
582 G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
583 G->SetValue(i,1,evalid->GetValue(i,0)*Pobs->GetValue(i,1));
584 G->SetValue(i,2,evalid->GetValue(i,0));
585 }
586 Matrix* M = NULL;
587 GSLsquarefit(&M,G,df);
588
589 for (int i=0;i<my;i++){
590 for(int j=0;j<mx;j++){
591 gcalgr->SetValue(i,j,gcalgr->GetValue(i,j)+xobs->GetValue(i,j)*M->GetValue(0,0)+yobs->GetValue(i,j)*M->GetValue(1,0)+M->GetValue(2,0));
592 }
593 }
594 reshape(gcalgr,g,mx,my);
595 for (int i=0;i<mx*my;i++){
596 a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
597 b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
598 }
599 e=2*a/(a+b);
600
601 delete m1;
602 delete m2;
603 delete g1;
604 delete g2;
605 delete g;
606 delete gcalgr;
607 delete gcalvec;
608 delete df;
609 delete G;
610 delete M;
611
612 return e;
613}/*}}}*/
614void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd){/*{{{*/
615 Matrix* m1gr=new Matrix(my,mx);
616 Matrix* m1grsm=new Matrix(my,mx);
617 Matrix* m1col=new Matrix(mx*my,1);
618 double u=0;
619 double y=0;
620 m1->MatrixEqual(m0);
621 for (int i=0;i<mx*my;i++){
622 if(icethick->GetValue(i,0)==0){
623 u=double (rand())/ double(RAND_MAX);
624 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
625 m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
626 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
627 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
628 }
629 if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
630 m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
631 }
632 }
633 }
634 m1->ExtractColumn(m1col,1);
635 vec2gridsimple(m1col,m1gr,mx,my);
636 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
637 reshape(m1grsm,m1col,mx,my);
638 for (int i=0;i<mx*my;i++){
639 if(icethick->GetValue(i,0)==0){
640 m1->SetValue(i,1,m1col->GetValue(i,0));
641 }
642 else{
643 m1->SetValue(i,1,m0->GetValue(i,1));
644 }
645 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
646 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
647 }
648 }
649
650 for (int i=0;i<mx*my;i++){
651 if(bathy->GetValue(i,0)==0){
652 u=double (rand())/ double(RAND_MAX);
653 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
654 m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
655 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
656 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
657 }
658 if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
659 m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
660 }
661 }
662 }
663 m1->ExtractColumn(m1col,2);
664 vec2gridsimple(m1col,m1gr,mx,my);
665 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
666 reshape(m1grsm,m1col,mx,my);
667 for (int i=0;i<mx*my;i++){
668 if(bathy->GetValue(i,0)==0){
669 m1->SetValue(i,2,m1col->GetValue(i,0));
670 }
671 else{
672 m1->SetValue(i,2,m0->GetValue(i,2));
673 }
674 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
675 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
676 }
677 }
678 delete m1gr;
679 delete m1grsm;
680 delete m1col;
681}/*}}}*/
682double signe(double a){/*{{{*/
683 if(a<0){return -1;}
684 else{return 1;}
685}/*}}}*/
686void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my){/*{{{*/
687 A->MatrixEqual(Ain);
688 for (int i=1;i<my-1;i++){
689 for(int j=1;j<mx-1;j++){
690 A->SetValue(i,j,(ctr*Ain->GetValue(i,j)+sd*(Ain->GetValue(i-1,j)+Ain->GetValue(i+1,j)+Ain->GetValue(i,j-1)+Ain->GetValue(i,j+1)))/(ctr+4*sd));
691 }
692 }
693}/*}}}*/
694double coolshed(double T0,double k,double c,double D){/*{{{*/
695 double T1=T0*exp(-c*pow(k,1/D));
696 return T1;
697}/*}}}*/
698void* plouffT(void* vpthread_handle){/*{{{*/
699
700 /*recover this thread info*/
701 pthread_handle *handle = (pthread_handle*)vpthread_handle;
702 int my_thread = handle->my_thread;
703 int num_threads = handle->num_threads;
704
705 /*Recover struct*/
706 AppStruct *usr = (AppStruct*)handle->usr;
707 Matrix *g = usr->g;
708 Matrix *Pobs = usr->Pobs;
709 Matrix *Pp = usr->Pp;
710 Matrix *mesh = usr->mesh;
711 Matrix *rho = usr->rho;
712 int dx =usr->dx;
713 int dy =usr->dy;
714 int dn =usr->dn;
715 int m =usr->m;
716 int n =usr->n;
717 int l =usr->l;
718
719 double gg=6.673e-11;
720 int si,sj,id,s;
721 double R,Q,P;
722 Matrix *xp= new Matrix(1,2);
723 Matrix *yp= new Matrix(1,2);
724 Matrix *xpp= new Matrix(1,2);
725 Matrix *ypp= new Matrix(1,2);
726 Matrix *U= new Matrix(l,4);
727 Matrix *U1=new Matrix(1,4);
728 Matrix *U2=new Matrix(1,4);
729 Matrix *gl= new Matrix(1,l-1);
730 bool test=true;
731 for(int c=my_thread;c<n;c+=num_threads){
732 g->SetValue(c,0,0);
733 for(int a=0;a<m;a++){
734 test=true;
735 xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
736 xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
737 if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
738 xpp->SetValue(0,0,xp->GetValue(0,1));
739 xpp->SetValue(0,1,xp->GetValue(0,0));
740 xp->MatrixAbs(xpp);
741 }
742 yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
743 yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
744 if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
745 ypp->SetValue(0,0,yp->GetValue(0,1));
746 ypp->SetValue(0,1,yp->GetValue(0,0));
747 yp->MatrixAbs(ypp);
748 }
749 P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
750 if(P>dn){
751 test=false;
752 for(int i=0;i<l-1;i++){
753 gl->SetValue(0,i,0);
754 }
755 }
756 if(test==true){
757 si=1;
758 sj=1;
759 id=0;
760 for(int i=0;i<2;i++){
761 si*=-1;
762 for(int j=0;j<2;j++){
763 si*=-1;
764 s=si*sj;
765 for(int k=0;k<l;k++){
766 R=sqrt(xp->GetValue(0,i)*xp->GetValue(0,i)+yp->GetValue(0,j)*yp->GetValue(0,j)+mesh->GetValue(a,k)*mesh->GetValue(a,k));
767 Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
768 U->SetValue(k,id,s*(mesh->GetValue(a,k)*Q-xp->GetValue(0,i)*log(R+yp->GetValue(0,j))-yp->GetValue(0,j)*log(R+xp->GetValue(0,i))));
769 }
770 id++;
771 }
772 }
773 for(int b=0;b<l-1;b++){
774 U->ExtractLine(U1,b);
775 U->ExtractLine(U2,b+1);
776 gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
777 }
778 }
779 g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
780 }
781 }
782 delete xp;
783 delete yp;
784 delete xpp;
785 delete ypp;
786 delete gl;
787 delete U;
788 delete U1;
789 delete U2;
790
791 return NULL;
792}/*}}}*/
793void LaunchThread(void* function(void*), void* usr,int num_threads){/*{{{*/
794
795 int i;
796 int *status = NULL;
797 pthread_t *threads = NULL;
798 pthread_handle *handles = NULL;
799
800 /*dynamically allocate: */
801 threads=(pthread_t*)malloc(num_threads*sizeof(pthread_t));
802 handles=(pthread_handle*)malloc(num_threads*sizeof(pthread_handle));
803
804 for(i=0;i<num_threads;i++){
805 handles[i].usr=usr;
806 handles[i].my_thread =i;
807 handles[i].num_threads=num_threads;
808 }
809 for(i=0;i<num_threads;i++){
810 if(pthread_create(threads+i,NULL,function,(void*)(handles+i))){
811 std::cerr<<"pthread_create error";
812 }
813 }
814 for(i=0;i<num_threads;i++){
815 if(pthread_join(threads[i],(void**)&status)){
816 std::cerr<<"pthread_join error";
817 }
818 }
819
820 /*Free ressources:*/
821 delete threads;
822 delete handles;
823}/*}}}*/
Note: See TracBrowser for help on using the repository browser.