Changeset 18355


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

CHG: trying to reintroduce ulmin ulmax wihtout success

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18353 r18355  
    194194
    195195}/*}}}*/
    196 void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble deltat){/*{{{*/
     196void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble* ulmin,IssmDouble* ulmax,IssmDouble deltat){/*{{{*/
    197197
    198198        /*Intermediaries*/
    199199        Vec         Ri_plus  = NULL;
    200200        Vec         Ri_minus = NULL;
    201         IssmDouble  uiLmax   = 3.;
    202         IssmDouble  uiLmin   = 2.;
    203201        int         ncols,ncols2,rstart,rend;
    204202        double      d;
    … …  
    234232
    235233                /*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]);
     234                double Qi_plus  = ml_serial[row]/deltat*(3. - u[row]);
     235                double Qi_minus = ml_serial[row]/deltat*(2. - u[row]);
     236                //double Qi_plus  = ml_serial[row]/deltat*(ulmax[row] - u[row]);
     237                //double Qi_minus = ml_serial[row]/deltat*(ulmin[row] - u[row]);
    238238                d = 1.;
    239239                if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
    … …  
    303303        *pFbar = Fbar;
    304304}/*}}}*/
     305void UpdateSolution(Vec u,Vec udot,Vec Ml,Vec Fbar,IssmDouble deltat){/*{{{*/
     306
     307        /*Intermediary*/
     308        Vec temp1 = NULL;
     309
     310        /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
     311        VecDuplicate(u,&temp1);
     312        VecPointwiseDivide(temp1,Fbar,Ml); //temp1 = Ml^-1 temp2
     313        VecAXPBY(udot,1.,1.,temp1);
     314        VecAXPY(u,deltat,temp1);
     315
     316        /*CLean up and return*/
     317        VecFree(&temp1);
     318}/*}}}*/
    305319#endif
    306         void solutionsequence_fct(FemModel* femmodel){
    307 
    308                 /*intermediary: */
    309                 Vector<IssmDouble>*  Ml = NULL;
    310                 Matrix<IssmDouble>*  K  = NULL;
    311                 Matrix<IssmDouble>*  Mc = NULL;
    312                 Vector<IssmDouble>*  ug = NULL;
    313                 Vector<IssmDouble>*  uf = NULL;
    314 
    315                 IssmDouble theta,deltat,dmax;
    316                 int        dof,ncols,ncols2,rstart,rend;
    317                 int        configuration_type;
    318                 double     d;
    319                 int*       cols  = NULL;
    320                 int*       cols2 = NULL;
    321                 double*    vals  = NULL;
    322                 double*    vals2 = NULL;
    323 
    324                 /*Create analysis*/
    325                 MasstransportAnalysis* analysis = new MasstransportAnalysis();
    326 
    327                 /*Recover parameters: */
    328                 femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
    329                 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    330                 femmodel->UpdateConstraintsx();
    331                 theta = 0.5;
    332 
    333                 /*Create lumped mass matrix*/
    334                 analysis->LumpedMassMatrix(&Ml,femmodel);
    335                 analysis->MassMatrix(&Mc,femmodel);
    336                 analysis->FctKMatrix(&K,NULL,femmodel);
    337 
    338                 /*Convert matrices to PETSc matrices*/
    339                 Mat D_petsc  = NULL;
    340                 Mat LHS      = NULL;
    341                 Vec RHS      = NULL;
    342                 Vec u        = NULL;
    343                 Vec udot     = NULL;
    344                 Mat K_petsc  = K->pmatrix->matrix;
    345                 Vec Ml_petsc = Ml->pvector->vector;
    346                 Mat Mc_petsc = Mc->pmatrix->matrix;
    347 
    348                 /*Create D Matrix*/
    349                 #ifdef _HAVE_PETSC_
    350                 CreateDMatrix(&D_petsc,K_petsc);
    351 
    352                 /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
    353                 CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
    354 
    355                 /*Get previous solution u^n*/
    356                 GetSolutionFromInputsx(&ug,femmodel);
    357                 Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
    358                 delete ug;
    359 
    360                 /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
    361                 CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
    362                 delete uf;
    363 
    364                 /*Go solve lower order solution*/
    365                 SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
    366                 MatFree(&LHS);
    367                 VecFree(&RHS);
    368 
    369                 /*Richardson to calculate udot*/
    370                 RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
    371 
    372                 /*Serialize u and udot*/
    373                 IssmDouble* udot_serial = NULL;
    374                 IssmDouble* u_serial    = NULL;
    375                 IssmDouble* ml_serial   = NULL;
    376                 VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
    377                 VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
    378                 VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
    379 
    380                 /*Anti diffusive fluxes*/
    381                 Vec         Fbar            = NULL;
    382                 IssmDouble* Ri_minus_serial = NULL;
    383                 IssmDouble* Ri_plus_serial  = NULL;
    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*/
    388                 MatFree(&D_petsc);
    389                 delete Mc;
    390                 xDelete<IssmDouble>(udot_serial);
    391                 xDelete<IssmDouble>(u_serial);
    392                 xDelete<IssmDouble>(ml_serial);
    393                 xDelete<IssmDouble>(Ri_plus_serial);
    394                 xDelete<IssmDouble>(Ri_minus_serial);
    395 
    396                 /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
    397                 Vec temp1 = NULL;
    398                 VecDuplicate(u,&temp1);
    399                 VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2
    400                 VecAXPBY(udot,1.,1.,temp1);
    401                 VecAXPY(u,deltat,temp1);
    402                 VecFree(&Fbar);
    403                 VecFree(&udot);
    404                 VecFree(&temp1);
    405 
    406                 uf =new Vector<IssmDouble>(u);
    407                 VecFree(&u);
    408 
    409                 InputUpdateFromSolutionx(femmodel,uf);
    410                 delete uf;
    411 
    412                 #else
    413                 _error_("PETSc needs to be installed");
    414                 #endif
    415 
    416                 delete Ml;
    417                 delete K;
    418                 delete analysis;
    419 
    420         }
     320void solutionsequence_fct(FemModel* femmodel){
     321
     322        /*intermediary: */
     323        IssmDouble           theta,deltat,dmax;
     324        int                  configuration_type;
     325        Vector<IssmDouble>*  Ml = NULL;
     326        Matrix<IssmDouble>*  K  = NULL;
     327        Matrix<IssmDouble>*  Mc = NULL;
     328        Vector<IssmDouble>*  ug = NULL;
     329        Vector<IssmDouble>*  uf = NULL;
     330
     331        /*Create analysis*/
     332        MasstransportAnalysis* analysis = new MasstransportAnalysis();
     333
     334        /*Recover parameters: */
     335        femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
     336        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     337        femmodel->UpdateConstraintsx();
     338        theta = 0.5;
     339
     340        /*Create lumped mass matrix*/
     341        analysis->LumpedMassMatrix(&Ml,femmodel);
     342        analysis->MassMatrix(&Mc,femmodel);
     343        analysis->FctKMatrix(&K,NULL,femmodel);
     344        delete analysis;
     345
     346        /*Convert matrices to PETSc matrices*/
     347        Mat D_petsc  = NULL;
     348        Mat LHS      = NULL;
     349        Vec RHS      = NULL;
     350        Vec u        = NULL;
     351        Vec udot     = NULL;
     352        Mat K_petsc  = K->pmatrix->matrix;
     353        Vec Ml_petsc = Ml->pvector->vector;
     354        Mat Mc_petsc = Mc->pmatrix->matrix;
     355
     356        /*Create D Matrix*/
     357        #ifdef _HAVE_PETSC_
     358        CreateDMatrix(&D_petsc,K_petsc);
     359
     360        /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
     361        CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
     362
     363        /*Get previous solution u^n*/
     364        GetSolutionFromInputsx(&ug,femmodel);
     365        Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
     366        delete ug;
     367
     368        /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
     369        CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
     370        delete uf;
     371
     372        /*Go solve lower order solution*/
     373        SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
     374        MatFree(&LHS);
     375        VecFree(&RHS);
     376
     377        /*Richardson to calculate udot*/
     378        RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
     379        delete K;
     380
     381        /*Serialize u and udot*/
     382        IssmDouble* udot_serial = NULL;
     383        IssmDouble* u_serial    = NULL;
     384        IssmDouble* ml_serial   = NULL;
     385        VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
     386        VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
     387        VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
     388
     389        /*Anti diffusive fluxes*/
     390        Vec         Fbar            = NULL;
     391        IssmDouble *Ri_minus_serial = NULL;
     392        IssmDouble *Ri_plus_serial  = NULL;
     393        IssmDouble *ulmin           = NULL;
     394        IssmDouble *ulmax           = NULL;
     395        femmodel->GetInputLocalMinMaxOnNodesx(&ulmin,&ulmax,ThicknessEnum);
     396        CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,ulmin,ulmax,deltat);
     397        CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
     398        xDelete<IssmDouble>(Ri_plus_serial);
     399        xDelete<IssmDouble>(Ri_minus_serial);
     400        xDelete<IssmDouble>(ulmin);
     401        xDelete<IssmDouble>(ulmax);
     402
     403        /*Clean up*/
     404        MatFree(&D_petsc);
     405        delete Mc;
     406        xDelete<IssmDouble>(udot_serial);
     407        xDelete<IssmDouble>(u_serial);
     408        xDelete<IssmDouble>(ml_serial);
     409
     410        /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
     411        UpdateSolution(u,udot,Ml_petsc,Fbar,deltat);
     412        uf =new Vector<IssmDouble>(u);
     413        VecFree(&u);
     414        VecFree(&Fbar);
     415        VecFree(&udot);
     416        delete Ml;
     417
     418        /*Update Element inputs*/
     419        InputUpdateFromSolutionx(femmodel,uf);
     420        delete uf;
     421
     422        #else
     423        _error_("PETSc needs to be installed");
     424        #endif
     425}
Note: See TracChangeset for help on using the changeset viewer.