Ignore:
Timestamp:
09/11/14 07:07:06 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on sea ice model

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18502 r18504  
    1010}/*}}}*/
    1111void SeaiceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     12        parameters->AddObject(iomodel->CopyConstantObject(SeaiceMinConcentrationEnum));
     13        parameters->AddObject(iomodel->CopyConstantObject(SeaiceMinThicknessEnum));
     14        parameters->AddObject(iomodel->CopyConstantObject(SeaiceMaxThicknessEnum));
    1215}/*}}}*/
    1316void SeaiceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    … …  
    3437        iomodel->FetchDataToInput(elements,SurfaceforcingsWindVxEnum);
    3538        iomodel->FetchDataToInput(elements,SurfaceforcingsWindVyEnum);
    36         iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
     39        iomodel->FetchDataToInput(elements,DamageEnum);
    3740        iomodel->FetchDataToInput(elements,VxStarEnum);
    3841        iomodel->FetchDataToInput(elements,VyStarEnum);
    … …  
    4043        iomodel->FetchDataToInput(elements,VyEnum);
    4144        iomodel->FetchDataToInput(elements,SeaiceCoriolisFactorEnum);
     45        iomodel->FetchDataToInput(elements,StressTensorPredictorxxEnum);
     46        iomodel->FetchDataToInput(elements,StressTensorPredictoryyEnum);
     47        iomodel->FetchDataToInput(elements,StressTensorPredictorxyEnum);
     48        iomodel->FetchDataToInput(elements,MeshXEnum);
     49        iomodel->FetchDataToInput(elements,MeshYEnum);
    4250}/*}}}*/
    4351void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    … …  
    7179        /* Check if there is ice in this element */
    7280        Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
    73         if(concentration_input->Max()==0.) return NULL;
     81        IssmDouble c_min; element->FindParam(&c_min,SeaiceMinConcentrationEnum);
     82        if(concentration_input->Max()<=c_min) return NULL;
    7483
    7584        /*Intermediaries */
    … …  
    96105        element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
    97106        element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
    98         IssmDouble rho_i           = element->GetMaterialParameter(MaterialsRhoIceEnum);
     107        IssmDouble rho_i                  = element->GetMaterialParameter(MaterialsRhoIceEnum);
     108        IssmDouble time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
    99109        Input*     thickness_input = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
    100110        Input*     vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
    … …  
    103113        Input*     oceanvy_input   = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
    104114
     115        /*Get minimum inertia to avoid 0 in the line, and time_ratio*/
     116        IssmDouble minimum_inertia = rho_i*0.01;
     117        IssmDouble time_ratio=(1-dt/time_relaxation_stress);
     118
    105119        /* Start  looping on the number of gaussian points: */
    106120        Gauss* gauss=element->NewGauss(2);
    … …  
    116130                /*Create first part of the stiffness matrix*/
    117131                this->GetM(M,element,xyz_list,gauss);
    118                 D_scalar=rho_i*thickness*gauss->weight*Jdet;
     132                D_scalar=(rho_i*thickness+minimum_inertia)/dt*gauss->weight*Jdet;
    119133                TripleMultiply(M,1,numdof,1,
    120134                                        &D_scalar,1,1,0,
    … …  
    122136                                        &Ke->values[0],1);
    123137
    124                 /*Create Elasitc part*/
     138                /*Create Elastic part*/
    125139                this->GetB(B,element,xyz_list,gauss);
    126140                this->CreateCTensor(&C[0][0],element,gauss);
    127141                for(int i=0;i<3;i++){
    128142                        for(int j=0;j<3;j++){
    129                                 C[i][j] = gauss->weight*Jdet*dt*C[i][j];
     143                                C[i][j] = dt*time_ratio*gauss->weight*Jdet*C[i][j];
    130144                        }
    131145                }
    … …  
    141155                vy_input->GetInputValue(&vy,gauss);
    142156                vnorm = sqrt(pow(oceanvx-vx,2) + pow(oceanvy-vy,2));
    143                 D_scalar = dt*concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
     157                D_scalar = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm)*cos(ocean_turning_angle)*gauss->weight*Jdet;
    144158                TripleMultiply(M,1,numdof,1,
    145159                                        &D_scalar,1,1,0,
    … …  
    160174        /* Check if there is ice in this element */
    161175        Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
    162         if(concentration_input->Max()==0.) return NULL;
     176        IssmDouble c_min; element->FindParam(&c_min,SeaiceMinConcentrationEnum);
     177        if(concentration_input->Max()<=c_min) return NULL;
    163178
    164179        /*Intermediaries */
    165         IssmDouble  air_coef,ocean_coef,constant_part;
    166         IssmDouble  rho_ice,rho_air,rho_ocean,gravity;
     180        IssmDouble  air_coef,ocean_coef,constant_part,time_relaxation_stress;
     181        IssmDouble  rho_ice,rho_air,rho_ocean,gravity,omega;
    167182        IssmDouble  vx,vy,vxstar,vystar,windvx,windvy,oceanvx,oceanvy,vnorm;
    168183        IssmDouble  air_lin_drag_coef,air_quad_drag_coef;
    169184        IssmDouble  ocean_lin_drag_coef,ocean_quad_drag_coef;
    170         IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet;
     185        IssmDouble  concentration,thickness,coriolis_factor,dt,Jdet,D_scalar;
     186        IssmDouble  sigma_xx,sigma_yy,sigma_xy,sigma_vec[3];
    171187        IssmDouble  ocean_turning_angle,dssh[2];
    172188        IssmDouble* xyz_list = NULL;
    … …  
    177193
    178194        /*Initialize Element vector and vectors*/
    179         ElementVector* pe    = element->NewElementVector();
    180         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     195        ElementVector* pe      = element->NewElementVector();
     196        IssmDouble*    B       = xNew<IssmDouble>(3*numdof);
     197        IssmDouble*    basis   = xNew<IssmDouble>(numnodes);
     198        IssmDouble     D[3][3] = {0.};
    181199
    182200        /*Retrieve all inputs and parameters*/
    … …  
    192210        element->FindParam(&ocean_quad_drag_coef,BasalforcingsOceanQuadDragCoefEnum);
    193211        element->FindParam(&ocean_turning_angle,BasalforcingsOceanTurningAngleEnum);
    194         rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
    195         gravity = element->GetMaterialParameter(ConstantsGEnum);
     212        rho_ice                = element->GetMaterialParameter(MaterialsRhoIceEnum);
     213        gravity                = element->GetMaterialParameter(ConstantsGEnum);
     214        omega                  = element->GetMaterialParameter(ConstantsOmegaEnum);
     215        time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
    196216        Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
    197217        Input* coriolis_fact_input = element->GetInput(SeaiceCoriolisFactorEnum);  _assert_(coriolis_fact_input);
    … …  
    205225        Input* oceanvy_input       = element->GetInput(BasalforcingsOceanVyEnum);  _assert_(oceanvy_input);
    206226        Input* oceanssh_input      = element->GetInput(BasalforcingsOceanSshEnum); _assert_(oceanssh_input);
     227        Input* sigma_xx_input      = element->GetInput(StressTensorxxEnum);        _assert_(sigma_xx_input);
     228        Input* sigma_yy_input      = element->GetInput(StressTensoryyEnum);        _assert_(sigma_yy_input);
     229        Input* sigma_xy_input      = element->GetInput(StressTensorxyEnum);        _assert_(sigma_xy_input);
     230
     231        /*Calculate time_ratio*/
     232        IssmDouble time_ratio=(1-dt/time_relaxation_stress);
    207233
    208234        /* Start  looping on the number of gaussian points: */
    … …  
    214240                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    215241                element->NodalFunctions(basis, gauss);
     242                this->GetB(B,element,xyz_list,gauss);
    216243
    217244                /*Get all inputs on gauss point*/
    … …  
    227254                oceanvx_input->GetInputValue(&oceanvx,gauss);
    228255                oceanvy_input->GetInputValue(&oceanvy,gauss);
     256                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
     257                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
     258                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
    229259                oceanssh_input->GetInputDerivativeValue(&dssh[0],xyz_list,gauss);
    230260
    … …  
    239269                constant_part = concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*vnorm);
    240270                for(int i=0;i<numnodes;i++){
    241                         pe->values[i*2+0] += dt*constant_part*windvx*Jdet*gauss->weight*basis[i];
    242                         pe->values[i*2+1] += dt*constant_part*windvy*Jdet*gauss->weight*basis[i];
     271                        pe->values[i*2+0] += constant_part*windvx*Jdet*gauss->weight*basis[i];
     272                        pe->values[i*2+1] += constant_part*windvy*Jdet*gauss->weight*basis[i];
    243273                }
    244274
    … …  
    247277                constant_part = concentration*ocean_coef*rho_ocean*(ocean_lin_drag_coef+ocean_quad_drag_coef*vnorm);
    248278                for(int i=0;i<numnodes;i++){
    249                         pe->values[i*2+0] += dt*constant_part*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
    250                         pe->values[i*2+1] += dt*constant_part*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
    251 
    252                         pe->values[i*2+0] += dt*constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
    253                         pe->values[i*2+1] += dt*constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
     279                        pe->values[i*2+0] += constant_part*omega*(oceanvy-vy)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
     280                        pe->values[i*2+1] += constant_part*omega*(vx-oceanvx)*sin(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
     281
     282                        pe->values[i*2+0] += constant_part*oceanvx*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
     283                        pe->values[i*2+1] += constant_part*oceanvy*cos(ocean_turning_angle)*Jdet*gauss->weight*basis[i];
     284
     285                        pe->values[i*2+0] += (-rho_ice*thickness*gravity*omega*dssh[0])*Jdet*gauss->weight*basis[i];
     286                        pe->values[i*2+1] += (-rho_ice*thickness*gravity*omega*dssh[1])*Jdet*gauss->weight*basis[i];
    254287                }
    255288
    256289                /*Coriolis forces (use ustar)*/
    257290                for(int i=0;i<numnodes;i++){
    258                         pe->values[i*2+0] += dt*(-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
    259                         pe->values[i*2+1] += dt*(+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
    260 
    261                         pe->values[i*2+0] += dt*(-rho_ice*thickness*gravity*dssh[0])*Jdet*gauss->weight*basis[i];
    262                         pe->values[i*2+1] += dt*(-rho_ice*thickness*gravity*dssh[1])*Jdet*gauss->weight*basis[i];
    263                 }
     291                        pe->values[i*2+0] += (-rho_ice*thickness*coriolis_factor*vystar)*Jdet*gauss->weight*basis[i];
     292                        pe->values[i*2+1] += (+rho_ice*thickness*coriolis_factor*vxstar)*Jdet*gauss->weight*basis[i];
     293                }
     294               
     295                /*Add elastic part of previous time step*/
     296                sigma_vec[0] = sigma_xx;
     297                sigma_vec[1] = sigma_yy;
     298                sigma_vec[2] = sigma_xy;
     299                D[0][0] = D[1][1] = D[2][2] = time_ratio*thickness*Jdet*gauss->weight;
     300                TripleMultiply(B,3,numdof,1,
     301                                        &D[0][0],3,3,0,
     302                                        &sigma_vec[0],3,1,0,
     303                                        &pe->values[0],1);
    264304        }
    265305
    … …  
    267307        xDelete<IssmDouble>(xyz_list);
    268308        xDelete<IssmDouble>(basis);
     309        xDelete<IssmDouble>(B);
    269310        delete gauss;
    270311        return pe;
    … …  
    323364
    324365        /*Intermediaries*/
     366        IssmDouble          c_min;
    325367        Vector<IssmDouble>* mask        = NULL;
    326368        IssmDouble*         serial_mask = NULL;
     369
     370        /*Get minimum concentration allowed*/
     371        femmodel->parameters->FindParam(&c_min,SeaiceMinConcentrationEnum);
    327372
    328373        /*Step 1: update mask of active nodes*/
    … …  
    334379                int    numnodes            = element->GetNumberOfNodes();
    335380                Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
    336                 if(concentration_input->Max()>0.){
     381                if(concentration_input->Max()>c_min){
    337382                        for(int in=0;in<numnodes;in++) mask->SetValue(element->nodes[in]->Sid(),1.,INS_VAL);
    338383                }
    … …  
    379424
    380425        /*Get damage input at this location*/
    381         Input* damage_input        = element->GetInput(MaterialsDamageEnum);     _assert_(damage_input);
     426        Input* damage_input        = element->GetInput(DamageEnum);              _assert_(damage_input);
    382427        Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
    383428        Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
    … …  
    518563
    519564                /*Get current stress state*/
    520                 Input*   sigma_xx_input = element->GetInput(StressTensorxxEnum);
    521                 Input*   sigma_yy_input = element->GetInput(StressTensoryyEnum);
    522                 Input*   sigma_xy_input = element->GetInput(StressTensorxyEnum);
    523                 if(sigma_xx_input) sigma_xx_input->GetInputAverage(&sigma_xx);
    524                 else               sigma_xx = 0.;
    525                 if(sigma_yy_input) sigma_yy_input->GetInputAverage(&sigma_yy);
    526                 else               sigma_yy = 0.;
    527                 if(sigma_xy_input) sigma_xy_input->GetInputAverage(&sigma_xy);
    528                 else               sigma_xy = 0.;
     565                Input*   sigma_xx_input = element->GetInput(StressTensorPredictorxxEnum); _assert_(sigma_xx_input);
     566                Input*   sigma_yy_input = element->GetInput(StressTensorPredictoryyEnum); _assert_(sigma_yy_input);
     567                Input*   sigma_xy_input = element->GetInput(StressTensorPredictorxyEnum); _assert_(sigma_xy_input);
     568                sigma_xx_input->GetInputAverage(&sigma_xx);
     569                sigma_yy_input->GetInputAverage(&sigma_yy);
     570                sigma_xy_input->GetInputAverage(&sigma_xy);
    529571
    530572                /* Compute the invariants of the elastic deformation and instantaneous deformation rate */
    … …  
    537579
    538580                /* estimate the internal constraints using the current elastic deformation */
    539                 Input* damage_input    = element->GetInput(MaterialsDamageEnum); _assert_(damage_input);
    540                 Input* damagenew_input = element->GetInput(MaterialsDamageEnum); _assert_(damagenew_input);//FIXME
     581                Input* damage_input    = element->GetInput(DamageEnum); _assert_(damage_input);
    541582                damage_input->GetInputAverage(&damage);
    542                 damagenew_input->GetInputAverage(&damage_new);
    543583                damage_test = damage;
     584                damage_new  = damage;
    544585                if(sigma_n>traction || sigma_n<-compression_max){
    545586                        if(sigma_n>traction){
    … …  
    564605                /* The damage variable is changed */
    565606                damage_new=damage_test;
    566                 element->AddInput(MaterialsDamageEnum,&damage_test,P0Enum);//FIXME
     607                element->AddInput(DamageEnum,&damage_new,P0Enum);
    567608
    568609                /* Recompute the internal stress*/
Note: See TracChangeset for help on using the changeset viewer.