Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 24384) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 24385) @@ -383,4 +383,31 @@ return pe; }/*}}}*/ +void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For node i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dN/dx ] + * [ dN/dy ] + * where N is the finiteelement function for node i. + * + * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(2*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +}/*}}}*/ void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum); @@ -437,4 +464,84 @@ /*Intermediaries*/ +IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ + IssmDouble epl_storing; + IssmDouble water_sheet,storing; + IssmDouble epl_thickness,prestep_head,base_elev; + IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); + IssmDouble epl_porosity = element->FindParam(HydrologydcEplPorosityEnum); + IssmDouble epl_compressibility = element->FindParam(HydrologydcEplCompressibilityEnum); + IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum); + + epl_thick_input->GetInputValue(&epl_thickness,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + base_input->GetInputValue(&base_elev,gauss); + water_sheet=max(0.0,(prestep_head-base_elev)); + storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); + + /* //porosity for unconfined region */ + /* if (water_sheet<=0.9*epl_thickness){ */ + /* epl_storing=epl_porosity; */ + /* } */ + /* //continuity ramp */ + /* else if((water_sheet0.9*epl_thickness)){ */ + /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */ + /* } */ + /* //storing coefficient for confined */ + /* else{ */ + /* epl_storing=storing; */ + /* } */ + /* return epl_storing; */ + return storing; +}/*}}}*/ +IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ + IssmDouble epl_transmitivity; + IssmDouble water_sheet; + IssmDouble epl_thickness,base_elev,prestep_head; + IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum); + epl_thick_input->GetInputValue(&epl_thickness,gauss); + epl_head_input->GetInputValue(&prestep_head,gauss); + base_input->GetInputValue(&base_elev,gauss); + + water_sheet=max(0.0,(prestep_head-base_elev)); + epl_transmitivity=epl_conductivity*epl_thickness; + //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness))); + return epl_transmitivity; +}/*}}}*/ +void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ + + int hmax_flag; + IssmDouble h_max; + IssmDouble rho_ice,rho_water; + IssmDouble thickness,bed; + /*Get the flag to the limitation method*/ + element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); + + /*Switch between the different cases*/ + switch(hmax_flag){ + case 0: + h_max=1.0e+10; + break; + case 1: + element->FindParam(&h_max,HydrologydcSedimentlimitEnum); + break; + case 2: + /*Compute max*/ + rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + rho_ice = element->FindParam(MaterialsRhoIceEnum); + element-> GetInputValue(&thickness,innode,ThicknessEnum); + element-> GetInputValue(&bed,innode,BaseEnum); + h_max=((rho_ice*thickness)/rho_water)+bed; + break; + case 3: + _error_("Using normal stress not supported yet"); + break; + default: + _error_("no case higher than 3 for SedimentlimitFlag"); + } + /*Assign output pointer*/ + *ph_max=h_max; +} +/*}}}*/ IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyKMatrixTransfer(Element* element){/*{{{*/ @@ -582,31 +689,4 @@ } }/*}}}*/ -void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. - * For node i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ dN/dx ] - * [ dN/dy ] - * where N is the finiteelement function for node i. - * - * We assume B has been allocated already, of size: 3x(NDOF2*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - - /*Build B: */ - for(int i=0;i(dbasis); -}/*}}}*/ void HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector* vec_mask, Vector* recurence, Element* element){/*{{{*/ @@ -703,49 +783,4 @@ } /*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ - IssmDouble epl_storing; - IssmDouble water_sheet,storing; - IssmDouble epl_thickness,prestep_head,base_elev; - IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum); - IssmDouble g = element->FindParam(ConstantsGEnum); - IssmDouble epl_porosity = element->FindParam(HydrologydcEplPorosityEnum); - IssmDouble epl_compressibility = element->FindParam(HydrologydcEplCompressibilityEnum); - IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum); - - epl_thick_input->GetInputValue(&epl_thickness,gauss); - epl_head_input->GetInputValue(&prestep_head,gauss); - base_input->GetInputValue(&base_elev,gauss); - water_sheet=max(0.0,(prestep_head-base_elev)); - storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); - - /* //porosity for unconfined region */ - /* if (water_sheet<=0.9*epl_thickness){ */ - /* epl_storing=epl_porosity; */ - /* } */ - /* //continuity ramp */ - /* else if((water_sheet0.9*epl_thickness)){ */ - /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */ - /* } */ - /* //storing coefficient for confined */ - /* else{ */ - /* epl_storing=storing; */ - /* } */ - /* return epl_storing; */ - return storing; -}/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ - IssmDouble epl_transmitivity; - IssmDouble water_sheet; - IssmDouble epl_thickness,base_elev,prestep_head; - IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum); - epl_thick_input->GetInputValue(&epl_thickness,gauss); - epl_head_input->GetInputValue(&prestep_head,gauss); - base_input->GetInputValue(&base_elev,gauss); - - water_sheet=max(0.0,(prestep_head-base_elev)); - epl_transmitivity=epl_conductivity*epl_thickness; - //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness))); - return epl_transmitivity; -}/*}}}*/ void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector* active_vec, Element* element){/*{{{*/ /*Constants*/ @@ -796,37 +831,2 @@ } /*}}}*/ -void HydrologyDCEfficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/ - - int hmax_flag; - IssmDouble h_max; - IssmDouble rho_ice,rho_water; - IssmDouble thickness,bed; - /*Get the flag to the limitation method*/ - element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum); - - /*Switch between the different cases*/ - switch(hmax_flag){ - case 0: - h_max=1.0e+10; - break; - case 1: - element->FindParam(&h_max,HydrologydcSedimentlimitEnum); - break; - case 2: - /*Compute max*/ - rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); - rho_ice = element->FindParam(MaterialsRhoIceEnum); - element-> GetInputValue(&thickness,innode,ThicknessEnum); - element-> GetInputValue(&bed,innode,BaseEnum); - h_max=((rho_ice*thickness)/rho_water)+bed; - break; - case 3: - _error_("Using normal stress not supported yet"); - break; - default: - _error_("no case higher than 3 for SedimentlimitFlag"); - } - /*Assign output pointer*/ - *ph_max=h_max; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 24384) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 24385) @@ -29,4 +29,5 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); + void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); @@ -35,13 +36,12 @@ /*Intermediaries*/ - void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + 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); - 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 ComputeEPLThickness(FemModel* femmodel); void HydrologyEPLGetMask(Vector* vec_mask, Vector* recurence, Element* element); void HydrologyEPLGetActive(Vector* active_vec, Element* element); - void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); - void ComputeEPLThickness(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24384) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 24385) @@ -587,4 +587,5 @@ return; }/*}}}*/ +/*Intermediaries*/ IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input){/*{{{*/ int unconf_scheme; @@ -776,5 +777,5 @@ } /*Intermediaries*/ - int numnodes = basalelement->GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); IssmDouble* meltingrate = xNew(numnodes); IssmDouble* groundedice = xNew(numnodes); Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24384) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24385) @@ -153,8 +153,5 @@ break; case SMBgradientscomponentsEnum: - iomodel->FetchDataToInput(inputs2,elements,"md.smb.accualti",SmbAccualtiEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.accugrad",SmbAccugradEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.runoffalti",SmbRunoffaltiEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.runoffgrad",SmbRunoffgradEnum); + /* Nothing to add to input */ break; case SMBsemicEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24384) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24385) @@ -664,10 +664,10 @@ IssmDouble accugrad, runoffgrad; //gradients from reference altitude IssmDouble rho_water, rho_ice; - IssmDouble time,yts; - - IssmDouble* smb = xNew(NUM_VERTICES); - IssmDouble* surf = xNew(NUM_VERTICES); - IssmDouble* accu = xNew(NUM_VERTICES); - IssmDouble* runoff = xNew(NUM_VERTICES); + IssmDouble time; + + IssmDouble* smb = xNew(NUM_VERTICES); + IssmDouble* surf = xNew(NUM_VERTICES); + IssmDouble* accu = xNew(NUM_VERTICES); + IssmDouble* runoff = xNew(NUM_VERTICES); /*Get material parameters :*/ @@ -677,5 +677,4 @@ /*Recover parameters*/ parameters->FindParam(&time,TimeEnum); - parameters->FindParam(&yts,ConstantsYtsEnum); parameters->FindParam(&accualti,SmbAccualtiEnum); parameters->FindParam(&accugrad,SmbAccugradEnum); Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 24384) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 24385) @@ -80,5 +80,5 @@ Var = md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] DimDict = CreateVar(NCData, Var, subfield, Listgroup, DimDict, md.__dict__[group], field, listindex) - # No subgroup, we directly treat the variable + # No subgroup, we directly treat the variable elif type(md.__dict__[group].__dict__[field]) in typelist or field == 'bamg': NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) @@ -87,6 +87,6 @@ elif md.__dict__[group].__dict__[field] is None: print('field md.{}.{} is None'.format(group, field)) - # do nothing - # if it is a masked array + # do nothing + # if it is a masked array elif type(md.__dict__[group].__dict__[field]) is np.ma.core.MaskedArray: NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__) @@ -175,8 +175,8 @@ print(('WARNING type "{}" is unknown for "{}.{}"'.format(val_type, Group.name, field))) return DimDict + + # ============================================================================ - # retriev the dimension tuple from a dictionnary - - +# retriev the dimension tuple from a dictionnary def GetDim(NCData, val_shape, val_type, DimDict, val_dim): output = [] Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24384) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 24385) @@ -94,4 +94,5 @@ if enveloppe: if dim == 3: + mesh_alti = '1' is_enveloppe = np.logical_or(md.mesh.vertexonbase, md.mesh.vertexonsurface) enveloppe_index = np.where(is_enveloppe)[0]