![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 24385
- Timestamp:
- 11/22/19 02:51:08 (5 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 7 edited
-
c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (5 diffs)
-
c/analyses/HydrologyDCEfficientAnalysis.h (modified) (2 diffs)
-
c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (2 diffs)
-
c/analyses/SmbAnalysis.cpp (modified) (1 diff)
-
c/classes/Elements/Element.cpp (modified) (2 diffs)
-
m/contrib/defleurian/netCDF/export_netCDF.py (modified) (3 diffs)
-
m/contrib/defleurian/paraview/exportVTK.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
r24344 r24385 383 383 return pe; 384 384 }/*}}}*/ 385 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 386 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 387 * For node i, Bi can be expressed in the actual coordinate system 388 * by: 389 * Bi=[ dN/dx ] 390 * [ dN/dy ] 391 * where N is the finiteelement function for node i. 392 * 393 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) 394 */ 395 396 /*Fetch number of nodes for this finite element*/ 397 int numnodes = element->GetNumberOfNodes(); 398 399 /*Get nodal functions derivatives*/ 400 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes); 401 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); 402 403 /*Build B: */ 404 for(int i=0;i<numnodes;i++){ 405 B[numnodes*0+i] = dbasis[0*numnodes+i]; 406 B[numnodes*1+i] = dbasis[1*numnodes+i]; 407 } 408 409 /*Clean-up*/ 410 xDelete<IssmDouble>(dbasis); 411 }/*}}}*/ 385 412 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/ 386 413 element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum); … … 437 464 438 465 /*Intermediaries*/ 466 IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ 467 IssmDouble epl_storing; 468 IssmDouble water_sheet,storing; 469 IssmDouble epl_thickness,prestep_head,base_elev; 470 IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum); 471 IssmDouble g = element->FindParam(ConstantsGEnum); 472 IssmDouble epl_porosity = element->FindParam(HydrologydcEplPorosityEnum); 473 IssmDouble epl_compressibility = element->FindParam(HydrologydcEplCompressibilityEnum); 474 IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum); 475 476 epl_thick_input->GetInputValue(&epl_thickness,gauss); 477 epl_head_input->GetInputValue(&prestep_head,gauss); 478 base_input->GetInputValue(&base_elev,gauss); 479 water_sheet=max(0.0,(prestep_head-base_elev)); 480 storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); 481 482 /* //porosity for unconfined region */ 483 /* if (water_sheet<=0.9*epl_thickness){ */ 484 /* epl_storing=epl_porosity; */ 485 /* } */ 486 /* //continuity ramp */ 487 /* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */ 488 /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */ 489 /* } */ 490 /* //storing coefficient for confined */ 491 /* else{ */ 492 /* epl_storing=storing; */ 493 /* } */ 494 /* return epl_storing; */ 495 return storing; 496 }/*}}}*/ 497 IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ 498 IssmDouble epl_transmitivity; 499 IssmDouble water_sheet; 500 IssmDouble epl_thickness,base_elev,prestep_head; 501 IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum); 502 epl_thick_input->GetInputValue(&epl_thickness,gauss); 503 epl_head_input->GetInputValue(&prestep_head,gauss); 504 base_input->GetInputValue(&base_elev,gauss); 505 506 water_sheet=max(0.0,(prestep_head-base_elev)); 507 epl_transmitivity=epl_conductivity*epl_thickness; 508 //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness))); 509 return epl_transmitivity; 510 }/*}}}*/ 511 void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ 512 513 int hmax_flag; 514 IssmDouble h_max; 515 IssmDouble rho_ice,rho_water; 516 IssmDouble thickness,bed; 517 /*Get the flag to the limitation method*/ 518 element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); 519 520 /*Switch between the different cases*/ 521 switch(hmax_flag){ 522 case 0: 523 h_max=1.0e+10; 524 break; 525 case 1: 526 element->FindParam(&h_max,HydrologydcSedimentlimitEnum); 527 break; 528 case 2: 529 /*Compute max*/ 530 rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); 531 rho_ice = element->FindParam(MaterialsRhoIceEnum); 532 element-> GetInputValue(&thickness,innode,ThicknessEnum); 533 element-> GetInputValue(&bed,innode,BaseEnum); 534 h_max=((rho_ice*thickness)/rho_water)+bed; 535 break; 536 case 3: 537 _error_("Using normal stress not supported yet"); 538 break; 539 default: 540 _error_("no case higher than 3 for SedimentlimitFlag"); 541 } 542 /*Assign output pointer*/ 543 *ph_max=h_max; 544 } 545 /*}}}*/ 439 546 IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/ 440 547 … … 582 689 } 583 690 }/*}}}*/ 584 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/585 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.586 * For node i, Bi can be expressed in the actual coordinate system587 * by:588 * Bi=[ dN/dx ]589 * [ dN/dy ]590 * where N is the finiteelement function for node i.591 *592 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)593 */594 595 /*Fetch number of nodes for this finite element*/596 int numnodes = element->GetNumberOfNodes();597 598 /*Get nodal functions derivatives*/599 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);600 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);601 602 /*Build B: */603 for(int i=0;i<numnodes;i++){604 B[numnodes*0+i] = dbasis[0*numnodes+i];605 B[numnodes*1+i] = dbasis[1*numnodes+i];606 }607 608 /*Clean-up*/609 xDelete<IssmDouble>(dbasis);610 }/*}}}*/611 691 void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, Element* element){/*{{{*/ 612 692 … … 703 783 } 704 784 /*}}}*/ 705 IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/706 IssmDouble epl_storing;707 IssmDouble water_sheet,storing;708 IssmDouble epl_thickness,prestep_head,base_elev;709 IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum);710 IssmDouble g = element->FindParam(ConstantsGEnum);711 IssmDouble epl_porosity = element->FindParam(HydrologydcEplPorosityEnum);712 IssmDouble epl_compressibility = element->FindParam(HydrologydcEplCompressibilityEnum);713 IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);714 715 epl_thick_input->GetInputValue(&epl_thickness,gauss);716 epl_head_input->GetInputValue(&prestep_head,gauss);717 base_input->GetInputValue(&base_elev,gauss);718 water_sheet=max(0.0,(prestep_head-base_elev));719 storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));720 721 /* //porosity for unconfined region */722 /* if (water_sheet<=0.9*epl_thickness){ */723 /* epl_storing=epl_porosity; */724 /* } */725 /* //continuity ramp */726 /* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */727 /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */728 /* } */729 /* //storing coefficient for confined */730 /* else{ */731 /* epl_storing=storing; */732 /* } */733 /* return epl_storing; */734 return storing;735 }/*}}}*/736 IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/737 IssmDouble epl_transmitivity;738 IssmDouble water_sheet;739 IssmDouble epl_thickness,base_elev,prestep_head;740 IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum);741 epl_thick_input->GetInputValue(&epl_thickness,gauss);742 epl_head_input->GetInputValue(&prestep_head,gauss);743 base_input->GetInputValue(&base_elev,gauss);744 745 water_sheet=max(0.0,(prestep_head-base_elev));746 epl_transmitivity=epl_conductivity*epl_thickness;747 //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));748 return epl_transmitivity;749 }/*}}}*/750 785 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){/*{{{*/ 751 786 /*Constants*/ … … 796 831 } 797 832 /*}}}*/ 798 void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/799 800 int hmax_flag;801 IssmDouble h_max;802 IssmDouble rho_ice,rho_water;803 IssmDouble thickness,bed;804 /*Get the flag to the limitation method*/805 element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);806 807 /*Switch between the different cases*/808 switch(hmax_flag){809 case 0:810 h_max=1.0e+10;811 break;812 case 1:813 element->FindParam(&h_max,HydrologydcSedimentlimitEnum);814 break;815 case 2:816 /*Compute max*/817 rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);818 rho_ice = element->FindParam(MaterialsRhoIceEnum);819 element-> GetInputValue(&thickness,innode,ThicknessEnum);820 element-> GetInputValue(&bed,innode,BaseEnum);821 h_max=((rho_ice*thickness)/rho_water)+bed;822 break;823 case 3:824 _error_("Using normal stress not supported yet");825 break;826 default:827 _error_("no case higher than 3 for SedimentlimitFlag");828 }829 /*Assign output pointer*/830 *ph_max=h_max;831 }832 /*}}}*/ -
issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h
r24335 r24385 29 29 ElementMatrix* CreateKMatrix(Element* element); 30 30 ElementVector* CreatePVector(Element* element); 31 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); 31 32 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element); 32 33 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index); … … 35 36 36 37 /*Intermediaries*/ 37 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); 38 IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); 39 IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); 40 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); 38 41 IssmDouble GetHydrologyKMatrixTransfer(Element* element); 39 42 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input2* sed_head_input); 40 IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); 41 IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); 43 void ComputeEPLThickness(FemModel* femmodel); 42 44 void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, Element* element); 43 45 void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element); 44 void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);45 void ComputeEPLThickness(FemModel* femmodel);46 46 }; 47 47 #endif -
issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
r24335 r24385 587 587 return; 588 588 }/*}}}*/ 589 /*Intermediaries*/ 589 590 IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input){/*{{{*/ 590 591 int unconf_scheme; … … 776 777 } 777 778 /*Intermediaries*/ 778 int 779 int numnodes = basalelement->GetNumberOfNodes(); 779 780 IssmDouble* meltingrate = xNew<IssmDouble>(numnodes); 780 781 IssmDouble* groundedice = xNew<IssmDouble>(numnodes); -
issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp
r24360 r24385 153 153 break; 154 154 case SMBgradientscomponentsEnum: 155 iomodel->FetchDataToInput(inputs2,elements,"md.smb.accualti",SmbAccualtiEnum); 156 iomodel->FetchDataToInput(inputs2,elements,"md.smb.accugrad",SmbAccugradEnum); 157 iomodel->FetchDataToInput(inputs2,elements,"md.smb.runoffalti",SmbRunoffaltiEnum); 158 iomodel->FetchDataToInput(inputs2,elements,"md.smb.runoffgrad",SmbRunoffgradEnum); 155 /* Nothing to add to input */ 159 156 break; 160 157 case SMBsemicEnum: -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r24379 r24385 664 664 IssmDouble accugrad, runoffgrad; //gradients from reference altitude 665 665 IssmDouble rho_water, rho_ice; 666 IssmDouble time ,yts;667 668 IssmDouble* smb 669 IssmDouble* surf = xNew<IssmDouble>(NUM_VERTICES);670 IssmDouble* accu = xNew<IssmDouble>(NUM_VERTICES);671 IssmDouble* runoff 666 IssmDouble time; 667 668 IssmDouble* smb = xNew<IssmDouble>(NUM_VERTICES); 669 IssmDouble* surf = xNew<IssmDouble>(NUM_VERTICES); 670 IssmDouble* accu = xNew<IssmDouble>(NUM_VERTICES); 671 IssmDouble* runoff = xNew<IssmDouble>(NUM_VERTICES); 672 672 673 673 /*Get material parameters :*/ … … 677 677 /*Recover parameters*/ 678 678 parameters->FindParam(&time,TimeEnum); 679 parameters->FindParam(&yts,ConstantsYtsEnum);680 679 parameters->FindParam(&accualti,SmbAccualtiEnum); 681 680 parameters->FindParam(&accugrad,SmbAccugradEnum); -
issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py
r24241 r24385 80 80 Var = md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] 81 81 DimDict = CreateVar(NCData, Var, subfield, Listgroup, DimDict, md.__dict__[group], field, listindex) 82 # No subgroup, we directly treat the variable82 # No subgroup, we directly treat the variable 83 83 elif type(md.__dict__[group].__dict__[field]) in typelist or field == 'bamg': 84 84 NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) … … 87 87 elif md.__dict__[group].__dict__[field] is None: 88 88 print('field md.{}.{} is None'.format(group, field)) 89 # do nothing90 # if it is a masked array89 # do nothing 90 # if it is a masked array 91 91 elif type(md.__dict__[group].__dict__[field]) is np.ma.core.MaskedArray: 92 92 NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) … … 175 175 print(('WARNING type "{}" is unknown for "{}.{}"'.format(val_type, Group.name, field))) 176 176 return DimDict 177 178 177 179 # ============================================================================ 178 # retriev the dimension tuple from a dictionnary 179 180 180 # retriev the dimension tuple from a dictionnary 181 181 def GetDim(NCData, val_shape, val_type, DimDict, val_dim): 182 182 output = [] -
issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py
r24275 r24385 94 94 if enveloppe: 95 95 if dim == 3: 96 mesh_alti = '1' 96 97 is_enveloppe = np.logical_or(md.mesh.vertexonbase, md.mesh.vertexonsurface) 97 98 enveloppe_index = np.where(is_enveloppe)[0]
Note:
See TracChangeset
for help on using the changeset viewer.
