![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21759 for issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
- Timestamp:
- 06/07/17 10:50:54 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
r21233 r21759 9 9 return 1; 10 10 }/*}}}*/ 11 11 12 void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ 12 13 … … 33 34 iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp"); 34 35 parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp)); 35 36 37 }/*}}}*/ 36 }/*}}}*/ 37 38 38 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ 39 39 … … 83 83 if(!isefficientlayer) return; 84 84 85 if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); 85 if(iomodel->domaintype!=Domain2DhorizontalEnum){ 86 iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); 87 } 86 88 ::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum); 87 89 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); … … 104 106 105 107 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ 106 /*Nothing for now*/ 107 108 /*Fetch parameters: */ 108 109 /*Do we really want DC?*/ 109 110 int hydrology_model; 110 111 iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); 111 112 if(hydrology_model!=HydrologydcEnum) return; 112 113 113 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase"); 114 115 //create penalties for nodes: no node can have water above the max 114 /*Do we want an efficient layer*/ 115 bool isefficientlayer; 116 iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); 117 if(!isefficientlayer) return; 118 119 /*Fetch parameters: */ 120 if(iomodel->domaintype==Domain3DEnum){ 121 iomodel->FetchData(1,"md.mesh.vertexonbase"); 122 } 123 124 //Add moulin inputs as loads 116 125 CreateSingleNodeToElementConnectivity(iomodel); 117 126 for(int i=0;i<iomodel->numberofvertices;i++){ … … 132 141 133 142 void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/ 143 134 144 int* eplzigzag_counter =NULL; 135 136 145 eplzigzag_counter=xNewZeroInit<int>(femmodel->nodes->Size()); 137 146 femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size())); … … 142 151 143 152 int* eplzigzag_counter=NULL; 144 Element* element=NULL;145 146 153 femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 147 154 for(int i=0;i<femmodel->nodes->Size();i++){ 148 149 155 eplzigzag_counter[i]=0; 150 156 } … … 201 207 /* Intermediaries */ 202 208 IssmDouble D_scalar,Jdet,dt; 203 IssmDouble epl_thickness;204 209 IssmDouble transfer; 210 IssmDouble epl_transmitivity; 211 IssmDouble epl_storing; 205 212 IssmDouble *xyz_list = NULL; 206 213 … … 219 226 220 227 Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input); 221 Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); 222 Input* epl_head_input = basalelement->GetInput(EplHeadEnum); 223 Input* thick_input = basalelement->GetInput(ThicknessEnum); 224 Input* base_input = basalelement->GetInput(BaseEnum); 225 226 IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); 227 IssmDouble epl_conductivity = basalelement->GetMaterialParameter(HydrologydcEplConductivityEnum); 228 Input* epl_head_input = basalelement->GetInput(EplHeadEnum); _assert_(epl_head_input); 229 Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); 228 230 229 231 /* Start looping on the number of gaussian points: */ 230 Gauss* gauss =basalelement->NewGauss(2);232 Gauss* gauss = basalelement->NewGauss(2); 231 233 for(int ig=gauss->begin();ig<gauss->end();ig++){ 232 234 gauss ->GaussPoint(ig); 233 235 basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); 234 epl_thick_input ->GetInputValue(&epl_thickness,gauss); 236 237 epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input); 238 epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); 235 239 236 240 /*Diffusivity*/ 237 D_scalar=epl_ conductivity*epl_thickness*gauss->weight*Jdet;241 D_scalar=epl_transmitivity*gauss->weight*Jdet; 238 242 if(dt!=0.) D_scalar=D_scalar*dt; 239 243 D[0][0]=D_scalar; … … 248 252 if(dt!=0.){ 249 253 basalelement->NodalFunctions(&basis[0],gauss); 250 D_scalar=epl_s pecificstoring*epl_thickness*gauss->weight*Jdet;254 D_scalar=epl_storing*gauss->weight*Jdet; 251 255 TripleMultiply(basis,numnodes,1,0, 252 256 &D_scalar,1,1,0, … … 256 260 257 261 /*Transfer EPL part*/ 258 transfer=GetHydrologyKMatrixTransfer(basalelement ,gauss,sed_head_input,epl_head_input,thick_input,base_input);262 transfer=GetHydrologyKMatrixTransfer(basalelement); 259 263 D_scalar=dt*transfer*gauss->weight*Jdet; 260 264 TripleMultiply(basis,numnodes,1,0, … … 299 303 /*Check that all nodes are active, else return empty matrix*/ 300 304 if(!active_element) { 301 if(domaintype!=Domain2DhorizontalEnum){305 if(domaintype!=Domain2DhorizontalEnum){ 302 306 basalelement->DeleteMaterials(); 303 307 delete basalelement; … … 308 312 IssmDouble dt,scalar,water_head; 309 313 IssmDouble water_load,transfer; 310 IssmDouble epl_ thickness;314 IssmDouble epl_storing; 311 315 IssmDouble Jdet; 312 316 IssmDouble residual,connectivity; … … 328 332 329 333 Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessEnum); _assert_(epl_thick_input); 330 Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); 331 Input* epl_head_input = basalelement->GetInput(EplHeadEnum); 332 Input* thick_input = basalelement->GetInput(ThicknessEnum); 333 Input* base_input = basalelement->GetInput(BaseEnum); 334 Input* sed_head_input = basalelement->GetInput(SedimentHeadEnum); _assert_(sed_head_input); 335 Input* epl_head_input = basalelement->GetInput(EplHeadEnum); _assert_(epl_head_input); 334 336 Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); 335 Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); 336 if(dt!= 0.){old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input);} 337 338 IssmDouble epl_specificstoring = EplSpecificStoring(basalelement); 339 337 Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); 338 Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); 339 340 if(dt!= 0.){ 341 old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input); 342 } 340 343 /* Start looping on the number of gaussian points: */ 341 Gauss* gauss =basalelement->NewGauss(2);344 Gauss* gauss = basalelement->NewGauss(2); 342 345 for(int ig=gauss->begin();ig<gauss->end();ig++){ 343 346 gauss->GaussPoint(ig); 344 345 347 basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); 346 348 basalelement ->NodalFunctions(basis,gauss); 347 349 epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); 348 350 /*Loading term*/ 349 351 water_input->GetInputValue(&water_load,gauss); 350 352 scalar = Jdet*gauss->weight*(water_load); 351 353 if(dt!=0.) scalar = scalar*dt; 352 for(int i=0;i<numnodes;i++){ 353 pe->values[i]+=scalar*basis[i]; 354 } 354 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i]; 355 355 356 356 /*Transient and transfer terms*/ 357 357 if(dt!=0.){ 358 old_wh_input ->GetInputValue(&water_head,gauss); 359 epl_thick_input ->GetInputValue(&epl_thickness,gauss); 360 358 old_wh_input->GetInputValue(&water_head,gauss); 361 359 /*Dealing with the epl part of the transfer term*/ 362 transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input ,epl_head_input,thick_input,base_input);363 scalar = Jdet*gauss->weight*((water_head*epl_s pecificstoring*epl_thickness)+(dt*transfer));360 transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input); 361 scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer)); 364 362 for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i]; 365 363 } … … 371 369 for(int iv=0;iv<numvertices;iv++){ 372 370 gauss->GaussVertex(iv); 373 374 371 connectivity = IssmDouble(basalelement->VertexConnectivity(iv)); 375 372 residual_input->GetInputValue(&residual,gauss); 376 373 pe->values[iv]+=residual/connectivity; 377 378 } 379 374 } 380 375 /*Clean up and return*/ 381 376 xDelete<IssmDouble>(xyz_list); … … 395 390 396 391 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ 397 398 bool active_element; 399 int domaintype,i; 400 Element* basalelement=NULL; 401 392 /*Intermediaries*/ 393 bool active_element; 394 int domaintype; 395 Element* basalelement=NULL; 396 397 /*Get basal element*/ 402 398 element->FindParam(&domaintype,DomainTypeEnum); 403 404 if(domaintype!=Domain2DhorizontalEnum){ 405 if(!element->IsOnBase()) return; 406 basalelement=element->SpawnBasalElement(); 407 } 408 else{ 409 basalelement = element; 410 } 411 399 switch(domaintype){ 400 case Domain2DhorizontalEnum: 401 basalelement = element; 402 break; 403 case Domain3DEnum: 404 if(!element->IsOnBase()) return; 405 basalelement = element->SpawnBasalElement(); 406 break; 407 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 408 } 412 409 /*Intermediary*/ 413 410 int* doflist = NULL; … … 417 414 418 415 /*Fetch dof list and allocate solution vector*/ 419 IssmDouble* sedhead = xNew<IssmDouble>(numnodes);420 416 IssmDouble* eplHeads = xNew<IssmDouble>(numnodes); 421 422 417 basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 423 basalelement->GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);424 418 425 419 Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); … … 427 421 428 422 /*Use the dof list to index into the solution vector: */ 429 for(i=0;i<numnodes;i++){ 430 eplHeads[i]=solution[doflist[i]]; 431 if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector"); 432 if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector"); 423 /*If the EPL is not active we revert to the bedrock elevation*/ 424 if(active_element){ 425 for(int i=0;i<numnodes;i++){ 426 eplHeads[i]=solution[doflist[i]]; 427 if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector"); 428 if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector"); 429 } 430 } 431 else{ 432 basalelement->GetInputListOnVertices(&eplHeads[0],BaseEnum); 433 for(int i=0;i<numnodes;i++){ 434 if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector"); 435 if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector"); 436 } 433 437 } 434 438 /*Add input to the element: */ 435 439 element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum); 436 437 440 /*Free ressources:*/ 438 441 xDelete<IssmDouble>(eplHeads); 439 xDelete<IssmDouble>(sedhead);440 442 xDelete<int>(doflist); 441 443 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; … … 448 450 449 451 /*Intermediaries*/ 450 IssmDouble HydrologyDCEfficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/ 451 IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); 452 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); 453 IssmDouble epl_porosity = element->GetMaterialParameter(HydrologydcEplPorosityEnum); 454 IssmDouble epl_compressibility = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum); 455 IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); 456 return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity)); 457 }/*}}}*/ 458 459 IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/ 460 IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); 461 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); 462 IssmDouble sediment_porosity = element->GetMaterialParameter(HydrologydcSedimentPorosityEnum); 463 IssmDouble sediment_thickness = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum); 464 IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum); 465 IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); 466 return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); 467 }/*}}}*/ 468 469 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ 452 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/ 470 453 471 454 int transfermethod; 472 IssmDouble hmax;473 IssmDouble epl_head,sediment_head;474 455 IssmDouble leakage,transfer; 475 IssmDouble continuum, factor;476 HydrologyDCInefficientAnalysis* inefanalysis = NULL;477 456 478 457 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); 479 458 /*Switch between the different transfer methods cases*/ 459 switch(transfermethod){ 460 case 0: 461 /*Just keepping the transfer to zero*/ 462 transfer=0.0; 463 break; 464 case 1: 465 element->FindParam(&leakage,HydrologydcLeakageFactorEnum); 466 transfer=leakage; 467 break; 468 default: 469 _error_("no case higher than 1 for the Transfer method"); 470 } 471 return transfer; 472 }/*}}}*/ 473 474 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input){/*{{{*/ 475 476 int transfermethod; 477 IssmDouble sediment_head; 478 IssmDouble leakage,transfer; 479 480 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); 480 481 /*Switch between the different transfer methods cases*/ 481 482 switch(transfermethod){ … … 486 487 case 1: 487 488 _assert_(sed_head_input); 488 _assert_(epl_head_input);489 490 inefanalysis = new HydrologyDCInefficientAnalysis();491 hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input);492 delete inefanalysis;493 494 489 sed_head_input->GetInputValue(&sediment_head,gauss); 495 epl_head_input->GetInputValue(&epl_head,gauss);496 490 element->FindParam(&leakage,HydrologydcLeakageFactorEnum); 497 498 //Computing continuum function to apply to transfer term, transfer is null only if 499 // epl_head>sediment_head AND sediment_head>h_max 500 continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); 501 factor=min(continuum,1.0); 502 transfer=leakage*factor; 503 break; 504 default: 505 _error_("no case higher than 1 for the Transfer method"); 506 } 507 return transfer; 508 }/*}}}*/ 509 510 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ 511 512 int transfermethod; 513 IssmDouble hmax; 514 IssmDouble epl_head,sediment_head; 515 IssmDouble leakage,transfer; 516 IssmDouble continuum, factor; 517 518 HydrologyDCInefficientAnalysis* inefanalysis = NULL; 519 520 element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); 521 522 /*Switch between the different transfer methods cases*/ 523 switch(transfermethod){ 524 case 0: 525 /*Just keepping the transfer to zero*/ 526 transfer=0.0; 527 break; 528 case 1: 529 _assert_(sed_head_input); 530 _assert_(epl_head_input); 531 532 inefanalysis = new HydrologyDCInefficientAnalysis(); 533 hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input); 534 delete inefanalysis; 535 536 sed_head_input->GetInputValue(&sediment_head,gauss); 537 epl_head_input->GetInputValue(&epl_head,gauss); 538 element->FindParam(&leakage,HydrologydcLeakageFactorEnum); 539 540 //Computing continuum function to apply to transfer term, transfer is null only if 541 // epl_head>sediment_head AND sediment_head>h_max 542 continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); 543 factor=min(continuum,1.0); 544 transfer=sediment_head*leakage*factor; 545 491 transfer=sediment_head*leakage; 546 492 break; 547 493 default: … … 560 506 IssmDouble EPLgrad2; 561 507 IssmDouble EPL_N; 562 563 508 564 509 femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); 510 femmodel->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum); 511 if(iseplthickcomp==0) return; 565 512 566 513 for(int j=0;j<femmodel->elements->Size();j++){ 567 514 568 515 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); 569 element->parameters->FindParam(&iseplthickcomp,HydrologydcEplThickCompEnum);570 if(iseplthickcomp==0) return;571 516 572 517 switch(domaintype){ … … 615 560 616 561 if(!active_element){ 617 618 562 /*Keeping thickness to initial value if EPL is not active*/ 619 563 for(int i=0;i<numnodes;i++){ … … 623 567 else{ 624 568 for(int i=0;i<numnodes;i++){ 625 626 569 /*Compute first the effective pressure in the EPL*/ 627 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water* max(0.0,(eplhead[i]-bed[i]))));570 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i]))); 628 571 if(EPL_N<0.0)EPL_N=0.0; 629 572 /*Get then the square of the gradient of EPL heads*/ 630 573 EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]); 631 632 574 /*And proceed to the real thing*/ 633 thickness[i] = old_thickness[i] *(1.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-634 (2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))));635 575 thickness[i] = old_thickness[i]/(1.0 576 -((rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat)) 577 +((2.0*A*dt*pow(EPL_N,n))/(pow(n,n)))); 636 578 /*Take care of otherthikening*/ 637 579 if(thickness[i]>max_thick){ … … 682 624 683 625 bool active_element; 684 int i,j;685 626 int domaintype; 686 627 IssmDouble h_max; … … 707 648 IssmDouble* eplhead =xNew<IssmDouble>(numnodes); 708 649 IssmDouble* residual =xNew<IssmDouble>(numnodes); 650 IssmDouble* base =xNew<IssmDouble>(numnodes); 709 651 710 652 IssmDouble init_thick =basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum); … … 719 661 basalelement-> GetInputListOnVertices(&eplhead[0],EplHeadEnum); 720 662 basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); 663 basalelement-> GetInputListOnVertices(&base[0],BaseEnum); 721 664 722 665 /*Get minimum sediment head of the element*/ 723 666 sedheadmin=sedhead[0]; 724 for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i]; 725 for(i=0;i<numnodes;i++){ 726 /*If node is now closed bring its thickness back to initial*/ 727 if (old_active[i]==0.){ 728 epl_thickness[i]=init_thick; 729 } 730 731 /*Now starting to look at the activations*/ 732 if(residual[i]>0.){ 667 for(int i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i]; 668 for(int i=0;i<numnodes;i++){ 669 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]); 670 /*If mask was already one, keep one or colapse*/ 671 if(old_active[i]>0.){ 733 672 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 734 if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 735 } 736 /*If mask was already one, keep one or colapse*/ 737 else if(old_active[i]>0.){ 738 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 739 /*If epl thickness gets under colapse thickness, close the layer*/ 673 /* If epl thickness gets under colapse thickness, close the layer */ 740 674 if(epl_thickness[i]<colapse_thick){ 741 675 vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); 742 676 recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 743 677 } 678 /* //If epl head gets under base elevation, close the layer */ 679 /* else if(eplhead[i]<(base[i]-1.0e-8)){ */ 680 /* vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); */ 681 /* recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); */ 682 /* } */ 683 } 684 /*If node is now closed bring its thickness back to initial*/ 685 if (old_active[i]==0.){ 686 epl_thickness[i]=init_thick; 687 /*Activate if we have a residual from sediment*/ 688 if(residual[i]>0.){ 689 vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 690 if(old_active[i]==0.){ 691 recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 692 } 693 } 744 694 } 745 695 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/ 746 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);747 696 if(eplhead[i]>=h_max && active_element){ 748 for( j=0;j<numnodes;j++){697 for(int j=0;j<numnodes;j++){ 749 698 /*Increase of the domain is on the downstream node in term of sediment head*/ 750 if( sedhead[j] == sedheadmin){699 if((sedhead[j] == sedheadmin) && (i!=j)){ 751 700 vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); 752 if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); 701 if(old_active[j]==0.){ 702 recurence->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); 703 } 753 704 } 754 705 } … … 763 714 xDelete<IssmDouble>(eplhead); 764 715 xDelete<IssmDouble>(residual); 716 xDelete<IssmDouble>(base); 765 717 } 766 718 /*}}}*/ 767 768 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){ 719 IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/ 720 IssmDouble epl_storing; 721 IssmDouble water_sheet,storing; 722 IssmDouble epl_thickness,prestep_head,base_elev; 723 IssmDouble rho_freshwater = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); 724 IssmDouble g = element->GetMaterialParameter(ConstantsGEnum); 725 IssmDouble epl_porosity = element->GetMaterialParameter(HydrologydcEplPorosityEnum); 726 IssmDouble epl_compressibility = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum); 727 IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum); 728 729 epl_thick_input->GetInputValue(&epl_thickness,gauss); 730 epl_head_input->GetInputValue(&prestep_head,gauss); 731 base_input->GetInputValue(&base_elev,gauss); 732 water_sheet=max(0.0,(prestep_head-base_elev)); 733 734 storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); 735 736 /* //porosity for unconfined region */ 737 /* if (water_sheet<=0.9*epl_thickness){ */ 738 /* epl_storing=epl_porosity; */ 739 /* } */ 740 /* //continuity ramp */ 741 /* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */ 742 /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */ 743 /* } */ 744 /* //storing coefficient for confined */ 745 /* else{ */ 746 /* epl_storing=storing; */ 747 /* } */ 748 /* return epl_storing; */ 749 return storing; 750 }/*}}}*/ 751 752 IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input* epl_thick_input, Input* epl_head_input, Input* base_input){/*{{{*/ 753 IssmDouble epl_transmitivity; 754 IssmDouble water_sheet; 755 IssmDouble epl_thickness,base_elev,prestep_head; 756 IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum); 757 epl_thick_input->GetInputValue(&epl_thickness,gauss); 758 epl_head_input->GetInputValue(&prestep_head,gauss); 759 base_input->GetInputValue(&base_elev,gauss); 760 761 water_sheet=max(0.0,(prestep_head-base_elev)); 762 763 epl_transmitivity=epl_conductivity*epl_thickness; 764 //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness))); 765 return epl_transmitivity; 766 }/*}}}*/ 767 768 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/ 769 769 /*Constants*/ 770 770 int domaintype; … … 787 787 IssmDouble flag = 0.; 788 788 IssmDouble* active = xNew<IssmDouble>(numnodes); 789 bool active_element; 789 790 790 791 /*Pass the activity mask from elements to nodes*/ 791 792 basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); 792 bool active_element;793 793 Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); 794 794 active_element_input->GetInputValue(&active_element); … … 814 814 xDelete<IssmDouble>(active); 815 815 } 816 /*}}}*/ 816 817 817 818 void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
Note:
See TracChangeset
for help on using the changeset viewer.
