 |
Ice Sheet System Model
4.18
Code documentation
|
#include <HydrologyDCEfficientAnalysis.h>
|
int | DofsPerNode (int **doflist, int domaintype, int approximation) |
|
void | UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum) |
|
void | UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type) |
|
void | CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false) |
|
void | CreateConstraints (Constraints *constraints, IoModel *iomodel) |
|
void | CreateLoads (Loads *loads, IoModel *iomodel) |
|
void | InitZigZagCounter (FemModel *femmodel) |
|
void | ResetCounter (FemModel *femmodel) |
|
void | Core (FemModel *femmodel) |
|
ElementVector * | CreateDVector (Element *element) |
|
ElementMatrix * | CreateJacobianMatrix (Element *element) |
|
ElementMatrix * | CreateKMatrix (Element *element) |
|
ElementVector * | CreatePVector (Element *element) |
|
void | GetB (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss) |
|
void | GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element) |
|
void | GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index) |
|
void | InputUpdateFromSolution (IssmDouble *solution, Element *element) |
|
void | UpdateConstraints (FemModel *femmodel) |
|
IssmDouble | EplStoring (Element *element, Gauss *gauss, Input2 *epl_thick_input, Input2 *epl_head_input, Input2 *base_input) |
|
IssmDouble | EplTransmitivity (Element *element, Gauss *gauss, Input2 *epl_thick_input, Input2 *epl_head_input, Input2 *base_input) |
|
void | GetHydrologyDCInefficientHmax (IssmDouble *ph_max, Element *element, Node *innode) |
|
IssmDouble | GetHydrologyKMatrixTransfer (Element *element) |
|
IssmDouble | GetHydrologyPVectorTransfer (Element *element, Gauss *gauss, Input2 *sed_head_input) |
|
void | ComputeEPLThickness (FemModel *femmodel) |
|
void | HydrologyEPLGetMask (Vector< IssmDouble > *vec_mask, Vector< IssmDouble > *recurence, Element *element) |
|
void | HydrologyEPLGetActive (Vector< IssmDouble > *active_vec, Element *element) |
|
virtual | ~Analysis () |
|
◆ DofsPerNode()
int HydrologyDCEfficientAnalysis::DofsPerNode |
( |
int ** |
doflist, |
|
|
int |
domaintype, |
|
|
int |
approximation |
|
) |
| |
|
virtual |
◆ UpdateParameters()
void HydrologyDCEfficientAnalysis::UpdateParameters |
( |
Parameters * |
parameters, |
|
|
IoModel * |
iomodel, |
|
|
int |
solution_enum, |
|
|
int |
analysis_enum |
|
) |
| |
|
virtual |
◆ UpdateElements()
void HydrologyDCEfficientAnalysis::UpdateElements |
( |
Elements * |
elements, |
|
|
Inputs2 * |
inputs2, |
|
|
IoModel * |
iomodel, |
|
|
int |
analysis_counter, |
|
|
int |
analysis_type |
|
) |
| |
|
virtual |
◆ CreateNodes()
void HydrologyDCEfficientAnalysis::CreateNodes |
( |
Nodes * |
nodes, |
|
|
IoModel * |
iomodel, |
|
|
bool |
isamr = false |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 70 of file HydrologyDCEfficientAnalysis.cpp.
74 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
78 bool isefficientlayer;
79 iomodel->
FindConstant(&isefficientlayer,
"md.hydrology.isefficientlayer");
80 if(!isefficientlayer)
return;
83 iomodel->
FetchData(2,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface");
86 iomodel->
DeleteData(2,
"md.mesh.vertexonbase",
"md.mesh.vertexonsurface");
◆ CreateConstraints()
void HydrologyDCEfficientAnalysis::CreateConstraints |
( |
Constraints * |
constraints, |
|
|
IoModel * |
iomodel |
|
) |
| |
|
virtual |
◆ CreateLoads()
void HydrologyDCEfficientAnalysis::CreateLoads |
( |
Loads * |
loads, |
|
|
IoModel * |
iomodel |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 102 of file HydrologyDCEfficientAnalysis.cpp.
106 iomodel->
FindConstant(&hydrology_model,
"md.hydrology.model");
110 bool isefficientlayer;
111 iomodel->
FindConstant(&isefficientlayer,
"md.hydrology.isefficientlayer");
112 if(!isefficientlayer)
return;
116 iomodel->
FetchData(1,
"md.mesh.vertexonbase");
128 else if(reCast<int>(iomodel->
Data(
"md.mesh.vertexonbase")[i])){
134 iomodel->
DeleteData(1,
"md.mesh.vertexonbase");
◆ InitZigZagCounter()
void HydrologyDCEfficientAnalysis::InitZigZagCounter |
( |
FemModel * |
femmodel | ) |
|
◆ ResetCounter()
void HydrologyDCEfficientAnalysis::ResetCounter |
( |
FemModel * |
femmodel | ) |
|
◆ Core()
void HydrologyDCEfficientAnalysis::Core |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
◆ CreateDVector()
◆ CreateJacobianMatrix()
◆ CreateKMatrix()
Implements Analysis.
Definition at line 165 of file HydrologyDCEfficientAnalysis.cpp.
176 basalelement = element;
179 if(!element->
IsOnBase())
return NULL;
188 if(!active_element) {
209 IssmDouble* basis = xNew<IssmDouble>(numnodes);
221 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
225 epl_transmitivity =
EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
226 epl_storing =
EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
229 D_scalar=epl_transmitivity*gauss->
weight*Jdet;
231 if(dt!=0.) D_scalar=D_scalar*dt;
234 GetB(B,basalelement,xyz_list,gauss);
243 D_scalar=epl_storing*gauss->
weight*Jdet;
252 D_scalar=dt*transfer*gauss->
weight*Jdet;
262 xDelete<IssmDouble>(xyz_list);
263 xDelete<IssmDouble>(basis);
264 xDelete<IssmDouble>(B);
◆ CreatePVector()
Implements Analysis.
Definition at line 269 of file HydrologyDCEfficientAnalysis.cpp.
280 basalelement = element;
283 if(!element->
IsOnBase())
return NULL;
292 if(!active_element) {
303 int hydrologysubstepping;
311 Input2 *old_wh_input = NULL;
312 Input2 *dummy_input = NULL;
313 Input2 *surface_runoff_input = NULL;
321 IssmDouble* basis = xNew<IssmDouble>(numnodes);
346 if(smbsubstepping==1){
350 else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
358 surface_runoff_input=xDynamicCast<Input2*>(dummy_input);
_assert_(surface_runoff_input);
363 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
367 epl_storing =
EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
368 epl_transmitivity =
EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
372 if(surface_runoff_input) surface_runoff_input->
GetInputValue(&runoff_value,gauss);
373 else runoff_value = 0.;
374 scalar = Jdet*gauss->
weight*(water_load+runoff_value);
376 if(dt!=0.) scalar = scalar*dt;
377 for(
int i=0;i<numnodes;i++)pe->
values[i]+=scalar*basis[i];
384 scalar = Jdet*gauss->
weight*((water_head*epl_storing)+(dt*transfer));
386 for(
int i=0;i<numnodes;i++)pe->
values[i]+=scalar*basis[i];
393 for(
int iv=0;iv<numvertices;iv++){
395 epl_transmitivity =
EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
398 pe->
values[iv]+=residual/connectivity;
402 xDelete<IssmDouble>(xyz_list);
403 xDelete<IssmDouble>(basis);
◆ GetB()
Definition at line 408 of file HydrologyDCEfficientAnalysis.cpp.
423 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
427 for(
int i=0;i<numnodes;i++){
428 B[numnodes*0+i] = dbasis[0*numnodes+i];
429 B[numnodes*1+i] = dbasis[1*numnodes+i];
433 xDelete<IssmDouble>(dbasis);
◆ GetSolutionFromInputs()
◆ GradientJ()
void HydrologyDCEfficientAnalysis::GradientJ |
( |
Vector< IssmDouble > * |
gradient, |
|
|
Element * |
element, |
|
|
int |
control_type, |
|
|
int |
control_index |
|
) |
| |
|
virtual |
◆ InputUpdateFromSolution()
void HydrologyDCEfficientAnalysis::InputUpdateFromSolution |
( |
IssmDouble * |
solution, |
|
|
Element * |
element |
|
) |
| |
|
virtual |
Implements Analysis.
Definition at line 441 of file HydrologyDCEfficientAnalysis.cpp.
450 basalelement = element;
466 IssmDouble* eplHeads = xNew<IssmDouble>(numnodes);
470 for(
int i=0;i<numnodes;i++){
471 eplHeads[i]=solution[doflist[i]];
472 if(xIsNan<IssmDouble>(eplHeads[i]))
_error_(
"NaN found in solution vector");
473 if(xIsInf<IssmDouble>(eplHeads[i]))
_error_(
"Inf found in solution vector");
479 xDelete<IssmDouble>(eplHeads);
480 xDelete<int>(doflist);
◆ UpdateConstraints()
void HydrologyDCEfficientAnalysis::UpdateConstraints |
( |
FemModel * |
femmodel | ) |
|
|
virtual |
◆ EplStoring()
Definition at line 489 of file HydrologyDCEfficientAnalysis.cpp.
492 IssmDouble epl_thickness,prestep_head,base_elev;
502 water_sheet=
max(0.0,(prestep_head-base_elev));
503 storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
◆ EplTransmitivity()
Definition at line 520 of file HydrologyDCEfficientAnalysis.cpp.
523 IssmDouble epl_thickness,base_elev,prestep_head;
529 water_sheet=
max(0.0,(prestep_head-base_elev));
530 epl_transmitivity=epl_conductivity*epl_thickness;
532 return epl_transmitivity;
◆ GetHydrologyDCInefficientHmax()
void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax |
( |
IssmDouble * |
ph_max, |
|
|
Element * |
element, |
|
|
Node * |
innode |
|
) |
| |
◆ GetHydrologyKMatrixTransfer()
IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer |
( |
Element * |
element | ) |
|
◆ GetHydrologyPVectorTransfer()
◆ ComputeEPLThickness()
void HydrologyDCEfficientAnalysis::ComputeEPLThickness |
( |
FemModel * |
femmodel | ) |
|
Definition at line 615 of file HydrologyDCEfficientAnalysis.cpp.
627 if(iseplthickcomp==0)
return;
637 IssmDouble* thickness = xNew<IssmDouble>(numnodes);
640 IssmDouble* eplhead = xNew<IssmDouble>(numnodes);
641 IssmDouble* epl_slopeX = xNew<IssmDouble>(numnodes);
642 IssmDouble* epl_slopeY = xNew<IssmDouble>(numnodes);
643 IssmDouble* old_thickness = xNew<IssmDouble>(numnodes);
644 IssmDouble* ice_thickness = xNew<IssmDouble>(numnodes);
662 default:
_error_(
"not Implemented Yet");
675 for(
int i=0;i<numnodes;i++){
676 thickness[i]=init_thick;
680 for(
int i=0;i<numnodes;i++){
683 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
684 if(EPL_N<0.0)EPL_N=0.0;
686 EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]);
688 opening=(rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat);
689 closing=(2.0*A*dt*pow(EPL_N,n[i]))/(pow(n[i],n[i]));
691 thickness[i] = old_thickness[i]/(1.0-opening+closing);
697 if(thickness[i]>max_thick){
698 thickness[i] = max_thick;
703 xDelete<IssmDouble>(thickness);
704 xDelete<IssmDouble>(eplhead);
705 xDelete<IssmDouble>(epl_slopeX);
706 xDelete<IssmDouble>(epl_slopeY);
707 xDelete<IssmDouble>(old_thickness);
708 xDelete<IssmDouble>(ice_thickness);
709 xDelete<IssmDouble>(bed);
710 xDelete<IssmDouble>(B);
711 xDelete<IssmDouble>(n);
◆ HydrologyEPLGetMask()
Definition at line 714 of file HydrologyDCEfficientAnalysis.cpp.
726 basalelement = element;
736 IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
737 IssmDouble* old_active =xNew<IssmDouble>(numnodes);
738 IssmDouble* sedhead =xNew<IssmDouble>(numnodes);
739 IssmDouble* eplhead =xNew<IssmDouble>(numnodes);
740 IssmDouble* residual =xNew<IssmDouble>(numnodes);
753 basalelement-> GetInputListOnVertices(&base[0],
BaseEnum);
756 sedheadmin=sedhead[0];
757 for(
int i=1;i<numnodes;i++)
if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
758 for(
int i=0;i<numnodes;i++){
761 if(old_active[i]>0.){
764 if(epl_thickness[i]<colapse_thick && residual[i]<=0.){
769 else if (old_active[i]==0.){
777 epl_thickness[i]=init_thick;
782 if(eplhead[i]>=h_max && active_element){
783 for(
int j=0;j<numnodes;j++){
785 if((sedhead[j] == sedheadmin) && (i!=j)){
787 if(old_active[j]==0.){
800 xDelete<IssmDouble>(epl_thickness);
801 xDelete<IssmDouble>(old_active);
802 xDelete<IssmDouble>(sedhead);
803 xDelete<IssmDouble>(eplhead);
804 xDelete<IssmDouble>(residual);
805 xDelete<IssmDouble>(base);
◆ HydrologyEPLGetActive()
Definition at line 808 of file HydrologyDCEfficientAnalysis.cpp.
817 basalelement = element;
828 IssmDouble* active = xNew<IssmDouble>(numnodes);
835 for(
int i=0;i<numnodes;i++) flag+=active[i];
839 for(
int i=0;i<numnodes;i++){
844 else if(active_element){
845 for(
int i=0;i<numnodes;i++){
853 xDelete<IssmDouble>(active);
The documentation for this class was generated from the following files:
virtual int GetElementType(void)=0
@ HydrologydcEplColapseThicknessEnum
virtual int GetNumberOfNodes(void)=0
void FindParam(bool *pvalue, int paramenum)
int AddObject(Object *object)
@ HydrologydcMaskEplactiveNodeEnum
@ HydrologydcEplInitialThicknessEnum
void GetHydrologyDCInefficientHmax(IssmDouble *ph_max, Element *element, Node *innode)
void GetB(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
@ MaterialsRhoFreshwaterEnum
@ TimesteppingTimeStepEnum
virtual Element * SpawnBasalElement(void)=0
@ HydrologydcEplThicknessSubstepEnum
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
@ MaterialsLatentheatEnum
IssmDouble GetHydrologyKMatrixTransfer(Element *element)
@ HydrologydcEplThicknessOldEnum
void AddObject(Param *newparam)
@ MaterialsRheologyBbarEnum
virtual Input2 * GetInput2(int inputenum)=0
void DeleteMaterials(void)
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
@ SedimentHeadResidualEnum
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
void SetParam(bool boolean, int enum_type)
void DeleteData(int num,...)
@ BasalforcingsGroundediceMeltingRateEnum
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
virtual Gauss * NewGauss(void)=0
const char * EnumToStringx(int enum_in)
@ SMBgradientscomponentsEnum
void FindConstant(bool *pvalue, const char *constant_name)
@ HydrologydcLeakageFactorEnum
void GetVerticesCoordinates(IssmDouble **xyz_list)
@ HydrologydcEplThickCompEnum
@ HydrologydcEplMaxThicknessEnum
@ HydrologydcEplConductivityEnum
@ HydrologydcSedimentlimitEnum
void FetchData(bool *pboolean, const char *data_name)
IssmDouble EplStoring(Element *element, Gauss *gauss, Input2 *epl_thick_input, Input2 *epl_head_input, Input2 *base_input)
@ HydrologyDCEfficientAnalysisEnum
void CreateSingleNodeToElementConnectivity(IoModel *iomodel)
@ HydrologydcWaterCompressibilityEnum
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
@ HydrologydcEplflipLockEnum
@ HydrologydcSedimentlimitFlagEnum
virtual void AddBasalInput2(int input_enum, IssmDouble *values, int interpolation_enum)
virtual void GaussVertex(int iv)=0
@ SedimentHeadSubstepEnum
IssmDouble * Data(const char *data_name)
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
IssmDouble GetHydrologyPVectorTransfer(Element *element, Gauss *gauss, Input2 *sed_head_input)
#define _error_(StreamArgs)
virtual int begin(void)=0
Object * GetObjectByOffset(int offset)
@ MeshVertexonsurfaceEnum
IssmDouble EplTransmitivity(Element *element, Gauss *gauss, Input2 *epl_thick_input, Input2 *epl_head_input, Input2 *base_input)
@ HydrologydcBasalMoulinInputEnum
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
virtual void GaussPoint(int ig)=0
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
void GetInputValue(bool *pvalue, int enum_type)
@ HydrologyStepsPerStepEnum
void FindParam(bool *pinteger, int enum_type)
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
@ HydrologydcTransferFlagEnum
virtual int VertexConnectivity(int vertexindex)=0
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
@ HydrologydcEplCompressibilityEnum
IssmDouble max(IssmDouble a, IssmDouble b)
virtual int GetNumberOfVertices(void)=0
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
@ HydrologydcEplPorosityEnum
@ HydrologydcMaskEplactiveEltEnum
void SetValue(int dof, doubletype value, InsMode mode)
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)