Changeset 18550


Ignore:
Timestamp:
09/26/14 15:31:42 (11 years ago)
Author:
abuzzi
Message:

CNG: adds landmask

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/contrib/gravity/vfsa_mpi.cpp

    r18491 r18550  
    133133double signe(double a);
    134134void filtergrav(Matrix* A,Matrix* Ain,double ctr,double sd,int mx,int my);
    135 void 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);
     135void 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, Matrix *landmask);
    136136double coolshed(double T0,double k,double c,double D);
    137137/*}}}*/
    … …  
    165165
    166166
    167         Matrix *Pp=new Matrix(mx*my,2); /* data positions */
     167        Matrix *Pp=new Matrix(mx*my,2); /* prisms positions */
    168168        makep(Pp,mx,my,dx,dy);
    169169        // Pp->Echo();
    170170
    171         double  rhoi = 917;           /* ice density     */
     171        double  rhoi = 890;           /* ice density     */
    172172        double  rhow = 1030;          /* water density   */
    173173        // double  rhos = 2013;               /* sediment density */
    … …  
    185185        rho2->SetValue(0,1,rhow-rhoc);
    186186
     187        double dlevel=860;         /* level of data acquisition */
     188
    187189        double ctr=1;            /* parameter for filtering */
    188190        double sd=0.5;
    … …  
    195197        //      yobs->Echo();
    196198
    197         double dlevel= 860;                /* mean level of data acquisition */
    198199
    199200        double mmax  = 1000;               /* max value for layer interfaces */
    … …  
    223224        /* load the data {{{*/
    224225
     226        /*landmask */
     227
     228        ifstream file("landmask.txt");
     229        Matrix * landmask= new Matrix(nx*ny,1);
     230        double inputnumber;
     231        for(int i=0;i<ny*nx; i++){
     232                file >> inputnumber;
     233                landmask->SetValue(i,0,inputnumber);
     234        }
     235        file.close();
     236
    225237        /* Observed gravity anomaly */
    226238
    227239        ifstream file1("store_fa500_36s.txt");
    228240        Matrix * gobs= new Matrix(nx*ny,1);
    229         double inputnumber;
    230241        for(int i=0;i<ny*nx; i++){
    231242                file1 >> inputnumber;
    … …  
    261272
    262273        ifstream file4("store_flag_eval500.txt");
    263         Matrix * evalid= new Matrix(mx*my,1);
    264         for(int s=0;s<mx*my; s++){
     274        Matrix * evalid= new Matrix(nx*ny,1);
     275        for(int s=0;s<nx*ny; s++){
    265276                file4 >> inputnumber;
    266277                evalid->SetValue(s,0,inputnumber);
    … …  
    346357                /* update model and calculate energy */
    347358
    348                 newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd);  /* new model */
     359                newmodelgen(m_old,m_new,bathy,icethick,mx,my,Tm,ptval,mmax,mmax2,ctr,sd, landmask);  /* new model */
    349360                E_new=misfit(m_new,evalid,gobs,dlevel,Pobs,xobs,yobs,Pp,rho1,rho2,dx,dy,dn,nx,ny,mx,my, my_rank, num_procs); /* new energy */
    350361                dE=E_new-E_old;                                        /* energy difference */
    … …  
    508519                                        m2->SetValue(i,j,i*1e-10);
    509520                                }
     521                                m1->SetValue(i,j,m1->GetValue(i,j));
     522                                m2->SetValue(i,j,m2->GetValue(i,j));
    510523                        }
    511524                        else{
    … …  
    619632        Matrix* g=new Matrix(nx*ny,1);
    620633        Matrix* gcalgr=new Matrix(ny,nx);
    621         Matrix* gcalvec=new Matrix(mx*my,1);
    622         Matrix* df=new Matrix(mx*my,1);
    623         Matrix* G=new Matrix(mx*my,3);
     634        Matrix* gcalvec=new Matrix(nx*ny,1);
     635        Matrix* df=new Matrix(nx*ny,1);
     636        Matrix* G=new Matrix(nx*ny,3);
    624637        double a=0;
    625638        double b=0;
    … …  
    630643        g->MatrixSum(g1,g2);
    631644        vec2gridsimple(g,gcalgr,nx,ny);
    632         reshape(gcalgr,gcalvec,mx,my);
    633         for (int i=0;i<mx*my;i++){
     645        reshape(gcalgr,gcalvec,nx,ny);
     646        for (int i=0;i<nx*ny;i++){
    634647                df->SetValue(i,0,evalid->GetValue(i,0)*(gobs->GetValue(i,0)-gcalvec->GetValue(i,0)));
    635648                G->SetValue(i,0,evalid->GetValue(i,0)*Pobs->GetValue(i,0));
    … …  
    640653        GSLsquarefit(&M,G,df);
    641654
    642         for (int i=0;i<my;i++){
    643                 for(int j=0;j<mx;j++){
     655        for (int i=0;i<ny;i++){
     656                for(int j=0;j<nx;j++){
    644657                        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));
    645658                }
    646659        }
    647         reshape(gcalgr,g,mx,my);
    648         for (int i=0;i<mx*my;i++){
     660        reshape(gcalgr,g,nx,ny);
     661        for (int i=0;i<nx*ny;i++){
    649662                a=a+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)-g->GetValue(i,0)));
    650663                b=b+fabs(evalid->GetValue(i,0)*(gobs->GetValue(i,0)+g->GetValue(i,0)));
    … …  
    665678        return e;
    666679}/*}}}*/
    667 void 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){/*{{{*/
     680void 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, Matrix *landmask){/*{{{*/
    668681        Matrix* m1gr=new Matrix(my,mx);
    669682        Matrix* m1grsm=new Matrix(my,mx);
    … …  
    679692        /* first layer: ice */
    680693        for (int i=0;i<mx*my;i++){
    681                 if(nptflag->GetValue(i,0)==0){
    682                         u=double(rand())/double(RAND_MAX);
    683                         y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
    684                         m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
    685                         m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
    686                         if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
    687                                 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
     694                if(landmask->GetValue(i,0)==2){
     695                        if(nptflag->GetValue(i,0)==0){
     696                                u=double(rand())/double(RAND_MAX);
     697                                y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
     698                                m1->SetValue(i,1,m0->GetValue(i,1)+y*ptval);
    688699                                m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
    689                         }
    690                         if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
    691                                 m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
    692                                 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
     700                                if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
     701                                        m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
     702                                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
     703                                }
     704                                if(m1->GetValue(i,1)>=m1->GetValue(i,0)+mmax){
     705                                        m1->SetValue(i,1,m1->GetValue(i,0)+mmax);
     706                                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
     707                                }
     708                        }
     709                }
     710                else if(landmask->GetValue(i,0)==0){
     711                        if(nptflag->GetValue(i,0)==0){
     712                                u=double(rand())/double(RAND_MAX);
     713                                y=signe(u-0.5)*T*(pow(1+1/T,fabs(2*u-1))-1);
     714                                m1->SetValue(i,2,m0->GetValue(i,2)+y*ptval);
     715                                if(m1->GetValue(i,2)<=m1->GetValue(i,0)){
     716                                        m1->SetValue(i,2,m1->GetValue(i,0)+1e-10);
     717                                }
     718                                if(m1->GetValue(i,2)>=m1->GetValue(i,0)+mmax2){
     719                                        m1->SetValue(i,2,m1->GetValue(i,0)+mmax2);
     720                                }
    693721                        }
    694722                }
    … …  
    705733
    706734        for (int i=0;i<mx*my;i++){
    707                 if(nptflag->GetValue(i,0)==0){
    708                         m1->SetValue(i,1,m1col->GetValue(i,0));
    709                         m1->SetValue(i,2,m1col2->GetValue(i,0));
    710                         if(m1->GetValue(i,1)<m1->GetValue(i,0)){
    711                                 m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
    712                                 m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
    713                         }
    714                 }
    715                 else{
    716                         m1->SetValue(i,1,m0->GetValue(i,1));
    717                         m1->SetValue(i,2,m0->GetValue(i,2));
     735                if(landmask->GetValue(i,0)==2){
     736                        if(nptflag->GetValue(i,0)==0){
     737                                m1->SetValue(i,1,m1col->GetValue(i,0));
     738                                m1->SetValue(i,2,m1col2->GetValue(i,0));
     739                                if(m1->GetValue(i,1)<=m1->GetValue(i,0)){
     740                                        m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
     741                                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
     742                                }
     743                                if(fabs(m1->GetValue(i,2)-m1->GetValue(i,1))>1){
     744                                        m1->SetValue(i,2,m1->GetValue(i,1)+1e-10);
     745                                }
     746                        }
     747                        else{
     748                                m1->SetValue(i,1,m0->GetValue(i,1));
     749                                m1->SetValue(i,2,m0->GetValue(i,2));
     750                        }
     751                }
     752                else if(landmask->GetValue(i,0)==0){
     753                        if(nptflag->GetValue(i,0)==0){
     754                                m1->SetValue(i,2,m1col2->GetValue(i,0));
     755                                if(m1->GetValue(i,2)<=m1->GetValue(i,0)){
     756                                        m1->SetValue(i,2,m1->GetValue(i,0)+1e-10);
     757                                }
     758                                if(fabs(m1->GetValue(i,0)-m1->GetValue(i,1))>1){
     759                                        m1->SetValue(i,1,m1->GetValue(i,0)+1e-10);
     760                                }
     761                        }
     762                        else{
     763                                m1->SetValue(i,1,m0->GetValue(i,1));
     764                                m1->SetValue(i,2,m0->GetValue(i,2));
     765                        }
     766                }
     767                else {
     768                        if(nptflag->GetValue(i,0)==0){
     769                                if(fabs(m1->GetValue(i,0)-m1->GetValue(i,1))>1){
     770                                        m1->SetValue(i,1,m1->GetValue(i,0));
     771                                }
     772                                if(fabs(m1->GetValue(i,0)-m1->GetValue(i,2))>1){
     773                                        m1->SetValue(i,2,m1->GetValue(i,0));
     774                                }
     775                        }
     776                        else{
     777                                m1->SetValue(i,1,m0->GetValue(i,1));
     778                                m1->SetValue(i,2,m0->GetValue(i,2));
     779                        }
    718780                }
    719781        }
Note: See TracChangeset for help on using the changeset viewer.