Changeset 14757
- Timestamp:
- 04/25/13 13:07:40 (12 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 1 added
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h
r14744 r14757 90 90 HydrologydcEnum, 91 91 SedimentHeadEnum, 92 EplHeadEnum, 92 93 HydrologydcSpcsedimentHeadEnum, 93 94 HydrologydcSedimentCompressibilityEnum, -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
r14746 r14757 5799 5799 ElementMatrix* Tria::CreateKMatrixHydrology(void){ 5800 5800 5801 int hydrology_model; 5801 int hydrology_model; 5802 bool isefficientlayer; 5802 5803 parameters->FindParam(&hydrology_model,HydrologyEnum); 5804 5803 5805 5804 5806 switch(hydrology_model){ … … 5806 5808 return CreateKMatrixHydrologyShreve(); 5807 5809 case HydrologydcEnum: 5808 return CreateKMatrixHydrologyDC(); 5810 parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 5811 if(isefficientlayer) 5812 return CreateKMatrixHydrologyDCsediment(); 5813 else 5814 return CreateKMatrixHydrologyDCepl(); 5809 5815 default: 5810 5816 _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); … … 5917 5923 } 5918 5924 /*}}}*/ 5919 /*FUNCTION Tria::CreateKMatrixHydrologyDC {{{*/5920 ElementMatrix* Tria::CreateKMatrixHydrologyDC (void){5925 /*FUNCTION Tria::CreateKMatrixHydrologyDCsediment{{{*/ 5926 ElementMatrix* Tria::CreateKMatrixHydrologyDCsediment(void){ 5921 5927 5922 5928 /*constants: */ … … 5978 5984 } 5979 5985 /*}}}*/ 5986 /*FUNCTION Tria::CreateKMatrixHydrologyDCepl{{{*/ 5987 ElementMatrix* Tria::CreateKMatrixHydrologyDCepl(void){ 5988 5989 /*constants: */ 5990 const int numdof=NDOF1*NUMVERTICES; 5991 5992 /* Intermediaries */ 5993 IssmDouble D_scalar,Jdet; 5994 IssmDouble epl_transmitivity,dt; 5995 IssmDouble epl_storing; 5996 IssmDouble xyz_list[NUMVERTICES][3]; 5997 IssmDouble B[2][numdof]; 5998 IssmDouble L[NUMVERTICES]; 5999 IssmDouble D[2][2]; 6000 GaussTria *gauss = NULL; 6001 6002 /*Initialize Element matrix*/ 6003 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); 6004 6005 /*Retrieve all inputs and parameters*/ 6006 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); 6007 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); 6008 epl_storing = matpar->GetEplStoring(); 6009 epl_transmitivity = matpar->GetEplTransmitivity(); 6010 6011 /* Start looping on the number of gaussian points: */ 6012 gauss=new GaussTria(2); 6013 for(int ig=gauss->begin();ig<gauss->end();ig++){ 6014 6015 gauss->GaussPoint(ig); 6016 6017 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 6018 6019 /*Diffusivity*/ 6020 D_scalar=epl_transmitivity*gauss->weight*Jdet; 6021 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt; 6022 D[0][0]=D_scalar; D[0][1]=0.; 6023 D[1][0]=0.; D[1][1]=D_scalar; 6024 GetBHydro(&B[0][0],&xyz_list[0][0],gauss); 6025 TripleMultiply(&B[0][0],2,numdof,1, 6026 &D[0][0],2,2,0, 6027 &B[0][0],2,numdof,0, 6028 &Ke->values[0],1); 6029 6030 /*Transient*/ 6031 if(reCast<bool,IssmDouble>(dt)){ 6032 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); 6033 D_scalar=epl_storing*gauss->weight*Jdet; 6034 6035 TripleMultiply(&L[0],numdof,1,0, 6036 &D_scalar,1,1,0, 6037 &L[0],1,numdof,0, 6038 &Ke->values[0],1); 6039 } 6040 } 6041 6042 /*Clean up and return*/ 6043 delete gauss; 6044 return Ke; 6045 } 6046 /*}}}*/ 5980 6047 /*FUNCTION Tria::CreatePVectorHydrology{{{*/ 5981 6048 ElementVector* Tria::CreatePVectorHydrology(void){ 5982 6049 5983 int hydrology_model; 6050 int hydrology_model; 6051 bool isefficientlayer; 5984 6052 parameters->FindParam(&hydrology_model,HydrologyEnum); 5985 6053 … … 5988 6056 return CreatePVectorHydrologyShreve(); 5989 6057 case HydrologydcEnum: 5990 return CreatePVectorHydrologyDC(); 6058 parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 6059 if(isefficientlayer) 6060 return CreatePVectorHydrologyDCsediment(); 6061 else 6062 return CreatePVectorHydrologyDCepl(); 5991 6063 default: 5992 6064 _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet"); … … 6044 6116 } 6045 6117 /*}}}*/ 6046 /*FUNCTION Tria::CreatePVectorHydrologyDC {{{*/6047 ElementVector* Tria::CreatePVectorHydrologyDC (void){6118 /*FUNCTION Tria::CreatePVectorHydrologyDCsediment {{{*/ 6119 ElementVector* Tria::CreatePVectorHydrologyDCsediment(void){ 6048 6120 6049 6121 /*Constants*/ … … 6092 6164 old_wh_input->GetInputValue(&water_head,gauss); 6093 6165 scalar = Jdet*gauss->weight*water_head*sediment_storing; 6166 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i]; 6167 } 6168 } 6169 6170 /*Clean up and return*/ 6171 delete gauss; 6172 return pe; 6173 } 6174 /*}}}*/ 6175 /*FUNCTION Tria::CreatePVectorHydrologyDCepl {{{*/ 6176 ElementVector* Tria::CreatePVectorHydrologyDCepl(void){ 6177 6178 /*Constants*/ 6179 const int numdof=NDOF1*NUMVERTICES; 6180 6181 /*Intermediaries */ 6182 IssmDouble Jdet; 6183 IssmDouble xyz_list[NUMVERTICES][3]; 6184 IssmDouble dt,scalar,water_head; 6185 IssmDouble water_load; 6186 IssmDouble epl_storing; 6187 IssmDouble basis[numdof]; 6188 GaussTria* gauss=NULL; 6189 6190 /*Initialize Element vector*/ 6191 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); 6192 6193 /*Retrieve all inputs and parameters*/ 6194 epl_storing = matpar->GetEplStoring(); 6195 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 6196 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); 6197 Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input); 6198 Input* old_wh_input=NULL; 6199 6200 if(reCast<bool,IssmDouble>(dt)){ 6201 old_wh_input=inputs->GetInput(EplHeadEnum); _assert_(old_wh_input); 6202 } 6203 6204 /* Start looping on the number of gaussian points: */ 6205 gauss=new GaussTria(2); 6206 for(int ig=gauss->begin();ig<gauss->end();ig++){ 6207 6208 gauss->GaussPoint(ig); 6209 6210 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 6211 GetNodalFunctions(basis, gauss); 6212 6213 /*Loading term*/ 6214 water_input->GetInputValue(&water_load,gauss); 6215 scalar = Jdet*gauss->weight*water_load; 6216 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt; 6217 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i]; 6218 6219 /*Transient term*/ 6220 if(reCast<bool,IssmDouble>(dt)){ 6221 old_wh_input->GetInputValue(&water_head,gauss); 6222 scalar = Jdet*gauss->weight*water_head*epl_storing; 6094 6223 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i]; 6095 6224 } -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h
r14695 r14757 240 240 ElementMatrix* CreateKMatrixHydrology(void); 241 241 ElementMatrix* CreateKMatrixHydrologyShreve(void); 242 ElementMatrix* CreateKMatrixHydrologyDC(void); 242 ElementMatrix* CreateKMatrixHydrologyDCsediment(void); 243 ElementMatrix* CreateKMatrixHydrologyDCepl(void); 243 244 ElementVector* CreatePVectorHydrology(void); 244 245 ElementVector* CreatePVectorHydrologyShreve(void); 245 ElementVector* CreatePVectorHydrologyDC(void); 246 ElementVector* CreatePVectorHydrologyDCsediment(void); 247 ElementVector* CreatePVectorHydrologyDCepl(void); 246 248 void CreateHydrologyWaterVelocityInput(void); 247 249 void GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution); -
issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp
r14743 r14757 63 63 iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum); 64 64 iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum); 65 66 iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum); 67 iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum); 68 iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum); 69 iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum); 65 70 } 66 71 else{ … … 355 360 } 356 361 /*}}}*/ 362 /*FUNCTION Matpar::GetEplStoring {{{*/ 363 IssmDouble Matpar::GetEplStoring(){ 364 return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness* 365 ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity)); 366 } 367 /*}}}*/ 357 368 /*FUNCTION Matpar::GetSedimentTransitivity {{{*/ 358 369 IssmDouble Matpar::GetSedimentTransmitivity(){ 359 370 return sediment_transmitivity; 371 } 372 /*}}}*/ 373 /*FUNCTION Matpar::GetEplTransitivity {{{*/ 374 IssmDouble Matpar::GetEplTransmitivity(){ 375 return epl_transmitivity; 360 376 } 361 377 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h
r14735 r14757 45 45 IssmDouble sediment_transmitivity; 46 46 IssmDouble water_compressibility; 47 48 IssmDouble epl_compressibility; 49 IssmDouble epl_porosity; 50 IssmDouble epl_thickness; 51 IssmDouble epl_transmitivity; 47 52 48 53 /*gia: */ … … 118 123 IssmDouble GetHydrologyN(); 119 124 IssmDouble GetSedimentStoring(); 125 IssmDouble GetEplStoring(); 120 126 IssmDouble GetSedimentTransmitivity(); 127 IssmDouble GetEplTransmitivity(); 121 128 IssmDouble TMeltingPoint(IssmDouble pressure); 122 129 IssmDouble PureIceEnthalpy(IssmDouble pressure); -
issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp
r14744 r14757 95 95 case HydrologydcEnum : return "Hydrologydc"; 96 96 case SedimentHeadEnum : return "SedimentHead"; 97 case EplHeadEnum : return "EplHead"; 97 98 case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead"; 98 99 case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility"; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp
r14562 r14757 16 16 Parameters *parameters = NULL; 17 17 int hydrology_model; 18 bool isefficientlayer; 18 19 19 20 /*Get parameters: */ … … 28 29 } 29 30 else if(hydrology_model==HydrologydcEnum){ 30 //No parameters for now 31 iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum); 32 parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); 31 33 } 32 34 else{ -
issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp
r14744 r14757 96 96 else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum; 97 97 else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; 98 else if (strcmp(name,"EplHead")==0) return EplHeadEnum; 98 99 else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum; 99 100 else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum; … … 137 138 else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum; 138 139 else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; 139 else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;140 140 else stage=2; 141 141 } 142 142 if(stage==2){ 143 if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; 143 if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum; 144 else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; 144 145 else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; 145 146 else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum; … … 260 261 else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; 261 262 else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; 262 else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;263 263 else stage=3; 264 264 } 265 265 if(stage==3){ 266 if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; 266 if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum; 267 else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; 267 268 else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum; 268 269 else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; … … 383 384 else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; 384 385 else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; 385 else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;386 386 else stage=4; 387 387 } 388 388 if(stage==4){ 389 if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; 389 if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; 390 else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; 390 391 else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; 391 392 else if (strcmp(name,"Boundary")==0) return BoundaryEnum; … … 506 507 else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; 507 508 else if (strcmp(name,"Incremental")==0) return IncrementalEnum; 508 else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;509 509 else stage=5; 510 510 } 511 511 if(stage==5){ 512 if (strcmp(name,"None")==0) return NoneEnum; 512 if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum; 513 else if (strcmp(name,"None")==0) return NoneEnum; 513 514 else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; 514 515 else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; -
issm/trunk-jpl/src/c/solutions/hydrology_core.cpp
r14591 r14757 78 78 } 79 79 } 80 80 81 else if (hydrology_model==HydrologydcEnum){ 81 82 if(VerboseSolution()) _pprintLine_(" computing water head"); -
issm/trunk-jpl/src/m/classes/initialization.m
r14643 r14757 6 6 classdef initialization 7 7 properties (SetAccess=public) 8 vx = NaN; 9 vy = NaN; 10 vz = NaN; 11 vel = NaN; 12 pressure = NaN; 13 temperature = NaN; 14 waterfraction = NaN; 15 sediment_head = NaN; 16 watercolumn = NaN; 8 vx = NaN; 9 vy = NaN; 10 vz = NaN; 11 vel = NaN; 12 pressure = NaN; 13 temperature = NaN; 14 waterfraction = NaN; 15 sediment_head = NaN; 16 epl_head = NaN; 17 watercolumn = NaN; 17 18 end 18 19 methods … … 59 60 if isa(md.hydrology,'hydrologydc'), 60 61 md = checkfield(md,'initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]); 62 md = checkfield(md,'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]); 61 63 else 62 64 md = checkfield(md,'initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices 1]); … … 75 77 fielddisplay(obj,'waterfraction','fraction of water in the ice'); 76 78 fielddisplay(obj,'sediment_head','sediment water head of subglacial system [m]'); 79 fielddisplay(obj,'epl_head','epl water head of subglacial system [m]'); 77 80 fielddisplay(obj,'watercolumn','thickness of subglacial water [m]'); 78 81 … … 86 89 WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum); 87 90 WriteData(fid,'data',obj.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum); 91 WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum); 88 92 WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum); 89 93 end % }}} -
issm/trunk-jpl/src/m/enum/EnumDefinitions.py
r14744 r14757 789 789 return StringToEnum('SedimentHead')[0] 790 790 791 def EplHeadEnum(): 792 """ 793 EPLHEADENUM - Enum of EplHead 794 795 Usage: 796 macro=EplHeadEnum() 797 """ 798 799 return StringToEnum('EplHead')[0] 800 791 801 def HydrologydcSpcsedimentHeadEnum(): 792 802 """ … … 5337 5347 """ 5338 5348 5339 return 53 25340 5349 return 533 5350 -
issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m
r14744 r14757 9 9 % macro=MaximumNumberOfEnums() 10 10 11 macro=53 2;11 macro=533;
Note:
See TracChangeset
for help on using the changeset viewer.