Changeset 18353


Ignore:
Timestamp:
08/08/14 14:34:03 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: one more subfunction

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

    r18351 r18353  
    194194
    195195}/*}}}*/
     196void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble deltat){/*{{{*/
     197
     198        /*Intermediaries*/
     199        Vec         Ri_plus  = NULL;
     200        Vec         Ri_minus = NULL;
     201        IssmDouble  uiLmax   = 3.;
     202        IssmDouble  uiLmin   = 2.;
     203        int         ncols,ncols2,rstart,rend;
     204        double      d;
     205        int        *cols     = NULL;
     206        int        *cols2    = NULL;
     207        double     *vals     = NULL;
     208        double     *vals2    = NULL;
     209
     210        /*Initialize vectors*/
     211        VecDuplicate(uvec,&Ri_plus);
     212        VecDuplicate(uvec,&Ri_minus);
     213
     214        /*Go through D and calculate Ris*/
     215        MatGetOwnershipRange(D,&rstart,&rend);
     216        for(int row=rstart; row<rend; row++){
     217                double Pi_plus  = 0.;
     218                double Pi_minus = 0.;
     219                MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
     220                MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     221                _assert_(ncols==ncols2);
     222                for(int j=0; j<ncols; j++) {
     223                        _assert_(cols[j]==cols2[j]);
     224                        d = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
     225                        if(row!=cols[j]){
     226                                if(d>0.){
     227                                        Pi_plus  += d;
     228                                }
     229                                else{
     230                                        Pi_minus += d;
     231                                }
     232                        }
     233                }
     234
     235                /*Compute Qis and Ris*/
     236                double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u[row]);
     237                double Qi_minus = ml_serial[row]/deltat*(uiLmin - u[row]);
     238                d = 1.;
     239                if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
     240                VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
     241                d = 1.;
     242                if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
     243                VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
     244
     245                MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
     246                MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     247        }
     248        VecAssemblyBegin(Ri_plus);
     249        VecAssemblyEnd(  Ri_plus);
     250        VecAssemblyBegin(Ri_minus);
     251        VecAssemblyEnd(  Ri_minus);
     252
     253        /*Serialize Ris*/
     254        IssmDouble* Ri_plus_serial  = NULL;
     255        IssmDouble* Ri_minus_serial = NULL;
     256        VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
     257        VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
     258        VecFree(&Ri_plus);
     259        VecFree(&Ri_minus);
     260
     261        /*Assign output pointer*/
     262        *pRi_plus_serial  = Ri_plus_serial;
     263        *pRi_minus_serial = Ri_minus_serial;
     264}/*}}}*/
     265void CreateFbar(Vec* pFbar,IssmDouble* Ri_plus,IssmDouble* Ri_minus,Mat Mc,Mat D,IssmDouble* udot,IssmDouble* u,Vec uvec){/*{{{*/
     266
     267        /*Intermediaries*/
     268        Vec         Fbar = NULL;
     269        int         ncols,ncols2,rstart,rend;
     270        double      d,f_ij;
     271        int        *cols     = NULL;
     272        int        *cols2    = NULL;
     273        double     *vals     = NULL;
     274        double     *vals2    = NULL;
     275
     276        /*Build fbar*/
     277        VecDuplicate(uvec,&Fbar);
     278        MatGetOwnershipRange(D,&rstart,&rend);
     279        for(int row=rstart; row<rend; row++){
     280                MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
     281                MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     282                _assert_(ncols==ncols2);
     283                d = 0.;
     284                for(int j=0; j<ncols; j++) {
     285                        _assert_(cols[j]==cols2[j]);
     286                        if(row==cols[j]) continue;
     287                        f_ij = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
     288                        if(f_ij>0){
     289                                d += min(Ri_plus[row],Ri_minus[cols[j]]) * f_ij;
     290                        }
     291                        else{
     292                                d += min(Ri_minus[row],Ri_plus[cols[j]]) * f_ij;
     293                        }
     294                }
     295                VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
     296                MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
     297                MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
     298        }
     299        VecAssemblyBegin(Fbar);
     300        VecAssemblyEnd(  Fbar);
     301
     302        /*Assign output pointer*/
     303        *pFbar = Fbar;
     304}/*}}}*/
    196305#endif
    197306        void solutionsequence_fct(FemModel* femmodel){
    … …  
    253362                delete uf;
    254363
    255                 /*Go solve!*/
     364                /*Go solve lower order solution*/
    256365                SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
    257366                MatFree(&LHS);
    … …  
    264373                IssmDouble* udot_serial = NULL;
    265374                IssmDouble* u_serial    = NULL;
    266                 IssmDouble* ml_serial    = NULL;
     375                IssmDouble* ml_serial   = NULL;
    267376                VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
    268377                VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
    … …  
    270379
    271380                /*Anti diffusive fluxes*/
    272                 Vec Ri_plus  = NULL;
    273                 Vec Ri_minus = NULL;
    274                 double uiLmax = 3.;
    275                 double uiLmin = 2.;
    276                 VecDuplicate(u,&Ri_plus);
    277                 VecDuplicate(u,&Ri_minus);
    278                 MatGetOwnershipRange(K_petsc,&rstart,&rend);
    279                 for(int row=rstart; row<rend; row++){
    280                         double Pi_plus  = 0.;
    281                         double Pi_minus = 0.;
    282                         MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
    283                         MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    284                         _assert_(ncols==ncols2);
    285                         for(int j=0; j<ncols; j++) {
    286                                 _assert_(cols[j]==cols2[j]);
    287                                 d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
    288                                 if(row!=cols[j]){
    289                                         if(d>0.){
    290                                                 Pi_plus  += d;
    291                                         }
    292                                         else{
    293                                                 Pi_minus += d;
    294                                         }
    295                                 }
    296                         }
    297 
    298                         /*Compute Qis and Ris*/
    299                         double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
    300                         double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
    301                         d = 1.;
    302                         if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
    303                         VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
    304                         d = 1.;
    305                         if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
    306                         VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
    307 
    308                         MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
    309                         MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    310                 }
    311                 VecAssemblyBegin(Ri_plus);
    312                 VecAssemblyEnd(  Ri_plus);
    313                 VecAssemblyBegin(Ri_minus);
    314                 VecAssemblyEnd(  Ri_minus);
    315 
    316                 /*Serialize Ris*/
     381                Vec         Fbar            = NULL;
     382                IssmDouble* Ri_minus_serial = NULL;
    317383                IssmDouble* Ri_plus_serial  = NULL;
    318                 IssmDouble* Ri_minus_serial = NULL;
    319                 VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
    320                 VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
    321                 VecFree(&Ri_plus);
    322                 VecFree(&Ri_minus);
    323 
    324                 /*Build fbar*/
    325                 Vec Fbar = NULL;
    326                 VecDuplicate(u,&Fbar);
    327                 for(int row=rstart; row<rend; row++){
    328                         MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
    329                         MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    330                         _assert_(ncols==ncols2);
    331                         d = 0.;
    332                         for(int j=0; j<ncols; j++) {
    333                                 _assert_(cols[j]==cols2[j]);
    334                                 if(row==cols[j]) continue;
    335                                 double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
    336                                 if(f_ij>0){
    337                                         d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
    338                                 }
    339                                 else{
    340                                         d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
    341                                 }
    342                         }
    343                         VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
    344                         MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
    345                         MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    346                 }
    347                 VecAssemblyBegin(Fbar);
    348                 VecAssemblyEnd(  Fbar);
    349 
     384                CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,deltat);
     385                CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
     386
     387                /*Clean up*/
    350388                MatFree(&D_petsc);
    351389                delete Mc;
    … …  
    357395
    358396                /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
    359                 Vec temp1 = NULL; VecDuplicate(u,&temp1);
     397                Vec temp1 = NULL;
    360398                VecDuplicate(u,&temp1);
    361399                VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
Note: See TracChangeset for help on using the changeset viewer.