Changeset 23900


Ignore:
Timestamp:
04/27/19 10:01:12 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added mass flux at the ice front vs where the levelset hits 0

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

Legend:

Unmodified
Added
Removed

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

    r23888 r23900  
    241241                virtual IssmDouble IceVolumeAboveFloatation(bool scaled)=0;
    242242                virtual IssmDouble IcefrontMassFlux(bool scaled){_error_("not implemented");};
     243                virtual IssmDouble IcefrontMassFluxLevelset(bool scaled){_error_("not implemented");};
    243244                virtual IssmDouble GroundinglineMassFlux(bool scaled){_error_("not implemented");};
    244245                virtual void       InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;

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

    r23893 r23900  
    15071507
    15081508        /* Intermediaries */
    1509         int i, dir,nrfrontnodes;
    15101509        IssmDouble  levelset[NUMVERTICES];
    15111510        int         indicesfront[NUMVERTICES];
    … …  
    15151514
    15161515        /* Get nodes where there is no ice */
    1517         nrfrontnodes=0;
    1518         for(i=0;i<NUMVERTICES;i++){
     1516        int num_frontnodes=0;
     1517        for(int i=0;i<NUMVERTICES;i++){
    15191518                if(levelset[i]>=0.){
    1520                         indicesfront[nrfrontnodes]=i;
    1521                         nrfrontnodes++;
    1522                 }
    1523         }
    1524 
    1525         _assert_(nrfrontnodes==2);
     1519                        indicesfront[num_frontnodes]=i;
     1520                        num_frontnodes++;
     1521                }
     1522        }
     1523        _assert_(num_frontnodes==2);
    15261524
    15271525        /* arrange order of frontnodes such that they are oriented counterclockwise */
    … …  
    15321530        }       
    15331531
    1534         IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
     1532        IssmDouble* xyz_front = xNew<IssmDouble>(3*2);
    15351533        /* Return nodes */
    1536         for(i=0;i<nrfrontnodes;i++){
    1537                 for(dir=0;dir<3;dir++){
    1538                         xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
     1534        for(int i=0;i<2;i++){
     1535                for(int j=0;j<3;j++){
     1536                        xyz_front[3*i+j]=xyz_list[3*indicesfront[i]+j];
    15391537                }
    15401538        }
    … …  
    19531951/*}}}*/
    19541952IssmDouble Tria::IcefrontMassFlux(bool scaled){/*{{{*/
     1953
     1954        /*Make sure there is an ice front here*/
     1955        if(!IsIceInElement() || !IsIcefront()) return 0;
     1956
     1957        /*Scaled not implemented yet...*/
     1958        _assert_(!scaled);
     1959
     1960        /*Get domain type*/
     1961        int domaintype;
     1962        parameters->FindParam(&domaintype,DomainTypeEnum);
     1963
     1964        /*Get ice front coordinates*/
     1965        IssmDouble *xyz_list = NULL;
     1966        IssmDouble* xyz_front = NULL;
     1967        this->GetVerticesCoordinates(&xyz_list);
     1968        this->GetIcefrontCoordinates(&xyz_front,xyz_list,MaskIceLevelsetEnum);
     1969
     1970        /*Get normal vector*/
     1971        IssmDouble normal[3];
     1972        this->NormalSection(&normal[0],xyz_front);
     1973        normal[0] = -normal[0];
     1974        normal[1] = -normal[1];
     1975
     1976        /*Get inputs*/
     1977        IssmDouble flux = 0.;
     1978        IssmDouble vx,vy,thickness,Jdet;
     1979        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     1980        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1981        Input* vx_input=NULL;
     1982        Input* vy_input=NULL;
     1983        if(domaintype==Domain2DhorizontalEnum){
     1984                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     1985                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     1986        }
     1987        else{
     1988                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     1989                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     1990        }
     1991
     1992        /*Start looping on Gaussian points*/
     1993        Gauss* gauss=this->NewGauss(xyz_list,xyz_front,3);
     1994        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1995
     1996                gauss->GaussPoint(ig);
     1997                thickness_input->GetInputValue(&thickness,gauss);
     1998                vx_input->GetInputValue(&vx,gauss);
     1999                vy_input->GetInputValue(&vy,gauss);
     2000                this->JacobianDeterminantSurface(&Jdet,xyz_front,gauss);
     2001
     2002                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
     2003        }
     2004
     2005        return flux;
     2006}
     2007/*}}}*/
     2008IssmDouble Tria::IcefrontMassFluxLevelset(bool scaled){/*{{{*/
    19552009
    19562010        /*Make sure there is an ice front here*/

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

    r23888 r23900  
    9696                IssmDouble  IceVolumeAboveFloatation(bool scaled);
    9797                IssmDouble  IcefrontMassFlux(bool scaled);
     98                IssmDouble  IcefrontMassFluxLevelset(bool scaled);
    9899                IssmDouble  GroundinglineMassFlux(bool scaled);
    99100                void        Ismip6FloatingiceMeltingRate(void);

  • TabularUnified issm/trunk-jpl/src/c/classes/FemModel.cpp

    r23888 r23900  
    15461546
    15471547}/*}}}*/
     1548void FemModel::IcefrontMassFluxLevelsetx(IssmDouble* pM, bool scaled){/*{{{*/
     1549
     1550        IssmDouble local_mass_flux = 0;
     1551        IssmDouble total_mass_flux;
     1552
     1553        for(int i=0;i<this->elements->Size();i++){
     1554                Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
     1555                local_mass_flux+=element->IcefrontMassFluxLevelset(scaled);
     1556        }
     1557        ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1558        ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1559
     1560        /*Assign output pointers: */
     1561        *pM=total_mass_flux;
     1562
     1563}/*}}}*/
    15481564void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
    15491565
    … …  
    22152231                                        case IceMassEnum:                        this->IceMassx(&double_result,false);                  break;
    22162232                                        case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(&double_result,false);         break;
     2233                                        case IcefrontMassFluxLevelsetEnum:       this->IcefrontMassFluxLevelsetx(&double_result,false);         break;
    22172234                                        case IceMassScaledEnum:                  this->IceMassx(&double_result,true);                   break;
    22182235                                        case IceVolumeEnum:                      this->IceVolumex(&double_result,false);                break;
    … …  
    24212438                case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break;
    24222439                case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, false); break;
     2440                case IcefrontMassFluxLevelsetEnum:       this->IcefrontMassFluxLevelsetx(responses, false); break;
    24232441                case GroundedAreaEnum:                   this->GroundedAreax(responses, false); break;
    24242442                case GroundedAreaScaledEnum:             this->GroundedAreax(responses, true); break;

  • TabularUnified issm/trunk-jpl/src/c/classes/FemModel.h

    r23888 r23900  
    102102                void IcefrontAreax();
    103103                void IcefrontMassFluxx(IssmDouble* presponse, bool scaled);
     104                void IcefrontMassFluxLevelsetx(IssmDouble* presponse, bool scaled);
    104105                void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled);
    105106                void IceMassx(IssmDouble* pV, bool scaled);

  • TabularUnified issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r23897 r23900  
    10561056        GroundinglineMassFluxEnum,
    10571057        IcefrontMassFluxEnum,
     1058        IcefrontMassFluxLevelsetEnum,
    10581059        MasstransportAnalysisEnum,
    10591060        MasstransportSolutionEnum,

  • TabularUnified issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r23897 r23900  
    10601060                case GroundinglineMassFluxEnum : return "GroundinglineMassFlux";
    10611061                case IcefrontMassFluxEnum : return "IcefrontMassFlux";
     1062                case IcefrontMassFluxLevelsetEnum : return "IcefrontMassFluxLevelset";
    10621063                case MasstransportAnalysisEnum : return "MasstransportAnalysis";
    10631064                case MasstransportSolutionEnum : return "MasstransportSolution";

  • TabularUnified issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r23897 r23900  
    10841084              else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
    10851085              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
     1086              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    10861087              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    10871088              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
    … …  
    11201121              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    11211122              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    1122               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
     1126              if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     1127              else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
    11271128              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
    11281129              else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    … …  
    12431244              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    12441245              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    1245               else if (strcmp(name,"Water")==0) return WaterEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
     1249              if (strcmp(name,"Water")==0) return WaterEnum;
     1250              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    12501251              else if (strcmp(name,"XY")==0) return XYEnum;
    12511252              else if (strcmp(name,"XYZ")==0) return XYZEnum;
Note: See TracChangeset for help on using the changeset viewer.