Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 27648) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 27649) @@ -53,8 +53,17 @@ int M,N; int *segments = NULL; - iomodel->FetchData(&segments,&M,&N,"md.mesh.segments"); + if(iomodel->domaintype==Domain3DEnum){ + iomodel->FetchData(&segments,&M,&N,"md.mesh.segments2d"); + } + else if(iomodel->domaintype==Domain2DhorizontalEnum){ + iomodel->FetchData(&segments,&M,&N,"md.mesh.segments"); + } + else{ + _error_("mesh type not supported yet"); + } /*Check that the size seem right*/ _assert_(N==3); _assert_(M>=3); + for(int i=0;imy_elements[segments[i*3+2]-1]){ @@ -169,4 +178,9 @@ ElementMatrix* HydrologyShaktiAnalysis::CreateKMatrix(Element* element){/*{{{*/ + /* Check if ice in element */ + if(element->IsAllFloating() || !element->IsIceInElement()) return NULL; + if(!element->IsOnBase()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + /*Intermediaries */ IssmDouble Jdet; @@ -175,57 +189,55 @@ /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(); + ElementMatrix* Ke = basalelement->NewElementMatrix(); IssmDouble* dbasis = xNew(2*numnodes); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); + basalelement->GetVerticesCoordinates(&xyz_list); /*Get conductivity from inputs*/ - IssmDouble conductivity = GetConductivity(element); + IssmDouble conductivity = GetConductivity(basalelement); /*Get englacial storage coefficient*/ IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); - - /*Get all inputs and parameters*/ - element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum); - element->FindParam(&rho_ice,MaterialsRhoIceEnum); - element->FindParam(&g,ConstantsGEnum); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); - Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + basalelement->FindParam(&storage,HydrologyStorageEnum); + basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + + /*Get all inputs and parameters*/ + basalelement->FindParam(&rho_water,MaterialsRhoFreshwaterEnum); + basalelement->FindParam(&rho_ice,MaterialsRhoIceEnum); + basalelement->FindParam(&g,ConstantsGEnum); + Input* B_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = basalelement->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* gap_input = basalelement->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* head_input = basalelement->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(1); + Gauss* gauss=basalelement->NewGauss(1); while(gauss->next()){ - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - element->NodalFunctions(basis,gauss); - - base_input->GetInputValue(&bed,gauss); - thickness_input->GetInputValue(&thickness,gauss); - gap_input->GetInputValue(&gap,gauss); - head_input->GetInputValue(&head,gauss); - - /*Get ice A parameter*/ - B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - A=pow(B,-n); - - /*Get water and ice pressures*/ - IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); - IssmDouble pressure_water = rho_water*g*(head-bed); - if(pressure_water>pressure_ice) pressure_water = pressure_ice; - - + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); + + base_input->GetInputValue(&bed,gauss); + thickness_input->GetInputValue(&thickness,gauss); + gap_input->GetInputValue(&gap,gauss); + head_input->GetInputValue(&head,gauss); + + /*Get ice A parameter*/ + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + A=pow(B,-n); + + /*Get water and ice pressures*/ + IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.); + IssmDouble pressure_water = rho_water*g*(head-bed); + if(pressure_water>pressure_ice) pressure_water = pressure_ice; for(int i=0;ivalues[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) + gauss->weight*Jdet*storage/dt*basis[i]*basis[j] - +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j]; + +gauss->weight*Jdet*A*(n)*(pow(fabs(pressure_ice-pressure_water),(n-1))*rho_water*g)*gap*basis[i]*basis[j]; } } @@ -243,4 +255,5 @@ xDelete(dbasis); delete gauss; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; }/*}}}*/ @@ -248,5 +261,7 @@ /*Skip if water or ice shelf element*/ - if(element->IsAllFloating()) return NULL; + if(element->IsAllFloating() || !element->IsIceInElement()) return NULL; + if(!element->IsOnBase()) return NULL; + Element* basalelement = element->SpawnBasalElement(); /*Intermediaries */ @@ -259,45 +274,45 @@ /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ - ElementVector* pe = element->NewElementVector(); + ElementVector* pe = basalelement->NewElementVector(); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum); - IssmDouble g = element->FindParam(ConstantsGEnum); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); - Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); - Input* head_input = element->GetInput(HydrologyHeadEnum); _assert_(head_input); - Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); - Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); - Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); - Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); + basalelement->GetVerticesCoordinates(&xyz_list); + IssmDouble latentheat = basalelement->FindParam(MaterialsLatentheatEnum); + IssmDouble g = basalelement->FindParam(ConstantsGEnum); + IssmDouble rho_ice = basalelement->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = basalelement->FindParam(MaterialsRhoFreshwaterEnum); + Input* geothermalflux_input = basalelement->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input); + Input* head_input = basalelement->GetInput(HydrologyHeadEnum); _assert_(head_input); + Input* gap_input = basalelement->GetInput(HydrologyGapHeightEnum); _assert_(gap_input); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); + Input* B_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = basalelement->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* englacial_input = basalelement->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); + Input* lr_input = basalelement->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); + Input* br_input = basalelement->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); + Input* headold_input = basalelement->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); /*Get conductivity from inputs*/ - IssmDouble conductivity = GetConductivity(element); + IssmDouble conductivity = GetConductivity(basalelement); /*Get englacial storage coefficient*/ IssmDouble storage,dt; - element->FindParam(&storage,HydrologyStorageEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); - - /*Build friction element, needed later: */ - Friction* friction=new Friction(element,2); + basalelement->FindParam(&storage,HydrologyStorageEnum); + basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + + /*Build friction basalelement, needed later: */ + Friction* friction=new Friction(basalelement,2); /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); + Gauss* gauss=basalelement->NewGauss(2); while(gauss->next()){ - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctions(basis,gauss); + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); geothermalflux_input->GetInputValue(&G,gauss); base_input->GetInputValue(&bed,gauss); @@ -338,21 +353,21 @@ /*Compute change in sensible heat due to changes in pressure melting point*/ - dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); + dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); - meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); - - for(int i=0;ivalues[i]+=Jdet*gauss->weight* - ( - meltrate*(1/rho_water-1/rho_ice) - +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap - +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap - -beta*sqrt(vx*vx+vy*vy) - +ieb - +storage*head_old/dt - )*basis[i]; - - - } + meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); + + for(int i=0;ivalues[i]+=Jdet*gauss->weight* + ( + meltrate*(1/rho_water-1/rho_ice) + +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice + rho_water*g*bed)*gap + +(n-1)*A*pow(fabs(pressure_ice - pressure_water),n-1)*(rho_water*g*head)*gap + -beta*sqrt(vx*vx+vy*vy) + +ieb + +storage*head_old/dt + )*basis[i]; + + } + /*Clean up and return*/ xDelete(xyz_list); @@ -360,4 +375,5 @@ delete friction; delete gauss; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; return pe; }/*}}}*/ @@ -369,4 +385,7 @@ }/*}}}*/ void HydrologyShaktiAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + /*Only update if on base*/ + if(!element->IsOnBase()) return; /*Intermediary*/ @@ -411,5 +430,5 @@ /*Add input to the element: */ - element->AddInput(HydrologyHeadEnum,values,element->GetElementType()); + element->AddBasalInput(HydrologyHeadEnum,values,element->GetElementType()); /*Update reynolds number according to new solution*/ @@ -424,5 +443,5 @@ IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU; - element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum); + element->AddBasalInput(HydrologyReynoldsEnum,&reynolds,P0Enum); /*Compute new effective pressure*/ @@ -438,5 +457,48 @@ }/*}}}*/ void HydrologyShaktiAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ - /*Default, do nothing*/ + /*Update active elements based on ice levelset and ocean levelset*/ + GetMaskOfIceVerticesLSMx(femmodel,true); + SetActiveNodesLSMx(femmodel,true); + + IssmDouble rho_ice = femmodel->parameters->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = femmodel->parameters->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble g = femmodel->parameters->FindParam(ConstantsGEnum); + + /*Constrain all nodes that are grounded and unconstrain the ones that float*/ + for(Object* & object : femmodel->elements->objects){ + + /*Get current element and return if not on base*/ + Element *element = xDynamicCast(object); + if(!element->IsOnBase()) continue; + + int numnodes = element->GetNumberOfNodes(); + IssmDouble *mask = xNew(numnodes); + IssmDouble *bed = xNew(numnodes); + IssmDouble *thickness = xNew(numnodes); + IssmDouble *ls_active = xNew(numnodes); + + element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum); + element->GetInputListOnNodes(&bed[0],BaseEnum); + element->GetInputListOnNodes(&thickness[0],ThicknessEnum); + element->GetInputListOnNodes(&ls_active[0],HydrologyMaskNodeActivationEnum); + + //for(int in=0;inGetNode(in); + if(mask[in]>0. && ls_active[in]==1.){ + node->Activate(); //Not sure if we need this! + } + else{ + IssmDouble phi = rho_ice*g*thickness[in] + rho_water*g*bed[in]; //FIXME this is correct! + node->Deactivate();// Not sure if we need this + node->ApplyConstraint(0,phi); + } + } + xDelete(mask); + xDelete(bed); + xDelete(thickness); + xDelete(ls_active); + } + return; }/*}}}*/ @@ -475,8 +537,9 @@ /*Skip if water or ice shelf element*/ - if(element->IsAllFloating()) return; + if(element->IsAllFloating() || !element->IsIceInElement()) return; + if(!element->IsOnBase()) return; /*Intermediaries */ - IssmDouble newgap = 0.; + IssmDouble newgap = 0.; IssmDouble Jdet,meltrate,G,dh[2],B,A,n,dt; IssmDouble gap,bed,thickness,head; @@ -484,7 +547,7 @@ IssmDouble alpha2,frictionheat; IssmDouble* xyz_list = NULL; - IssmDouble dpressure_water[2],dbed[2],PMPheat,dissipation; + IssmDouble dpressure_water[2],dbed[2],PMPheat,dissipation; IssmDouble q = 0.; - IssmDouble channelization = 0.; + IssmDouble channelization = 0.; /*Retrieve all inputs and parameters*/ @@ -551,6 +614,6 @@ if(pressure_water>pressure_ice) pressure_water = pressure_ice; - /* Compute change in sensible heat due to changes in pressure melting point*/ - dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); + /* Compute change in sensible heat due to changes in pressure melting point*/ + dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]); dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]); dissipation=rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]); @@ -558,15 +621,14 @@ meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])); - element->AddInput(DummyEnum,&meltrate,P0Enum); - element->AddInput(EsaEmotionEnum,&frictionheat,P0Enum); - element->AddInput(EsaNmotionEnum,&dissipation,P0Enum); - element->AddInput(EsaUmotionEnum,&PMPheat,P0Enum); - + element->AddBasalInput(DummyEnum,&meltrate,P0Enum); + element->AddBasalInput(EsaEmotionEnum,&frictionheat,P0Enum); + element->AddBasalInput(EsaNmotionEnum,&dissipation,P0Enum); + element->AddBasalInput(EsaUmotionEnum,&PMPheat,P0Enum); newgap += gauss->weight*Jdet*(gap+dt*( - meltrate/rho_ice - -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap - +beta*sqrt(vx*vx+vy*vy) - )); + meltrate/rho_ice + -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap + +beta*sqrt(vx*vx+vy*vy) + )); @@ -574,5 +636,5 @@ /* Compute basal water flux */ - q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])); + q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])); /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */ @@ -590,13 +652,13 @@ /*Add new gap as an input*/ - element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum); + element->AddBasalInput(HydrologyGapHeightEnum,&newgap,P0Enum); /*Divide by connectivity, add basal flux as an input*/ q = q/totalweights; - element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum); + element->AddBasalInput(HydrologyBasalFluxEnum,&q,P0Enum); /* Divide by connectivity, add degree of channelization as an input */ channelization = channelization/totalweights; - element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum); + element->AddBasalInput(DegreeOfChannelizationEnum,&channelization,P0Enum); /*Clean up and return*/ @@ -616,5 +678,6 @@ /*Skip if water or ice shelf element*/ - if(element->IsAllFloating()) return; + if(element->IsAllFloating() || !element->IsIceInElement()) return; + if(!element->IsOnBase()) return; /*Intermediaries*/ @@ -633,5 +696,4 @@ Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); - Gauss* gauss=element->NewGauss(); for (int i=0;iAddInput(EffectivePressureEnum,N,element->GetElementType()); + element->AddBasalInput(EffectivePressureEnum,N,element->GetElementType()); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 27648) +++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 27649) @@ -352,5 +352,11 @@ /*Initialize Load Vector and return if necessary*/ - Tria* tria=(Tria*)element; + Tria* tria=NULL; + if(element->ObjectEnum()==TriaEnum){ + tria = (Tria*)this->element; + } + else if(element->ObjectEnum()==PentaEnum){ + tria = (Tria*)this->element->SpawnBasalElement(); + } _assert_(tria->FiniteElement()==P1Enum); if(!tria->IsIceInElement() || tria->IsAllFloating()) return NULL; @@ -380,4 +386,5 @@ /*Clean up and return*/ delete gauss; + if(tria->IsSpawnedElement()){tria->DeleteMaterials(); delete tria;}; return pe; } Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 27648) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 27649) @@ -147,4 +147,6 @@ analysis_enum==HydrologyDCInefficientAnalysisEnum || analysis_enum==HydrologyDCEfficientAnalysisEnum || + analysis_enum==HydrologyShaktiAnalysisEnum || + analysis_enum==HydrologyGlaDSAnalysisEnum || analysis_enum==GLheightadvectionAnalysisEnum || analysis_enum==LevelsetAnalysisEnum Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27648) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27649) @@ -65,7 +65,5 @@ /*solid earth considerations:*/ SolidEarthWaterUpdates(femmodel); - delete ug; - } Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 27648) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 27649) @@ -117,5 +117,15 @@ if(ishydrology){ /*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/ - femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum); + int hydrology_model; + femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum); + if(hydrology_model==HydrologyshaktiEnum){ + femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum); + } + else if(hydrology_model==HydrologyGlaDSEnum){ + femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum); + } + else{ + _error_("hydrology model not supported yet"); + } }else if(isdebris){ femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);