Changeset 27572


Ignore:
Timestamp:
02/14/23 05:09:07 (2 years ago)
Author:
Mathieu Morlighem
Message:

CHG: minor changes to age analysis

Location:
issm/trunk-jpl/src/c/analyses
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r27154 r27572  
    9292}/*}}}*/
    9393ElementMatrix* AgeAnalysis::CreateKMatrix(Element* element){/*{{{*/
    94 
    95         _error_("STOP");
    96         /* Check if ice in element */
    97         if(!element->IsIceInElement()) return NULL;
    98 
    99         /*compute all stiffness matrices for this element*/
    100         ElementMatrix* Ke1=CreateKMatrixVolume(element);
    101         ElementMatrix* Ke2=CreateKMatrixShelf(element);
    102         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    103 
    104         /*clean-up and return*/
    105         delete Ke1;
    106         delete Ke2;
    107         return Ke;
    108 }/*}}}*/
    109 ElementMatrix* AgeAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
    110 
    111         /* Check if ice in element */
    112         if(!element->IsIceInElement()) return NULL;
    113 
    114         /*Initialize Element matrix and return if necessary*/
    115         if(!element->IsOnBase() || !element->IsAllFloating()) return NULL;
    116 
    117         IssmDouble  dt,Jdet,D;
    118         IssmDouble *xyz_list_base = NULL;
    119 
    120         /*Get basal element*/
    121         if(!element->IsOnBase() || !element->IsAllFloating()) return NULL;
    122 
    123         /*Fetch number of nodes for this finite element*/
    124         int numnodes = element->GetNumberOfNodes();
    125 
    126         /*Initialize vectors*/
    127         ElementMatrix* Ke    = element->NewElementMatrix();
    128         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    129 
    130         /*Retrieve all inputs and parameters*/
    131         element->GetVerticesCoordinatesBase(&xyz_list_base);
    132         element->FindParam(&dt,TimesteppingTimeStepEnum);
    133         IssmDouble gravity             = element->FindParam(ConstantsGEnum);
    134         IssmDouble rho_water           = element->FindParam(MaterialsRhoSeawaterEnum);
    135         IssmDouble rho_ice             = element->FindParam(MaterialsRhoIceEnum);
    136         IssmDouble heatcapacity        = element->FindParam(MaterialsHeatcapacityEnum);
    137         IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
    138 
    139         /* Start  looping on the number of gaussian points: */
    140         Gauss* gauss=element->NewGaussBase(4);
    141         while(gauss->next()){
    142 
    143                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    144                 element->NodalFunctions(basis,gauss);
    145 
    146                 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity/(heatcapacity*rho_ice);
    147                 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
    148                 for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
    149         }
    150 
    151         /*Clean up and return*/
    152         delete gauss;
    153         xDelete<IssmDouble>(basis);
    154         xDelete<IssmDouble>(xyz_list_base);
    155         return Ke;
    156 }/*}}}*/
    157 ElementMatrix* AgeAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
    15894
    15995        /* Check if ice in element */
    … …  
    303239        if(!element->IsIceInElement()) return NULL;
    304240
    305         /*compute all load vectors for this element*/
    306         ElementVector* pe1=CreatePVectorVolume(element);
    307         ElementVector* pe2=CreatePVectorSheet(element);
    308         ElementVector* pe3=CreatePVectorShelf(element);
    309         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    310 
    311         /*clean-up and return*/
    312         delete pe1;
    313         delete pe2;
    314         delete pe3;
    315         return pe;
    316 }/*}}}*/
    317 ElementVector* AgeAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
    318 
    319         /* Check if ice in element */
    320         if(!element->IsIceInElement()) return NULL;
    321 
    322         /* Geothermal flux on ice sheet base and basal friction */
    323         if(!element->IsOnBase() || element->IsAllFloating()) return NULL;
    324 
    325         IssmDouble  dt,Jdet,geothermalflux,vx,vy,vz;
    326         IssmDouble  alpha2,scalar,basalfriction,heatflux;
    327         IssmDouble *xyz_list_base = NULL;
    328 
    329         /*Fetch number of nodes for this finite element*/
    330         int numnodes = element->GetNumberOfNodes();
    331 
    332         /*Initialize vectors*/
    333         ElementVector* pe    = element->NewElementVector();
    334         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    335 
    336         /*Retrieve all inputs and parameters*/
    337         element->GetVerticesCoordinatesBase(&xyz_list_base);
    338         element->FindParam(&dt,TimesteppingTimeStepEnum);
    339         Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
    340         Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
    341         Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
    342         Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    343         IssmDouble  rho_ice             = element->FindParam(MaterialsRhoIceEnum);
    344         IssmDouble  heatcapacity        = element->FindParam(MaterialsHeatcapacityEnum);
    345 
    346         /*Build friction element, needed later: */
    347         Friction* friction=new Friction(element,3);
    348 
    349         /* Start  looping on the number of gaussian points: */
    350         Gauss* gauss   = element->NewGaussBase(4);
    351         while(gauss->next()){
    352 
    353                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    354                 element->NodalFunctions(basis,gauss);
    355 
    356                 geothermalflux_input->GetInputValue(&geothermalflux,gauss);
    357                 friction->GetAlpha2(&alpha2,gauss);
    358                 vx_input->GetInputValue(&vx,gauss);
    359                 vy_input->GetInputValue(&vy,gauss);
    360                 vz_input->GetInputValue(&vz,gauss);
    361                 vz = 0.;//FIXME
    362                 basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
    363                 heatflux      = (basalfriction+geothermalflux)/(rho_ice*heatcapacity);
    364 
    365                 scalar = gauss->weight*Jdet*heatflux;
    366                 if(dt!=0.) scalar=dt*scalar;
    367 
    368                 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    369         }
    370 
    371         /*Clean up and return*/
    372         delete gauss;
    373         delete friction;
    374         xDelete<IssmDouble>(basis);
    375         xDelete<IssmDouble>(xyz_list_base);
    376         return pe;
    377 }/*}}}*/
    378 ElementVector* AgeAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
    379 
    380         /* Check if ice in element */
    381         if(!element->IsIceInElement()) return NULL;
    382 
    383         IssmDouble  t_pmp,dt,Jdet,scalar_ocean,pressure;
    384         IssmDouble *xyz_list_base = NULL;
    385 
    386         /*Get basal element*/
    387         if(!element->IsOnBase() || !element->IsAllFloating()) return NULL;
    388 
    389         /*Fetch number of nodes for this finite element*/
    390         int numnodes = element->GetNumberOfNodes();
    391 
    392         /*Initialize vectors*/
    393         ElementVector* pe    = element->NewElementVector();
    394         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    395 
    396         /*Retrieve all inputs and parameters*/
    397         element->GetVerticesCoordinatesBase(&xyz_list_base);
    398         element->FindParam(&dt,TimesteppingTimeStepEnum);
    399         Input*      pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
    400         IssmDouble  gravity             = element->FindParam(ConstantsGEnum);
    401         IssmDouble  rho_water           = element->FindParam(MaterialsRhoSeawaterEnum);
    402         IssmDouble  rho_ice             = element->FindParam(MaterialsRhoIceEnum);
    403         IssmDouble  heatcapacity        = element->FindParam(MaterialsHeatcapacityEnum);
    404         IssmDouble  mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
    405 
    406         /* Start  looping on the number of gaussian points: */
    407         Gauss* gauss=element->NewGaussBase(4);
    408         while(gauss->next()){
    409 
    410                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    411                 element->NodalFunctions(basis,gauss);
    412 
    413                 pressure_input->GetInputValue(&pressure,gauss);
    414                 t_pmp=element->TMeltingPoint(pressure);
    415 
    416                 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*(t_pmp)/(heatcapacity*rho_ice);
    417                 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
    418 
    419                 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
    420         }
    421 
    422         /*Clean up and return*/
    423         delete gauss;
    424         xDelete<IssmDouble>(basis);
    425         xDelete<IssmDouble>(xyz_list_base);
    426         return pe;
    427 }/*}}}*/
    428 ElementVector* AgeAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
    429 
    430         /* Check if ice in element */
    431         if(!element->IsIceInElement()) return NULL;
    432 
    433241        /*Intermediaries*/
    434242        int         stabilization;
    435         IssmDouble  Jdet,phi,dt;
     243        IssmDouble  Jdet,dt;
    436244        IssmDouble  temperature;
    437245        IssmDouble  tau_parameter,diameter,hx,hy,hz;
    … …  
    451259        /*Retrieve all inputs and parameters*/
    452260        element->GetVerticesCoordinates(&xyz_list);
    453         IssmDouble  rho_ice             = element->FindParam(MaterialsRhoIceEnum);
    454         IssmDouble  heatcapacity        = element->FindParam(MaterialsHeatcapacityEnum);
    455         IssmDouble  thermalconductivity = 1.;
    456         IssmDouble  kappa = thermalconductivity/(rho_ice*heatcapacity);
    457261        element->FindParam(&dt,TimesteppingTimeStepEnum);
    458262        element->FindParam(&stabilization,AgeStabilizationEnum);
    … …  
    469273                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    470274                element->NodalFunctions(basis,gauss);
    471                 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
    472 
    473                 scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
     275
     276                scalar_def=1.*Jdet*gauss->weight;
    474277                if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
    475278
    … …  
    490293                        vz_input->GetInputValue(&w,gauss);
    491294
    492                         tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa);
     295                        tau_parameter=element->StabilizationParameter(u,v,w,diameter,1.e-15); //assume very small conductivity to get tau
    493296
    494297                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
    … …  
    504307                        vy_input->GetInputValue(&v,gauss);
    505308                        vz_input->GetInputValue(&w,gauss);
    506                         element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa);
     309                        element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,1.e-15); //assume very small conductivity to get tau
    507310                        tau_parameter_hor=tau_parameter_anisotropic[0];
    508311                        tau_parameter_ver=tau_parameter_anisotropic[1];
    … …  
    518321        delete gauss;
    519322        return pe;
    520 
    521323}/*}}}*/
    522324void           AgeAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    533335void           AgeAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    534336        SetActiveNodesLSMx(femmodel);
    535 
    536337        _error_("Should also automatically constrain surface/basal nodes where we have inflow");
    537 }/*}}}*/
     338
     339        /*Constrain all nodes that are grounded and unconstrain the ones that float*/
     340        for(Object* & object : femmodel->elements->objects){
     341                Element    *element  = xDynamicCast<Element*>(object);
     342
     343                if(element->IsOnSurface()){
     344                        element = element->SpawnTopElement();
     345                        int         numnodes  = element->GetNumberOfNodes();
     346                        IssmDouble *mask      = xNew<IssmDouble>(numnodes);
     347                        IssmDouble *bed       = xNew<IssmDouble>(numnodes);
     348                        IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
     349
     350                        element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
     351                        element->GetInputListOnNodes(&bed[0],BaseEnum);
     352                        element->GetInputListOnNodes(&ls_active[0],IceMaskNodeActivationEnum);
     353
     354                        for(int in=0;in<numnodes;in++){
     355                                Node* node=element->GetNode(in);
     356                                if(mask[in]<0. && ls_active[in]==1.){
     357                                        node->Activate();
     358                                }
     359                                else{
     360                                        node->Deactivate();
     361                                        node->ApplyConstraint(0,bed[in]);
     362                                }
     363                        }
     364                        xDelete<IssmDouble>(mask);
     365                        xDelete<IssmDouble>(bed);
     366                        xDelete<IssmDouble>(ls_active);
     367                }
     368                else if(element->IsOnBase()){
     369                        element = element->SpawnBasalElement();
     370                        _error_("not supported");
     371                }
     372        }
     373}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/AgeAnalysis.h

    r27154 r27572  
    2626                ElementMatrix* CreateJacobianMatrix(Element* element);
    2727                ElementMatrix* CreateKMatrix(Element* element);
    28                 ElementMatrix* CreateKMatrixShelf(Element* element);
    29                 ElementMatrix* CreateKMatrixVolume(Element* element);
    3028                ElementVector* CreatePVector(Element* element);
    31                 ElementVector* CreatePVectorSheet(Element* element);
    32                 ElementVector* CreatePVectorShelf(Element* element);
    33                 ElementVector* CreatePVectorVolume(Element* element);
    3429                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3530                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
Note: See TracChangeset for help on using the changeset viewer.