Changeset 19101


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

TEMP: removed levelset code temporarily

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19099 r19101  
    123123}/*}}}*/
    124124ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
    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;
     125        _error_("Not implemented yet");
    399126}/*}}}*/
    400127ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
    401128       
    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;
     129        _error_("Not implemented yet");
    448130}/*}}}*/
    449131void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.