Changeset 26164


Ignore:
Timestamp:
03/30/21 18:44:37 (4 years ago)
Author:
Eric.Larour
Message:

CHG: test2002 replicated between old and new code.

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r26152 r26164  
    248248                virtual void       GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
    249249                virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
     250                virtual void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl)=0;
     251                virtual IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl)=0;
    250252                virtual IssmDouble GetIcefrontArea(){_error_("not implemented");};
    251253                virtual void       GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
    … …  
    295297                virtual void       JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
    296298                virtual void       JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
     299                virtual void       LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum)=0;
    297300                virtual void       Marshall(MarshallHandle* marshallhandle)=0;
    298301                virtual IssmDouble Masscon(IssmDouble* levelset)=0;
    … …  
    394397                virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    395398                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
    396                 virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0;
    397                 virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0;
    398                 virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,  SealevelMasks* masks)=0;
    399                 virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
     399                virtual void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0;
     400                virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0;
    400401                virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
    401402                #endif

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

    r26152 r26164  
    8383                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    8484                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
     85                void           GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
     86                IssmDouble     GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
    8587                IssmDouble              GetIcefrontArea();
    8688                void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    … …  
    125127                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    126128                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
     129                void           LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");};
    127130                IssmDouble     Masscon(IssmDouble* levelset){_error_("not implemented yet");};
    128131                IssmDouble     MassFlux(IssmDouble* segment);
    … …  
    231234                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    232235                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    233                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    234                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,SealevelMasks* masks){_error_("not implemented yet");};
    235                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    236                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
     236                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
     237                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    237238                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
    238239                #endif

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r26152 r26164  
    6363                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
    6464                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
     65                void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
     66                IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
    6567                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    6668                Input*     GetInput(int enumtype);
    … …  
    9799                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    98100                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
     101                void        LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");};
    99102                IssmDouble  Masscon(IssmDouble* levelset){_error_("not implemented yet");};
    100103                IssmDouble  MassFlux(IssmDouble* segment){_error_("not implemented yet");};
    … …  
    186189                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    187190                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    188                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    189                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    190                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    191                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
     191                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
     192                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    192193                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
    193194

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r26152 r26164  
    6767                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
    6868                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
     69                void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
     70                IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
    6971                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
    7072                Input*     GetInput(int enumtype);
    … …  
    103105                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    104106                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
     107                void        LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented yet");};
    105108                IssmDouble  Masscon(IssmDouble* levelset){_error_("not implemented yet");};
    106109                IssmDouble  MassFlux(IssmDouble* segment){_error_("not implemented yet");};
    … …  
    192195                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    193196                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    194                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    195                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    196                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    197                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
     197                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
     198                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    198199                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
    199200

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

    r26158 r26164  
    16001600}
    16011601/*}}}*/
    1602 void       Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
     1602void       Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlyfloating){/*{{{*/
    16031603        /*Computeportion of the element that is grounded*/
    16041604
    … …  
    16511651        *fraction1=f1;
    16521652        *fraction2=f2;
    1653         *mainlyfloating=floating;
     1653        *pmainlyfloating=floating;
    16541654}
    16551655/*}}}*/
    … …  
    17211721        _assert_(phi<=1. && phi>=0.);
    17221722        return phi;
     1723}
     1724/*}}}*/
     1725IssmDouble Tria::GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){/*{{{*/
     1726       
     1727        /*Computeportion of the element that is on a levelset: */
     1728        bool              mainlyfloating = true;
     1729        int               domaintype,index1,index2;
     1730        const IssmPDouble epsilon        = 1.e-15;
     1731        IssmDouble        phi,s1,s2;
     1732
     1733        /*Recover parameters and values*/
     1734        parameters->FindParam(&domaintype,DomainTypeEnum);
     1735
     1736        /*Be sure that values are not zero*/
     1737        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     1738        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     1739        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     1740
     1741        if(domaintype==Domain2DverticalEnum){
     1742                this->EdgeOnBaseIndices(&index1,&index2);
     1743                if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded
     1744                else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating
     1745                else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded
     1746                        phi=1./(1.-gl[index1]/gl[index2]);
     1747                }
     1748                else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded
     1749                        phi=1./(1.-gl[index2]/gl[index1]);
     1750                }
     1751
     1752        }
     1753        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     1754                /*Check that not all nodes are grounded or floating*/
     1755                if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
     1756                        phi=1;
     1757                }
     1758                else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
     1759                        phi=0;
     1760                }
     1761                else{
     1762                        /*Figure out if two nodes are floating or grounded*/
     1763                        if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
     1764
     1765                        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1766                                s1=gl[2]/(gl[2]-gl[1]);
     1767                                s2=gl[2]/(gl[2]-gl[0]);
     1768                        }
     1769                        else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     1770                                s1=gl[0]/(gl[0]-gl[1]);
     1771                                s2=gl[0]/(gl[0]-gl[2]);
     1772                        }
     1773                        else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     1774                                s1=gl[1]/(gl[1]-gl[0]);
     1775                                s2=gl[1]/(gl[1]-gl[2]);
     1776                        }
     1777                        else _error_("case not possible");
     1778                        if(mainlyfloating){
     1779                                phi = (1-s1*s2);
     1780                        }
     1781                        else{
     1782                                phi = s1*s2;
     1783                        }
     1784                }
     1785        }
     1786        else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     1787
     1788        _assert_(phi<=1. && phi>=0.);
     1789        return 1-phi;
     1790}
     1791/*}}}*/
     1792void       Tria::GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){/*{{{*/
     1793       
     1794        /*Computeportion of the element that is grounded*/
     1795        bool               negative=false;
     1796        int                point;
     1797        const IssmPDouble  epsilon= 1.e-15;
     1798        IssmDouble         f1,f2;
     1799
     1800        /*Be sure that values are not zero*/
     1801        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     1802        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     1803        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     1804
     1805        /*Check that not all nodes are positive or negative: */
     1806        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
     1807                point=0;
     1808                f1=1.;
     1809                f2=1.;
     1810        }
     1811        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
     1812                point=0;
     1813                f1=0.;
     1814                f2=0.;
     1815        }
     1816        else{
     1817                if(gl[0]*gl[1]*gl[2]<0) negative=true;
     1818
     1819                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1820                        point=2;
     1821                        f1=gl[2]/(gl[2]-gl[0]);
     1822                        f2=gl[2]/(gl[2]-gl[1]);
     1823                }
     1824                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     1825                        point=0;
     1826                        f1=gl[0]/(gl[0]-gl[1]);
     1827                        f2=gl[0]/(gl[0]-gl[2]);
     1828                }
     1829                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     1830                        point=1;
     1831                        f1=gl[1]/(gl[1]-gl[2]);
     1832                        f2=gl[1]/(gl[1]-gl[0]);
     1833                }
     1834                else _error_("case not possible");
     1835        }
     1836        *point1=point;
     1837        *fraction1=f1;
     1838        *fraction2=f2;
     1839        *pmainlynegative=negative;
    17231840}
    17241841/*}}}*/
    … …  
    55455662                IssmDouble phi=1.0;
    55465663                /*{{{*/
    5547        
    5548                
    5549 
    5550                 if (masks->isfullyfloating[this->lid]){
    5551                         I=0;
    5552                         this->AddInput(EffectivePressureEnum,&I,P0Enum);
    5553                         return;
    5554                 }
    5555                        
    5556 
    5557        
    55585664                if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    55595665
    … …  
    56055711                                total_weight+=gauss->weight;
    56065712                        }
    5607                         I=I/total_weight;
     5713                        if(total_weight) I=I/total_weight;
    56085714                        delete gauss;
    56095715                }
    … …  
    65516657                        Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
    65526658                        Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
    6553                         Grotm3[i]= - pre* (1/6.0 - 0.5*cos(2.0*lati));
     6659                        Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati));
    65546660                }
    65556661                this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
    … …  
    66876793
    66886794        /*Deal with ice loads if we are on grounded ice:*/
    6689         if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){
     6795        if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid] && isice){
    66906796
    66916797                /*Compute fraction of the element that is grounded: {{{*/
    … …  
    67346840                                total_weight+=gauss->weight;
    67356841                        }
    6736                         I=I/total_weight;
     6842                        if(total_weight) I=I/total_weight;
    67376843                        delete gauss;
    67386844                }
    … …  
    67486854
    67496855        /*Deal with water loads if we are on ground:*/
    6750         if(!masks->isfullyfloating[this->lid]){
     6856        if(!masks->isfullyfloating[this->lid] && ishydro){
    67516857
    67526858                constant+=10; //1 for water.
    … …  
    67696875
    67706876        /*Deal with ocean bottom pressures:*/
    6771         if(masks->isoceanin[this->lid]){
     6877        if(masks->isoceanin[this->lid] && isocean){
    67726878
    67736879                constant+=1; //1 for bottom pressure.
    … …  
    67966902}
    67976903/*}}}*/
    6798 void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
     6904void       Tria::SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationvector){ /*{{{*/
    67996905
    68006906        /*sal green function:*/
    68016907        IssmDouble* G=NULL;
    6802         IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
     6908        IssmDouble SealevelGRD[NUMVERTICES]={0};
    68036909        IssmDouble oceanaverage,oceanarea=0;
    68046910        IssmDouble rho_water;
    … …  
    68286934                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    68296935               
    6830                 for(int i=0;i<NUMVERTICES;i++) {
    6831                         for (int e=0;e<nel;e++){
    6832                                 SealevelGRD[i]+=G[i*nel+e]*loads[e];
     6936                if(allsealevelloads){
     6937                        for(int i=0;i<NUMVERTICES;i++) {
     6938                                for (int e=0;e<nel;e++){
     6939                                        SealevelGRD[i]+=G[i*nel+e]*(allsealevelloads[e]+allloads[e]);
     6940                                }
     6941                        }
     6942                }
     6943                else{  //this is the initial convolution where only loads are provided
     6944                        for(int i=0;i<NUMVERTICES;i++) {
     6945                                for (int e=0;e<nel;e++){
     6946                                        SealevelGRD[i]+=G[i*nel+e]*(allloads[e]);
     6947                                }
    68336948                        }
    68346949                }
    … …  
    68366951
    68376952        /*compute ocean average over element:*/
    6838         OceanAverageOptim(&oceanaverage,&oceanarea,&SealevelGRD[0],masks);
     6953        LevelsetAverage(&oceanaverage,&oceanarea,&SealevelGRD[0],MaskOceanLevelsetEnum);
    68396954       
    68406955        /*add ocean average in the global sealevelloads vector:*/
    68416956        sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    6842         oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     6957        if(!allsealevelloads)oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    68436958       
    68446959        return;
    68456960} /*}}}*/
    6846 void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
    6847 
    6848         IssmDouble SealevelGRD[3]={0,0,0};
    6849         IssmDouble oceanaverage,oceanarea=0;
    6850         IssmDouble rho_water;
    6851         int nel;
    6852         bool sal = false;
    6853         IssmDouble* G=NULL;
    6854         int size;
    6855         IssmDouble Grotm1[3];
    6856         IssmDouble Grotm2[3];
    6857         IssmDouble Grotm3[3];
    6858         bool rotation= false;
    6859 
    6860         this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    6861         this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
    6862         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    6863         this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    6864        
    6865         if(rotation){
    6866                 Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
    6867                 Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    6868                 Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    6869                
    6870                 for(int i=0;i<NUMVERTICES;i++) SealevelGRD[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
    6871         }       
    6872 
    6873         if(sal){
    6874                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    6875 
    6876                 for(int i=0;i<NUMVERTICES;i++) {
    6877                         for (int e=0;e<nel;e++){
    6878                                 SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
    6879                         }
    6880                 }
    6881         }
    6882         OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
    6883         newsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    6884 
    6885 } /*}}}*/
    6886 void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
     6961void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){ /*{{{*/
    68876962
    68886963        IssmDouble Sealevel[3]={0,0,0};
    … …  
    69587033
    69597034} /*}}}*/
    6960 void       Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
     7035void       Tria::LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){ /*{{{*/
    69617036
    69627037        IssmDouble phi=1.0;
    6963         bool iscoastline=false;
    69647038        IssmDouble area;
    6965         IssmDouble Sg_avg=0; //output
     7039        IssmDouble average = 0;
     7040        IssmDouble total_weight=0;
     7041        int        point1;
     7042        IssmDouble fraction1,fraction2;
     7043        IssmDouble levelsetvalues[NUMVERTICES];
     7044        bool       fractiongeometryflag=true;
    69667045       
    6967         /*Do we have an ocean?:*/
    6968         if(!masks->isoceanin[this->lid]){
    6969                 *poceanarea=0;
    6970                 *poceanaverage=0;
    6971         }
    6972 
    6973         /*Do we have a coastline?:*/
    6974         if(!masks->isfullyfloating[this->lid])iscoastline=true;
    6975 
    6976         if(iscoastline){
    6977                 IssmDouble xyz_list[NUMVERTICES][3];
    6978                 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6979                 phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]);
     7046
     7047
     7048        /*retrieve value of levelset:*/
     7049        Input *levelset= this->GetInput(levelsetenum);
     7050        this->Element::GetInputListOnVertices(&levelsetvalues[0],levelsetenum);
     7051       
     7052        /*Early return if no vertices are inside the desired area:*/
     7053        if(levelset->GetInputMin()>=0){
     7054                *paverage=0;
     7055                *parea=0;
     7056                return;
    69807057        }
    69817058
    … …  
    69837060        this->Element::GetInputValue(&area,AreaEnum);
    69847061
    6985         /*Average over ocean if there is no coastline, area of the ocean
    6986          *is the areaa of the element:*/
    6987         if(!iscoastline){
    6988 
    6989                 /*Average Sg over vertices:*/
    6990                 for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[i]/NUMVERTICES;
    6991 
    6992                 *poceanaverage=Sg_avg;
    6993                 *poceanarea=area;
     7062        /*Are we fully in the desired area, in which case average over the entire area?:*/
     7063        if(levelset->GetInputMax()<=0){
     7064                for(int i=0;i<NUMVERTICES;i++) average+=field_on_localvertices[i]/NUMVERTICES;
     7065               
     7066                *parea=area;
     7067                *paverage=average;
    69947068                return;
    69957069        }
    69967070
    6997         /*Average over  the ocean only if there is a coastline. Area of the
    6998          * ocean will be the fraction times the area of the element:*/
     7071        /*What fraction of the triangle is in the desired area?:*/
     7072        IssmDouble xyz_list[NUMVERTICES][3];
     7073        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7074        phi=this->GetFractionArea(&xyz_list[0][0],levelsetvalues);
    69997075        area=phi*area;
    7000    
    7001         IssmDouble total_weight=0;
    7002         bool mainlyfloating = true;
    7003         int         point1;
    7004         IssmDouble  fraction1,fraction2;
    7005 
    7006         /*Recover portion of element that is grounded*/
    7007         this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    7008         //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part
    7009         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2);
    7010 
    7011         /* Start  looping on the number of gaussian points and average over these gaussian points: */
     7076
     7077        /*Average over  the fraction area only, using the right gaussian points: */
     7078        this->GetFractionGeometry(&point1,&fraction1,&fraction2,&fractiongeometryflag,levelsetvalues);
     7079        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,fractiongeometryflag,2);
     7080
    70127081        total_weight=0;
    7013         Sg_avg=0;
     7082        average=0;
    70147083        while(gauss->next()){
    7015                 IssmDouble Sg_gauss=0;
    7016                 TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum);
    7017                 Sg_avg+=Sg_gauss*gauss->weight;
     7084                IssmDouble field_gauss=0;
     7085                TriaRef::GetInputValue(&field_gauss, field_on_localvertices, gauss,P1Enum);
     7086                average+=field_gauss*gauss->weight;
    70187087                total_weight+=gauss->weight;
    70197088        }
    7020         if(total_weight) Sg_avg=Sg_avg/total_weight;
     7089        if(total_weight) average=average/total_weight;
     7090       
     7091        /*free ressources:*/
    70217092        delete gauss;
    70227093
    7023         *poceanaverage=Sg_avg;
    7024         *poceanarea=area;
     7094        *paverage=average;
     7095        *parea=area;
    70257096        return;
    70267097

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

    r26152 r26164  
    8484                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    8585                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
     86                void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl);
     87                IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl);
    8688                IssmDouble  GetIcefrontArea();
    8789                void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    … …  
    111113                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
    112114                bool        IsZeroLevelset(int levelset_enum);
     115                void        LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum);
    113116                IssmDouble  Masscon(IssmDouble* levelset);
    114117                IssmDouble  MassFlux(IssmDouble* segment);
    … …  
    176179                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    177180                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
    178                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
    179                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
    180                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
    181                 void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
     181                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector);
     182                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector);
    182183                void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks);
    183184                #endif

  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r26157 r26164  
    405405        SealevelMasks* masks=NULL;
    406406        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
    407         IssmDouble rotationaxismotionvector[3];
     407        IssmDouble rotationaxismotionvector[3]={0};
    408408       
    409409        Vector<IssmDouble>*    loads=NULL;
    … …  
    468468        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    469469        loads=new Vector<IssmDouble>(nel);
     470        oceanareas=new Vector<IssmDouble>(nel);
    470471        sealevelloads=new Vector<IssmDouble>(nel);
    471         oceanareas=new Vector<IssmDouble>(nel);
     472        sealevelloads->Set(0);sealevelloads->Assemble();
    472473
    473474        /*call masks core: */
    … …  
    494495        }
    495496
    496         /*convolve loads:*/
    497         for(Object* & object : femmodel->elements->objects){
    498                 Element* element = xDynamicCast<Element*>(object);
    499                 element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,rotationaxismotionvector,masks);
    500         }
    501         sealevelloads->Assemble();
    502         oceanareas->Assemble();
    503        
    504         //Get ocean area:
    505         oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
    506         if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    507 
    508         //conserve ocean mass:
    509         oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
    510         ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
    511 
    512         //broadcast sea level loads
    513         allsealevelloads=sealevelloads->ToMPISerial();
    514 
    515        
    516         //compute rotation axis motion:
    517         RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
    518        
    519         if(VerboseSolution()) _printf0_("         converging ocean GRD convolutions\n");
     497        if(VerboseSolution()) _printf0_("         converging GRD convolutions\n");
    520498        for(;;){
    521499                       
    … …  
    525503                for(Object* & object : femmodel->elements->objects){
    526504                        Element* element = xDynamicCast<Element*>(object);
    527                         element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,rotationaxismotionvector,masks);
     505                        element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector);
    528506                }
    529507                sealevelloads->Assemble();
    530                        
     508
     509                /*compute ocean areas:*/
     510                if(!allsealevelloads){ //first time in the loop
     511                        oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
     512                        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
     513                }
     514               
    531515                //Conserve ocean mass:
    532516                oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
    … …  
    544528
    545529                //early return?
    546                 if(iterations>max_nonlinear_iterations)break;
     530                if(iterations>=max_nonlinear_iterations)break;
    547531                else iterations++;
    548532        }
    … …  
    555539        for(Object* & object : femmodel->elements->objects){
    556540                Element* element = xDynamicCast<Element*>(object);
    557                 element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector,masks);
     541                element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector);
    558542        }
    559543       
    … …  
    12271211        ndS=dRSLg->Norm(NORM_TWO);
    12281212
    1229         if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
     1213        if (xIsNan<IssmDouble>(ndS)){
     1214                _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")");
     1215        }
    12301216
    12311217        if(!xIsNan<IssmDouble>(eps_rel)){
    12321218                nS=RSLg_old->Norm(NORM_TWO);
    1233                 if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
     1219                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN! (check the initial RSL)");
    12341220        }
    12351221
Note: See TracChangeset for help on using the changeset viewer.