Changeset 26653


Ignore:
Timestamp:
11/22/21 08:51:54 (3 years ago)
Author:
vverjans
Message:

CHG: correction of weights computation in Tria::GetFractionGeometry and Penta::GetFractionGeometry2D, TotalSmb can account for partially ice-covered elements (feature currently disabled with if(false) condition)

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r26626 r26653  
    12961296       
    12971297        /*Compute weights*/
    1298         gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
     1298        gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2);
    12991299
    13001300        total_weight = 0.0;
    … …  
    13061306        }
    13071307
    1308         /*Normalize to phi*/
    1309         if(total_weight>0.){
    1310                 for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
    1311         }
     1308        /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
     1309   if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++) weights[i] = weights[i]*phi/total_weight;
    13121310        else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
    13131311
    … …  
    22072205                IssmDouble basetot;
    22082206                height = 0.0;
    2209                 for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
     2207                for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]/phi*heights[i];
    22102208                basetot = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    22112209                base    = basetot*phi; 
    … …  
    22172215                        /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
    22182216                        scalefactor = 0.0;
    2219                         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
     2217                        for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    22202218                        base = base*scalefactor;
    22212219                        xDelete<IssmDouble>(scalefactor_vertices);
    … …  
    43894387        IssmDouble base,smb,rho_ice,scalefactor;
    43904388        IssmDouble Total_Smb=0;
     4389        IssmDouble lsf[NUMVERTICES];
    43914390        IssmDouble xyz_list[NUMVERTICES][3];
    43924391
    … …  
    44044403
    44054404        /*Now get the average SMB over the element*/
    4406         Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    4407 
    4408         smb_input->GetInputAverage(&smb);
    4409         if(scaled==true){
    4410                 Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    4411                 scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
    4412         }
     4405        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     4406   if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
     4407                /*Partially ice-covered element*/
     4408                bool mainlyice;
     4409      int point;
     4410      IssmDouble* smb_vertices = xNew<IssmDouble>(NUMVERTICES);
     4411                IssmDouble  weights[NUMVERTICES2D];
     4412                IssmDouble  lsf2d[NUMVERTICES2D];
     4413      IssmDouble f1,f2,phi;
     4414      Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
     4415                for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
     4416                GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
     4417                smb = 0.0;
     4418                for(int i=0;i<NUMVERTICES2D;i++) smb += weights[i]*smb_vertices[i];
     4419
     4420                if(scaled==true){
     4421         IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
     4422         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     4423         /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
     4424         scalefactor = 0.0;
     4425         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     4426         xDelete<IssmDouble>(scalefactor_vertices);
     4427                }
     4428                else scalefactor = 1.0;
     4429
     4430                /*Cleanup*/
     4431      xDelete<IssmDouble>(smb_vertices);
     4432        }
     4433
    44134434        else{
    4414                 scalefactor=1.;
    4415         }
     4435                /*Fully ice-covered element*/
     4436                Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
     4437                smb_input->GetInputAverage(&smb);
     4438
     4439                if(scaled==true){
     4440                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     4441                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     4442                }
     4443                else scalefactor=1.0;
     4444        }
     4445
    44164446        Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
    44174447

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r26626 r26653  
    18171817
    18181818        /*Compute weights:*/
    1819         gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
     1819        gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2); //VV correction (16Nov2021)
    18201820
    18211821        total_weight=0;
    … …  
    18261826                total_weight+=gauss->weight;
    18271827        }
    1828         //normalize to phi.
     1828        /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
    18291829        if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi;
    18301830        else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
    … …  
    35463546        }*/
    35473547
    3548    IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
    3549    Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     3548   IssmDouble lsf[NUMVERTICES];
     3549        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    35503550   /*Deal with partially ice-covered elements*/
    35513551        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
    … …  
    35633563      for(int i=0;i<NUMVERTICES;i++) Hice[i] = surfaces[i]-bases[i];
    35643564      Haverage = 0.0;
    3565       for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]*Hice[i];
     3565                /*Use weights[i]/phi to get average thickness over subelement*/
     3566      for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]/phi*Hice[i];
    35663567                /*Get back area of ice-covered base*/
    35673568                area_basetot = this->GetArea();
    … …  
    35733574                        Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    35743575                        scalefactor = 0.0;
    3575                         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
     3576                        for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    35763577                        area_base = area_base*scalefactor;
    35773578                        xDelete<IssmDouble>(scalefactor_vertices);
    … …  
    36063607        xDelete<IssmDouble>(H);
    36073608        xDelete<IssmDouble>(SF);
    3608         xDelete<IssmDouble>(lsf);
    36093609
    36103610        if(domaintype==Domain2DverticalEnum){
    … …  
    54795479        IssmDouble base,smb,rho_ice,scalefactor;
    54805480        IssmDouble Total_Smb=0;
     5481        IssmDouble lsf[NUMVERTICES];
    54815482        IssmDouble xyz_list[NUMVERTICES][3];
    54825483
    … …  
    54925493         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    54935494        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
    5494 
     5495       
    54955496        /*Now get the average SMB over the element*/
    5496         Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    5497         smb_input->GetInputAverage(&smb);       // average smb on element in m ice s-1
    5498         if(scaled==true){
    5499                 Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    5500                 scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5497        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     5498        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
     5499                /*Partially ice-covered element*/
     5500                bool mainlyice;
     5501      int point;
     5502      IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
     5503      IssmDouble* smb_vertices  = xNew<IssmDouble>(NUMVERTICES);
     5504      IssmDouble f1,f2,phi;
     5505               
     5506                Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
     5507                GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
     5508                smb = 0.0;
     5509                for(int i=0;i<NUMVERTICES;i++) smb += weights[i]*smb_vertices[i];
     5510       
     5511                if(scaled==true){
     5512         IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
     5513         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     5514         scalefactor = 0.0;
     5515         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     5516         xDelete<IssmDouble>(scalefactor_vertices);
     5517      }
     5518                else scalefactor = 1.0;
     5519
     5520                /*Cleanup*/
     5521      xDelete<IssmDouble>(weights);
     5522      xDelete<IssmDouble>(smb_vertices);
    55015523        }
    55025524        else{
    5503                 scalefactor=1.;
    5504         }
     5525                /*Fully ice-covered element*/
     5526                Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
     5527                smb_input->GetInputAverage(&smb);   // average smb on element in m ice s-1
     5528
     5529                if(scaled==true){
     5530                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     5531                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5532                }
     5533                else scalefactor=1.0;
     5534        }
     5535       
    55055536   Total_Smb=rho_ice*base*smb*scalefactor;      // smb on element in kg s-1
    55065537
Note: See TracChangeset for help on using the changeset viewer.