Changeset 19107


Ignore:
Timestamp:
02/12/15 16:50:45 (10 years ago)
Author:
Mathieu Morlighem
Message:

CHG: reverted LevelsetAnalysis.cpp back

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r19101 r19107  
    123123}/*}}}*/
    124124ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
    125         _error_("Not implemented yet");
     125
     126        if(!element->IsOnBase()) return NULL;
     127        Element* basalelement = element->SpawnBasalElement();
     128
     129        /*Intermediaries */
     130        int  stabilization,dim, domaintype, calvinglaw;
     131        bool iscalving;
     132        int i, row, col;
     133        IssmDouble kappa;
     134        IssmDouble Jdet, dt, D_scalar;
     135        IssmDouble h,hx,hy,hz;
     136        IssmDouble vel;
     137        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate;
     138        IssmDouble* xyz_list = NULL;
     139
     140        /*Get problem dimension and whether there is calving or not*/
     141        basalelement->FindParam(&iscalving,TransientIscalvingEnum);
     142        basalelement->FindParam(&domaintype,DomainTypeEnum);
     143        basalelement->FindParam(&calvinglaw,CalvingLawEnum);
     144        basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
     145        switch(domaintype){
     146                case Domain2DverticalEnum:   dim = 1; break;
     147                case Domain2DhorizontalEnum: dim = 2; break;
     148                case Domain3DEnum:           dim = 2; break;
     149                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     150        }
     151
     152        /*Fetch number of nodes and dof for this finite element*/
     153        int numnodes    = basalelement->GetNumberOfNodes();
     154
     155        /*Initialize Element vector and other vectors*/
     156        ElementMatrix* Ke       = basalelement->NewElementMatrix();
     157        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
     158        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
     159        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
     160        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
     161        IssmDouble*    v        = xNew<IssmDouble>(dim);
     162        IssmDouble*    w        = xNew<IssmDouble>(dim);
     163        IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
     164        IssmDouble*    m        = xNewZeroInit<IssmDouble>(dim);
     165        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
     166
     167        /*Retrieve all inputs and parameters*/
     168        basalelement->GetVerticesCoordinates(&xyz_list);
     169        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     170        Input* vx_input           = NULL;
     171        Input* vy_input           = NULL;
     172        Input* calvingratex_input = NULL;
     173        Input* calvingratey_input = NULL;
     174        Input* lsf_slopex_input   = NULL;
     175        Input* lsf_slopey_input   = NULL;
     176        Input* calvingrate_input  = NULL;
     177        Input* meltingrate_input  = NULL;
     178
     179        /*Load velocities*/
     180        switch(domaintype){
     181                case Domain2DverticalEnum:
     182                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     183                        break;
     184                case Domain2DhorizontalEnum:
     185                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     186                        vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     187                        break;
     188                case Domain3DEnum:
     189                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
     190                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
     191                        break;
     192                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     193        }
     194
     195        /*Load calving inputs*/
     196        if(iscalving){
     197                switch(calvinglaw){
     198                        case DefaultCalvingEnum:
     199                                lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     200                                if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     201                                calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     202                                meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     203                                break;
     204                        case CalvingLevermannEnum:
     205                                switch(domaintype){
     206                                        case Domain2DverticalEnum:
     207                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     208                                                break;
     209                                        case Domain2DhorizontalEnum:
     210                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     211                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     212                                                break;
     213                                        case Domain3DEnum:
     214                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     215                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     216                                                break;
     217                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     218                                }
     219                                meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
     220                                break;
     221                        case CalvingPiEnum:
     222                                switch(domaintype){
     223                                        case Domain2DverticalEnum:
     224                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     225                                                break;
     226                                        case Domain2DhorizontalEnum:
     227                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     228                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     229                                                break;
     230                                        case Domain3DEnum:
     231                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     232                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     233                                                break;
     234                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     235                                }
     236                                meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum);     _assert_(meltingrate_input);
     237                                break;
     238                        case CalvingDevEnum:
     239                                switch(domaintype){
     240                                        case Domain2DverticalEnum:
     241                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     242                                                break;
     243                                        case Domain2DhorizontalEnum:
     244                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     245                                                calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     246                                                break;
     247                                        case Domain3DEnum:
     248                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     249                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     250                                                break;
     251                                        default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     252                                }
     253                                meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     254                                break;
     255                        default:
     256                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     257                }
     258        }
     259
     260        /* Start  looping on the number of gaussian points: */
     261        Gauss* gauss=basalelement->NewGauss(2);
     262        for(int ig=gauss->begin();ig<gauss->end();ig++){
     263                gauss->GaussPoint(ig);
     264
     265                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     266                D_scalar=gauss->weight*Jdet;
     267
     268                /* Transient */
     269                if(dt!=0.){
     270                        basalelement->NodalFunctions(basis,gauss);
     271                        TripleMultiply(basis,numnodes,1,0,
     272                                                &D_scalar,1,1,0,
     273                                                basis,1,numnodes,0,
     274                                                &Ke->values[0],1);
     275                        D_scalar*=dt;
     276                }
     277
     278                /* Advection */
     279                GetB(B,basalelement,xyz_list,gauss);
     280                GetBprime(Bprime,basalelement,xyz_list,gauss);
     281                vx_input->GetInputValue(&v[0],gauss);
     282                vy_input->GetInputValue(&v[1],gauss);
     283
     284                /*Get calving speed*/
     285                if(iscalving){
     286                        switch(calvinglaw){
     287                                case DefaultCalvingEnum:
     288                                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     289                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     290                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
     291                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
     292
     293                                        norm_dlsf=0.;
     294                                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     295                                        norm_dlsf=sqrt(norm_dlsf);
     296
     297                                        if(norm_dlsf>1.e-10)
     298                                         for(i=0;i<dim;i++){
     299                                          c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
     300                                         }
     301                                        else
     302                                         for(i=0;i<dim;i++){
     303                                                 c[i]=0.; m[i]=0.;
     304                                         }
     305                                        break;
     306
     307                                case CalvingLevermannEnum:
     308                                case CalvingPiEnum:
     309                                case CalvingDevEnum:
     310                                        calvingratex_input->GetInputValue(&c[0],gauss);
     311                                        if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
     312                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
     313                                        norm_calving=0.;
     314                                        for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
     315                                        norm_calving=sqrt(norm_calving)+1.e-14;
     316                                        for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     317                                        break;
     318
     319                                default:
     320                                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     321                        }
     322                }
     323
     324                /*Levelset speed is ice velocity - calving rate*/
     325                for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
     326
     327                /*Compute D*/
     328                for(row=0;row<dim;row++){
     329                        for(col=0;col<dim;col++){
     330                                if(row==col)
     331                                 D[row*dim+col]=D_scalar*w[row];
     332                                else
     333                                 D[row*dim+col]=0.;
     334                        }
     335                }
     336
     337                TripleMultiply(B,dim,numnodes,1,
     338                                        D,dim,dim,0,
     339                                        Bprime,dim,numnodes,0,
     340                                        &Ke->values[0],1);
     341
     342                /* Stabilization */
     343                vel=0.;
     344                for(i=0;i<dim;i++) vel+=w[i]*w[i];
     345                vel=sqrt(vel)+1.e-14;
     346                switch(stabilization){
     347                        case 0:
     348                                // no stabilization, do nothing
     349                                break;
     350                        case 1:
     351                                /* Artificial Diffusion */
     352                                basalelement->ElementSizes(&hx,&hy,&hz);
     353                                h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
     354                                kappa=h*vel/2.;
     355                                for(row=0;row<dim;row++)
     356                                        for(col=0;col<dim;col++)
     357                                        if(row==col)
     358                                                D[row*dim+col]=D_scalar*kappa;
     359                                        else
     360                                                D[row*dim+col]=0.;
     361
     362                                TripleMultiply(Bprime,dim,numnodes,1,
     363                                                        D,dim,dim,0,
     364                                                        Bprime,dim,numnodes,0,
     365                                                        &Ke->values[0],1);
     366                                break; 
     367                        case 2:
     368                                /* Streamline Upwinding */
     369                                basalelement->ElementSizes(&hx,&hy,&hz);
     370                                h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
     371                                for(row=0;row<dim;row++)
     372                                        for(col=0;col<dim;col++)
     373                                                D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
     374
     375                                TripleMultiply(Bprime,dim,numnodes,1,
     376                                                        D,dim,dim,0,
     377                                                        Bprime,dim,numnodes,0,
     378                                                        &Ke->values[0],1);
     379                                break;
     380                        default:
     381                                _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
     382                }
     383        }
     384
     385        /*Clean up and return*/
     386        xDelete<IssmDouble>(xyz_list);
     387        xDelete<IssmDouble>(basis);
     388        xDelete<IssmDouble>(B);
     389        xDelete<IssmDouble>(D);
     390        xDelete<IssmDouble>(Bprime);
     391        xDelete<IssmDouble>(v);
     392        xDelete<IssmDouble>(w);
     393        xDelete<IssmDouble>(c);
     394        xDelete<IssmDouble>(m);
     395        xDelete<IssmDouble>(dlsf);
     396        delete gauss;
     397        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     398        return Ke;
    126399}/*}}}*/
    127400ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
    128401       
    129         _error_("Not implemented yet");
     402        if(!element->IsOnBase()) return NULL;
     403        Element* basalelement = element->SpawnBasalElement();
     404
     405        /*Intermediaries */
     406        int i, ig, domaintype;
     407        IssmDouble  Jdet,dt;
     408        IssmDouble  lsf;
     409        IssmDouble* xyz_list = NULL;
     410       
     411        /*Fetch number of nodes and dof for this finite element*/
     412        int numnodes = basalelement->GetNumberOfNodes();
     413
     414        /*Initialize Element vector*/
     415        ElementVector* pe = basalelement->NewElementVector();
     416        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     417       
     418        if(dt!=0.){
     419                /*Initialize basis vector*/
     420                IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     421
     422                /*Retrieve all inputs and parameters*/
     423                basalelement->GetVerticesCoordinates(&xyz_list);
     424                Input* levelset_input     = basalelement->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
     425
     426                /* Start  looping on the number of gaussian points: */
     427                Gauss* gauss=basalelement->NewGauss(2);
     428                for(ig=gauss->begin();ig<gauss->end();ig++){
     429                        gauss->GaussPoint(ig);
     430
     431                        basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     432                        basalelement->NodalFunctions(basis,gauss);
     433
     434                        /* old function value */
     435                        levelset_input->GetInputValue(&lsf,gauss);
     436                        for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
     437                }
     438
     439                /*Clean up and return*/
     440                xDelete<IssmDouble>(xyz_list);
     441                xDelete<IssmDouble>(basis);
     442                basalelement->FindParam(&domaintype,DomainTypeEnum);
     443                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     444                delete gauss;
     445        }
     446
     447        return pe;
    130448}/*}}}*/
    131449void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.