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

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

MODIF: details
File size: 22.2 KB
RevLine 
[18434]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>
[18457]10#include <pthread.h>
[18434]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);
[18457]129void* plouffT(void* vpthread_handle);
[18434]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,int 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);
[18435]134double signe(double a);
135void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
[18439]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);
[18446]137double coolshed(double T0,int k,double c,double D);
[18457]138void LaunchThread(void* function(void*), void* usr,int num_threads);
[18434]139/*}}}*/
140
[18457]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
[18434]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 */
[18446]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
[18434]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 /*}}}*/
[18446]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++){
[18454]335 savefile2 << m_old->GetValue(i,0)<<" "<< m_old->GetValue(i,1)<<" "<< m_old->GetValue(i,2)<<endl;
[18446]336 }
[18454]337 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
338
[18446]339 /* beginning of the loop */
340
341 while(finished==0){
[18454]342
[18446]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 }
[18434]362 }
[18446]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 */
[18454]367 E_new=misfit(m_new,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my); /* new energy */
[18446]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 }
[18454]395 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
[18446]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 }
[18454]410 savefile2 << "111111111111111111111111111111111111111111111111111111111111111111111111111"<< endl;
[18446]411 }
412 }
413 else{
414 consec++;
415 }
416 }
[18457]417 cout<<totaliter<<endl;
[18434]418 }
419
[18446]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 }
[18457]427 savefile1.close();
428 savefile2.close();
[18461]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
[18446]445 /*}}}*/
[18434]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_matrix_free(cov);
488
489}/*}}}*/
490void makep(Matrix *Pobs,int nx,int ny, int dx, int dy){/*{{{*/
491 for(int i=0;i<ny;i++){
492 for(int j=0;j<nx;j++){
493 Pobs->SetValue(j+nx*i,0,j*dx);
494 Pobs->SetValue(j+nx*i,1,i*dy);
495 }
496 }
497}/*}}}*/
498void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny){/*{{{*/
499 for(int i=0;i<ny;i++){
500 for (int j=0;j<nx;j++){
501 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
502 V2->SetValue(i,j, V->GetValue(j+nx*i,1));
503 }
504 }
505}/*}}}*/
506void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel){/*{{{*/
507 int sizem1,sizem2;
508 m->GetSize(&sizem1,&sizem2);
509 for(int i=0;i<sizem1;i++){
510 for(int j=0;j<sizem2+1;j++){
511 if(j<sizem2){
512 m1->SetValue(i,j,1e-10*(sizem2+1-j));
513 m2->SetValue(i,j,m->GetValue(i,j));
514 if(m->GetValue(i,j)<0){
515 m1->SetValue(i,j,m->GetValue(i,j));
516 m2->SetValue(i,j,i*1e-10);
517 }
518 }
519 else{
520 m1->SetValue(i,j,1e-10);
521 }
522 }
523 }
524 m1->AddNumber(dlevel);
525 m2->AddNumber(dlevel);
526}/*}}}*/
527void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny){/*{{{*/
528 for(int i=0;i<ny;i++){
529 for (int j=0;j<nx;j++){
530 V1->SetValue(i,j, V->GetValue(j+nx*i,0));
531 }
532 }
533}/*}}}*/
534void reshape(Matrix* V,Matrix* V1,int nx,int ny){/*{{{*/
535 for (int i=0;i<ny;i++){
536 for(int j=0;j<nx;j++){
537 V1->SetValue(j+nx*i,0,V->GetValue(i,j));
538 }
539 }
540}/*}}}*/
541double misfit(Matrix* m0,Matrix* evalid,Matrix* gobs,int 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){/*{{{*/
542 Matrix* m1=new Matrix(mx*my,4);
543 Matrix* m2=new Matrix(mx*my,3);
544 Matrix* g1=new Matrix(nx*ny,1);
545 Matrix* g2=new Matrix(nx*ny,1);
546 Matrix* g=new Matrix(nx*ny,1);
547 Matrix* gcalgr=new Matrix(ny,nx);
548 Matrix* gcalvec=new Matrix(mx*my,1);
549 Matrix* df=new Matrix(mx*my,1);
550 Matrix* G=new Matrix(mx*my,3);
551 double a=0;
552 double b=0;
553 double e=0;
554 msplit(m0,m1,m2,dlevel);
[18457]555
556 /*Multithreaded core*/
[18461]557 int num_threads = 8;
[18457]558 AppStruct usr;
559 usr.g = g1;
560 usr.Pobs = Pobs;
561 usr.Pp=Pp;
562 usr.mesh= m1;
563 usr.rho= rho1;
564 usr.dx=dx;
565 usr.dy=dy;
566 usr.dn=dn;
567 usr.m=mx*my;
568 usr.n=nx*ny;
569 usr.l=4;
570 LaunchThread(plouffT,(void*)&usr,num_threads);
571 usr.g = g2;
572 usr.mesh= m2;
573 usr.rho= rho2;
574 usr.l=3;
575 LaunchThread(plouffT,(void*)&usr,num_threads);
[18434]576 g->MatrixSum(g1,g2);
577 vec2gridsimple(g,gcalgr,nx,ny);
578 reshape(gcalgr,gcalvec,mx,my);
579 for (int i=0;i<mx*my;i++){
580 df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
581 G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
582 G->SetValue(i,1,evalid->GetValue(i,0)*Pobs->GetValue(i,1));
583 G->SetValue(i,2,evalid->GetValue(i,0));
584 }
585 Matrix* M = NULL;
586 GSLsquarefit(&M,G,df);
587
588 for (int i=0;i<my;i++){
589 for(int j=0;j<mx;j++){
590 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));
591 }
592 }
593 reshape(gcalgr,g,mx,my);
594 for (int i=0;i<mx*my;i++){
595 a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
596 b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
597 }
598 e=2*a/(a+b);
[18461]599
600 delete m1;
601 delete m2;
602 delete g1;
603 delete g2;
604 delete g;
605 delete gcalgr;
606 delete gcalvec;
607 delete df;
608 delete G;
609 delete M;
610
[18434]611 return e;
612}/*}}}*/
[18435]613void 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){/*{{{*/
614 Matrix* m1gr=new Matrix(my,mx);
615 Matrix* m1grsm=new Matrix(my,mx);
616 Matrix* m1col=new Matrix(mx*my,1);
617 double u=0;
618 double y=0;
619 m1->MatrixEqual(m0);
620 for (int i=0;i<mx*my;i++){
621 if(icethick->GetValue(i,0)==0){
622 u=double (rand())/ double(RAND_MAX);
623 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
624 m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
625 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
626 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
627 }
628 if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
629 m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
630 }
631 }
632 }
633 m1->ExtractColumn(m1col,1);
634 vec2gridsimple(m1col,m1gr,mx,my);
635 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
636 reshape(m1grsm,m1col,mx,my);
637 for (int i=0;i<mx*my;i++){
638 if(icethick->GetValue(i,0)==0){
639 m1->SetValue(i,1,m1col->GetValue(i,0));
640 }
641 else{
642 m1->SetValue(i,1,m0->GetValue(i,1));
643 }
644 if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
645 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
646 }
647 }
648
649 for (int i=0;i<mx*my;i++){
650 if(bathy->GetValue(i,0)==0){
651 u=double (rand())/ double(RAND_MAX);
652 y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
653 m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
654 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
655 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
656 }
657 if(m1->GetValue(i,2)>=m1->GetValue(i,1)+mmax2){
658 m1->SetValue(i,2,m1->GetValue(i,1)+mmax2);
659 }
660 }
661 }
662 m1->ExtractColumn(m1col,2);
663 vec2gridsimple(m1col,m1gr,mx,my);
664 filtergrav(m1grsm,m1gr,ctr,sd,mx,my);
665 reshape(m1grsm,m1col,mx,my);
666 for (int i=0;i<mx*my;i++){
667 if(bathy->GetValue(i,0)==0){
668 m1->SetValue(i,2,m1col->GetValue(i,0));
669 }
670 else{
671 m1->SetValue(i,2,m0->GetValue(i,2));
672 }
673 if(m1->GetValue(i,2)<=m1->GetValue(i,1)){
674 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
675 }
676 }
677 delete m1gr;
678 delete m1grsm;
679 delete m1col;
680}/*}}}*/
681double signe(double a){/*{{{*/
682 if(a<0){return -1;}
683 else{return 1;}
684}/*}}}*/
685void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my){/*{{{*/
[18439]686 A->MatrixEqual(Ain);
[18435]687 for (int i=1;i<my-1;i++){
688 for(int j=1;j<mx-1;j++){
689 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));
690 }
691 }
692}/*}}}*/
[18446]693double coolshed(double T0,int k,double c,double D){/*{{{*/
694 double T1=T0*exp(-c*pow(k,1/D));
695 return T1;
696}/*}}}*/
[18457]697void* plouffT(void* vpthread_handle){/*{{{*/
698
699 /*recover this thread info*/
700 pthread_handle *handle = (pthread_handle*)vpthread_handle;
701 int my_thread = handle->my_thread;
702 int num_threads = handle->num_threads;
703
704 /*Recover struct*/
705 AppStruct *usr = (AppStruct*)handle->usr;
706 Matrix *g = usr->g;
707 Matrix *Pobs = usr->Pobs;
708 Matrix *Pp = usr->Pp;
709 Matrix *mesh = usr->mesh;
710 Matrix *rho = usr->rho;
711 int dx =usr->dx;
712 int dy =usr->dy;
713 int dn =usr->dn;
714 int m =usr->m;
715 int n =usr->n;
716 int l =usr->l;
717
718 double gg=6.673e-11;
719 int si,sj,id,s;
720 double R,Q,P;
721 Matrix *xp= new Matrix(1,2);
722 Matrix *yp= new Matrix(1,2);
723 Matrix *xpp= new Matrix(1,2);
724 Matrix *ypp= new Matrix(1,2);
725 Matrix *U= new Matrix(l,4);
726 Matrix *U1=new Matrix(1,4);
727 Matrix *U2=new Matrix(1,4);
728 Matrix *gl= new Matrix(1,l-1);
729 bool test=true;
730 for(int c=my_thread;c<n;c+=num_threads){
731 g->SetValue(c,0,0);
732 for(int a=0;a<m;a++){
733 test=true;
734 xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
735 xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
736 if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
737 xpp->SetValue(0,0,xp->GetValue(0,1));
738 xpp->SetValue(0,1,xp->GetValue(0,0));
739 xp->MatrixAbs(xpp);
740 }
741 yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
742 yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
743 if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
744 ypp->SetValue(0,0,yp->GetValue(0,1));
745 ypp->SetValue(0,1,yp->GetValue(0,0));
746 yp->MatrixAbs(ypp);
747 }
748 P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
749 if(P>dn){
750 test=false;
751 for(int i=0;i<l-1;i++){
752 gl->SetValue(0,i,0);
753 }
754 }
755 if(test==true){
756 si=1;
757 sj=1;
758 id=0;
759 for(int i=0;i<2;i++){
760 si*=-1;
761 for(int j=0;j<2;j++){
762 si*=-1;
763 s=si*sj;
764 for(int k=0;k<l;k++){
765 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));
766 Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
767 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))));
768 }
769 id++;
770 }
771 }
772 for(int b=0;b<l-1;b++){
773 U->ExtractLine(U1,b);
774 U->ExtractLine(U2,b+1);
775 gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
776 }
777 }
778 g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
779 }
780 }
781 delete xp;
782 delete yp;
783 delete xpp;
784 delete ypp;
785 delete gl;
786 delete U;
787 delete U1;
788 delete U2;
789
790 return NULL;
791}/*}}}*/
792void LaunchThread(void* function(void*), void* usr,int num_threads){/*{{{*/
793
794 int i;
795 int *status = NULL;
796 pthread_t *threads = NULL;
797 pthread_handle *handles = NULL;
798
799 /*dynamically allocate: */
800 threads=(pthread_t*)malloc(num_threads*sizeof(pthread_t));
801 handles=(pthread_handle*)malloc(num_threads*sizeof(pthread_handle));
802
803 for(i=0;i<num_threads;i++){
804 handles[i].usr=usr;
805 handles[i].my_thread =i;
806 handles[i].num_threads=num_threads;
807 }
808 for(i=0;i<num_threads;i++){
809 if(pthread_create(threads+i,NULL,function,(void*)(handles+i))){
810 std::cerr<<"pthread_create error";
811 }
812 }
813 for(i=0;i<num_threads;i++){
814 if(pthread_join(threads[i],(void**)&status)){
815 std::cerr<<"pthread_join error";
816 }
817 }
818
819 /*Free ressources:*/
820 delete threads;
821 delete handles;
822}/*}}}*/
Note: See TracBrowser for help on using the repository browser.