Changeset 18351


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

CHG: added Richardson subroutine

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18350 r18351  
    159159
    160160}/*}}}*/
    161 #endif
    162 void solutionsequence_fct(FemModel* femmodel){
    163 
    164         /*intermediary: */
    165         Vector<IssmDouble>*  Ml = NULL;
    166         Matrix<IssmDouble>*  K  = NULL;
    167         Matrix<IssmDouble>*  Mc = NULL;
    168         Vector<IssmDouble>*  ug = NULL;
    169         Vector<IssmDouble>*  uf = NULL;
    170 
    171         IssmDouble theta,deltat,dmax;
    172         int        dof,ncols,ncols2,rstart,rend;
    173         int        configuration_type;
    174         double     d;
    175         int*       cols  = NULL;
    176         int*       cols2 = NULL;
    177         double*    vals  = NULL;
    178         double*    vals2 = NULL;
    179 
    180         /*Create analysis*/
    181         MasstransportAnalysis* analysis = new MasstransportAnalysis();
    182 
    183         /*Recover parameters: */
    184         femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
    185         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    186         femmodel->UpdateConstraintsx();
    187         theta = 0.5;
    188 
    189         /*Create lumped mass matrix*/
    190         analysis->LumpedMassMatrix(&Ml,femmodel);
    191         analysis->MassMatrix(&Mc,femmodel);
    192         analysis->FctKMatrix(&K,NULL,femmodel);
    193 
    194         /*Convert matrices to PETSc matrices*/
    195         Mat D_petsc  = NULL;
    196         Mat LHS      = NULL;
    197         Vec RHS      = NULL;
    198         Vec u        = NULL;
    199         Mat K_petsc  = K->pmatrix->matrix;
    200         Vec Ml_petsc = Ml->pvector->vector;
    201         Mat Mc_petsc = Mc->pmatrix->matrix;
    202 
    203         /*Create D Matrix*/
    204         #ifdef _HAVE_PETSC_
    205         CreateDMatrix(&D_petsc,K_petsc);
    206 
    207         /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
    208         CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
    209 
    210         /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
    211         GetSolutionFromInputsx(&ug,femmodel);
    212         Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
    213         delete ug;
    214         CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
    215         delete uf;
    216 
    217         /*Go solve!*/
    218         SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
    219         MatFree(&LHS);
    220         VecFree(&RHS);
    221 
    222         /*Richardson to calculate udot*/
    223         Vec udot = NULL;
     161void RichardsonUdot(Vec* pudot,Vec u,Vec Ml,Mat K,Mat Mc){/*{{{*/
     162        /*Use Richardson's formulato get udot using 5 steps and an initial guess of 0
     163         *
     164         * udot_new = udot_old + Ml^-1 (K^(n+1) u - Mc udot_old)
     165         *
     166         * */
     167
     168        /*Intermediaries*/
     169        Vec udot  = NULL;
     170        Vec temp1 = NULL;
     171        Vec temp2 = NULL;
     172
     173        /*Initialize vectors*/
    224174        VecDuplicate(u,&udot);
     175        VecDuplicate(u,&temp1);
     176        VecDuplicate(u,&temp2);
     177
     178        /*Initial guess*/
    225179        VecZeroEntries(udot);
    226         Vec temp1 = NULL; VecDuplicate(u,&temp1);
    227         Vec temp2 = NULL; VecDuplicate(u,&temp2);
     180       
     181        /*Richardson's iterations*/
    228182        for(int i=0;i<5;i++){
    229                 /*udot_new = udot_old + Ml^-1 (K^(n+1) u -- Mc udot_old)*/
    230                 MatMult(Mc_petsc,udot,temp1);
    231                 MatMult(K_petsc, u,   temp2);
    232                 VecAXPBY(temp2,-1.,1.,temp1); // temp2 = (K^(n+1) u -- Mc udot_old)
    233                 VecPointwiseDivide(temp1,temp2,Ml_petsc); //temp1 = Ml^-1 temp2
     183                MatMult(Mc,udot,temp1);
     184                MatMult(K, u,   temp2);
     185                VecAXPBY(temp2,-1.,1.,temp1);       // temp2 = (K^(n+1) u -- Mc udot_old)
     186                VecPointwiseDivide(temp1,temp2,Ml); // temp1 = Ml^-1 temp2
    234187                VecAXPBY(udot,1.,1.,temp1);
    235188        }
     189
     190        /*Clean up and assign output pointer*/
    236191        VecFree(&temp1);
    237192        VecFree(&temp2);
    238 
    239         /*Serialize u and udot*/
    240         IssmDouble* udot_serial = NULL;
    241         IssmDouble* u_serial    = NULL;
    242         IssmDouble* ml_serial    = NULL;
    243         VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
    244         VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
    245         VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
    246 
    247         /*Anti diffusive fluxes*/
    248         Vec Ri_plus  = NULL;
    249         Vec Ri_minus = NULL;
    250         double uiLmax = 3.;
    251         double uiLmin = 2.;
    252         VecDuplicate(u,&Ri_plus);
    253         VecDuplicate(u,&Ri_minus);
    254         MatGetOwnershipRange(K_petsc,&rstart,&rend);
    255         for(int row=rstart; row<rend; row++){
    256                 double Pi_plus  = 0.;
    257                 double Pi_minus = 0.;
    258                 MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
    259                 MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    260                 _assert_(ncols==ncols2);
    261                 for(int j=0; j<ncols; j++) {
    262                         _assert_(cols[j]==cols2[j]);
    263                         d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
    264                         if(row!=cols[j]){
    265                                 if(d>0.){
    266                                         Pi_plus  += d;
     193        *pudot=udot;
     194
     195}/*}}}*/
     196#endif
     197        void solutionsequence_fct(FemModel* femmodel){
     198
     199                /*intermediary: */
     200                Vector<IssmDouble>*  Ml = NULL;
     201                Matrix<IssmDouble>*  K  = NULL;
     202                Matrix<IssmDouble>*  Mc = NULL;
     203                Vector<IssmDouble>*  ug = NULL;
     204                Vector<IssmDouble>*  uf = NULL;
     205
     206                IssmDouble theta,deltat,dmax;
     207                int        dof,ncols,ncols2,rstart,rend;
     208                int        configuration_type;
     209                double     d;
     210                int*       cols  = NULL;
     211                int*       cols2 = NULL;
     212                double*    vals  = NULL;
     213                double*    vals2 = NULL;
     214
     215                /*Create analysis*/
     216                MasstransportAnalysis* analysis = new MasstransportAnalysis();
     217
     218                /*Recover parameters: */
     219                femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
     220                femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     221                femmodel->UpdateConstraintsx();
     222                theta = 0.5;
     223
     224                /*Create lumped mass matrix*/
     225                analysis->LumpedMassMatrix(&Ml,femmodel);
     226                analysis->MassMatrix(&Mc,femmodel);
     227                analysis->FctKMatrix(&K,NULL,femmodel);
     228
     229                /*Convert matrices to PETSc matrices*/
     230                Mat D_petsc  = NULL;
     231                Mat LHS      = NULL;
     232                Vec RHS      = NULL;
     233                Vec u        = NULL;
     234                Vec udot     = NULL;
     235                Mat K_petsc  = K->pmatrix->matrix;
     236                Vec Ml_petsc = Ml->pvector->vector;
     237                Mat Mc_petsc = Mc->pmatrix->matrix;
     238
     239                /*Create D Matrix*/
     240                #ifdef _HAVE_PETSC_
     241                CreateDMatrix(&D_petsc,K_petsc);
     242
     243                /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
     244                CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
     245
     246                /*Get previous solution u^n*/
     247                GetSolutionFromInputsx(&ug,femmodel);
     248                Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
     249                delete ug;
     250
     251                /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
     252                CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
     253                delete uf;
     254
     255                /*Go solve!*/
     256                SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
     257                MatFree(&LHS);
     258                VecFree(&RHS);
     259
     260                /*Richardson to calculate udot*/
     261                RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
     262
     263                /*Serialize u and udot*/
     264                IssmDouble* udot_serial = NULL;
     265                IssmDouble* u_serial    = NULL;
     266                IssmDouble* ml_serial    = NULL;
     267                VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
     268                VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
     269                VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
     270
     271                /*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*/
     317                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;
    267338                                }
    268339                                else{
    269                                         Pi_minus += d;
     340                                        d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
    270341                                }
    271342                        }
    272                 }
    273 
    274                 /*Compute Qis and Ris*/
    275                 double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
    276                 double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
    277                 d = 1.;
    278                 if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
    279                 VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
    280                 d = 1.;
    281                 if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
    282                 VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
    283 
    284                 MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
    285                 MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    286         }
    287         VecAssemblyBegin(Ri_plus);
    288         VecAssemblyEnd(  Ri_plus);
    289         VecAssemblyBegin(Ri_minus);
    290         VecAssemblyEnd(  Ri_minus);
    291 
    292         /*Serialize Ris*/
    293         IssmDouble* Ri_plus_serial  = NULL;
    294         IssmDouble* Ri_minus_serial = NULL;
    295         VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
    296         VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
    297         VecFree(&Ri_plus);
    298         VecFree(&Ri_minus);
    299 
    300         /*Build fbar*/
    301         Vec Fbar = NULL;
    302         VecDuplicate(u,&Fbar);
    303         for(int row=rstart; row<rend; row++){
    304                 MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
    305                 MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    306                 _assert_(ncols==ncols2);
    307                 d = 0.;
    308                 for(int j=0; j<ncols; j++) {
    309                         _assert_(cols[j]==cols2[j]);
    310                         if(row==cols[j]) continue;
    311                         double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
    312                         if(f_ij>0){
    313                                 d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
    314                         }
    315                         else{
    316                                 d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
    317                         }
    318                 }
    319                 VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
    320                 MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
    321                 MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
    322         }
    323         VecAssemblyBegin(Fbar);
    324         VecAssemblyEnd(  Fbar);
    325 
    326         MatFree(&D_petsc);
    327         delete Mc;
    328         xDelete<IssmDouble>(udot_serial);
    329         xDelete<IssmDouble>(u_serial);
    330         xDelete<IssmDouble>(ml_serial);
    331         xDelete<IssmDouble>(Ri_plus_serial);
    332         xDelete<IssmDouble>(Ri_minus_serial);
    333 
    334         /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
    335         VecDuplicate(u,&temp1);
    336         VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
    337         VecAXPBY(udot,1.,1.,temp1);
    338         VecAXPY(u,deltat,temp1);
    339         VecFree(&Fbar);
    340         VecFree(&udot);
    341         VecFree(&temp1);
    342 
    343         uf =new Vector<IssmDouble>(u);
    344         VecFree(&u);
    345 
    346         InputUpdateFromSolutionx(femmodel,uf);
    347         delete uf;
    348 
    349         #else
    350         _error_("PETSc needs to be installed");
    351         #endif
    352 
    353         delete Ml;
    354         delete K;
    355         delete analysis;
    356 
    357 }
     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
     350                MatFree(&D_petsc);
     351                delete Mc;
     352                xDelete<IssmDouble>(udot_serial);
     353                xDelete<IssmDouble>(u_serial);
     354                xDelete<IssmDouble>(ml_serial);
     355                xDelete<IssmDouble>(Ri_plus_serial);
     356                xDelete<IssmDouble>(Ri_minus_serial);
     357
     358                /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
     359                Vec temp1 = NULL; VecDuplicate(u,&temp1);
     360                VecDuplicate(u,&temp1);
     361                VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
     362                VecAXPBY(udot,1.,1.,temp1);
     363                VecAXPY(u,deltat,temp1);
     364                VecFree(&Fbar);
     365                VecFree(&udot);
     366                VecFree(&temp1);
     367
     368                uf =new Vector<IssmDouble>(u);
     369                VecFree(&u);
     370
     371                InputUpdateFromSolutionx(femmodel,uf);
     372                delete uf;
     373
     374                #else
     375                _error_("PETSc needs to be installed");
     376                #endif
     377
     378                delete Ml;
     379                delete K;
     380                delete analysis;
     381
     382        }
Note: See TracChangeset for help on using the changeset viewer.