Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21429) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21430) @@ -9,4 +9,5 @@ return 1; }/*}}}*/ + void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -33,7 +34,6 @@ iomodel->FetchData(&eplthickcomp,"md.hydrology.epl_thick_comp"); parameters->AddObject(new IntParam(HydrologydcEplThickCompEnum,eplthickcomp)); - - -}/*}}}*/ +}/*}}}*/ + void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -83,5 +83,7 @@ if(!isefficientlayer) return; - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + } ::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum); iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); @@ -106,12 +108,20 @@ /*Nothing for now*/ - /*Fetch parameters: */ + /*Do we really want DC?*/ int hydrology_model; iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); if(hydrology_model!=HydrologydcEnum) return; - if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase"); - - //create penalties for nodes: no node can have water above the max + /*Do we want an efficient layer*/ + bool isefficientlayer; + iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); + if(!isefficientlayer) return; + + /*Fetch parameters: */ + if(iomodel->domaintype==Domain3DEnum){ + iomodel->FetchData(1,"md.mesh.vertexonbase"); + } + + //Add moulin inputs as loads CreateSingleNodeToElementConnectivity(iomodel); for(int i=0;inumberofvertices;i++){ @@ -132,6 +142,6 @@ void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/ + int* eplzigzag_counter =NULL; - eplzigzag_counter=xNewZeroInit(femmodel->nodes->Size()); femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size())); @@ -142,9 +152,6 @@ int* eplzigzag_counter=NULL; - Element* element=NULL; - femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); for(int i=0;inodes->Size();i++){ - eplzigzag_counter[i]=0; } @@ -299,5 +306,5 @@ /*Check that all nodes are active, else return empty matrix*/ if(!active_element) { - if(domaintype!=Domain2DhorizontalEnum){ + if(domaintype!=Domain2DhorizontalEnum){ basalelement->DeleteMaterials(); delete basalelement; @@ -342,5 +349,4 @@ for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement ->NodalFunctions(basis,gauss); @@ -396,18 +402,22 @@ void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - bool active_element; - int domaintype,i; - Element* basalelement=NULL; - + /*Intermediaries*/ + bool active_element; + int domaintype; + Element* basalelement; + + /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); - - if(domaintype!=Domain2DhorizontalEnum){ - if(!element->IsOnBase()) return; - basalelement=element->SpawnBasalElement(); - } - else{ - basalelement = element; - } - + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<FindParam(&transfermethod,HydrologydcTransferFlagEnum); + + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just keepping the transfer to zero*/ + transfer=0.0; + break; + case 1: + element->FindParam(&leakage,HydrologydcLeakageFactorEnum); + transfer=leakage; + break; + default: + _error_("no case higher than 1 for the Transfer method"); + } + return transfer; +}/*}}}*/ + +IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ + + int transfermethod; + IssmDouble sediment_head; + IssmDouble leakage,transfer; element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); @@ -486,66 +516,7 @@ case 1: _assert_(sed_head_input); - _assert_(epl_head_input); - - inefanalysis = new HydrologyDCInefficientAnalysis(); - hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input); - delete inefanalysis; - sed_head_input->GetInputValue(&sediment_head,gauss); - epl_head_input->GetInputValue(&epl_head,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - //Computing continuum function to apply to transfer term, transfer is null only if - // epl_head>sediment_head AND sediment_head>h_max - //let's try without that for a while - /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */ - /* factor=min(continuum,1.0); */ - /* transfer=leakage*factor; */ - transfer=leakage; - break; - default: - _error_("no case higher than 1 for the Transfer method"); - } - return transfer; -}/*}}}*/ - -IssmDouble HydrologyDCEfficientAnalysis::GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input){/*{{{*/ - - int transfermethod; - IssmDouble hmax; - IssmDouble epl_head,sediment_head; - IssmDouble leakage,transfer; - IssmDouble continuum, factor; - - HydrologyDCInefficientAnalysis* inefanalysis = NULL; - - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - - /*Switch between the different transfer methods cases*/ - switch(transfermethod){ - case 0: - /*Just keepping the transfer to zero*/ - transfer=0.0; - break; - case 1: - _assert_(sed_head_input); - _assert_(epl_head_input); - - inefanalysis = new HydrologyDCInefficientAnalysis(); - hmax = inefanalysis->GetHydrologyDCInefficientHmax(element,gauss,thick_input,base_input); - delete inefanalysis; - - sed_head_input->GetInputValue(&sediment_head,gauss); - epl_head_input->GetInputValue(&epl_head,gauss); - element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - //Computing continuum function to apply to transfer term, transfer is null only if - // epl_head>sediment_head AND sediment_head>h_max - //let's try without that for a while - /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */ - /* factor=min(continuum,1.0); */ - /* transfer=sediment_head*leakage*factor; */ transfer=sediment_head*leakage; - break; default: @@ -564,5 +535,4 @@ IssmDouble EPLgrad2; IssmDouble EPL_N; - femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); @@ -619,5 +589,4 @@ if(!active_element){ - /*Keeping thickness to initial value if EPL is not active*/ for(int i=0;imax_thick){ @@ -695,5 +657,5 @@ IssmDouble h_max; IssmDouble sedheadmin; - Element* basalelement=NULL; + Element* basalelement; /*Get basal element*/ @@ -704,5 +666,5 @@ break; case Domain3DEnum: - if(!element->IsOnBase()) return; + if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; @@ -737,9 +699,10 @@ epl_thickness[i]=init_thick; } - /*Now starting to look at the activations*/ if(residual[i]>0.){ vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); - if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + if(old_active[i]==0.){ + recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } } /*If mask was already one, keep one or colapse*/ @@ -759,5 +722,7 @@ if(sedhead[j] == sedheadmin){ vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL); - if(old_active[i]==0.) recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + if(old_active[i]==0.){ + recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } } } @@ -778,5 +743,5 @@ /*Constants*/ int domaintype; - Element* basalelement=NULL; + Element* basalelement; /*Get basal element*/ @@ -787,5 +752,5 @@ break; case Domain3DEnum: - if(!element->IsOnBase()) return; + if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; @@ -796,8 +761,8 @@ IssmDouble flag = 0.; IssmDouble* active = xNew(numnodes); + bool active_element; /*Pass the activity mask from elements to nodes*/ basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum); - bool active_element; Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); active_element_input->GetInputValue(&active_element); Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21429) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21430) @@ -19,8 +19,8 @@ int hydro_maxiter; bool isefficientlayer; + IssmDouble penalty_factor; + IssmDouble rel_tol; + IssmDouble leakagefactor; IssmDouble sedimentlimit; - IssmDouble penalty_factor; - IssmDouble leakagefactor; - IssmDouble rel_tol; /*retrieve some parameters: */ @@ -30,30 +30,28 @@ if(hydrology_model!=HydrologydcEnum) return; - iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer"); iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" ); iomodel->FetchData(&transfer_flag, "md.hydrology.transfer_flag" ); iomodel->FetchData(&penalty_factor, "md.hydrology.penalty_factor" ); + iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer"); + iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" ); + iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" ); iomodel->FetchData(&rel_tol, "md.hydrology.rel_tol" ); - iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" ); - iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" ); - + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag)); + parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); + parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); + parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); + parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); + parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor)); + parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); + if(transfer_flag==1){ + iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor"); + parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor)); + } if(sedimentlimit_flag==1){ iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit"); parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit)); } - - if(transfer_flag==1){ - iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor"); - parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor)); - } - - parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor)); - parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); - parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); - parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag)); - parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); - parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); - parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); - parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter)); }/*}}}*/ @@ -109,5 +107,7 @@ if(hydrology_model!=HydrologydcEnum) return; - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + if(iomodel->domaintype!=Domain2DhorizontalEnum){ + iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); + } ::CreateNodes(nodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum); iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); @@ -131,6 +131,7 @@ if(hydrology_model!=HydrologydcEnum) return; - if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(1,"md.mesh.vertexonbase"); - + if(iomodel->domaintype==Domain3DEnum){ + iomodel->FetchData(1,"md.mesh.vertexonbase"); + } //create penalties for nodes: no node can have water above the max CreateSingleNodeToElementConnectivity(iomodel); @@ -190,5 +191,5 @@ IssmDouble D_scalar,Jdet,dt; IssmDouble sediment_transmitivity; - IssmDouble prestep_head,base_elev; + IssmDouble base_elev; IssmDouble transfer,sediment_storing; IssmDouble *xyz_list = NULL; @@ -301,5 +302,4 @@ IssmDouble *xyz_list = NULL; - Input* old_wh_input = NULL; Input* active_element_input = NULL; @@ -321,5 +321,4 @@ Input* base_input = basalelement->GetInput(BaseEnum); Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); - if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);} /*Transfer related Inputs*/ @@ -327,5 +326,4 @@ active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); } - /* Start looping on the number of gaussian points: */ @@ -333,5 +331,4 @@ for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); @@ -359,9 +356,7 @@ } - /*Transient and transfer terms*/ if(dt!=0.){ - old_wh_input ->GetInputValue(&water_head,gauss); - sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);//sed_head_input + sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input); if(isefficientlayer){ /*Dealing with the sediment part of the transfer term*/ @@ -428,17 +423,23 @@ void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - int domaintype; - bool converged; - int* doflist=NULL; - Element* basalelement=NULL; - + /*Intermediaries*/ + int domaintype; + Element* basalelement; + bool converged; + int* doflist = NULL; + + /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DhorizontalEnum){ - if(!element->IsOnBase()) return; - basalelement=element->SpawnBasalElement(); - } - else{ - basalelement = element; - } + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNode(i)); if(values[i]>h_max) { @@ -527,16 +527,20 @@ water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness))); specific_porosity=sediment_porosity-rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); - if (water_sheet<=0.9*sediment_thickness){//porosity for unconfined region + //porosity for unconfined region + if (water_sheet<=0.9*sediment_thickness){ sediment_storing=sediment_porosity; } - else if((water_sheet0.9*sediment_thickness)){//continuity ramp + //continuity ramp + else if((water_sheet0.9*sediment_thickness)){ sediment_storing=(sediment_thickness-water_sheet)*specific_porosity/(0.1*sediment_thickness)+ rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); } - else{//storing coefficient for confined + //storing coefficient for confined + else{ sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); } return sediment_storing; }/*}}}*/ + IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input){/*{{{*/ IssmDouble sediment_transmitivity; @@ -574,11 +578,8 @@ break; case 2: - rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - _assert_(thick_input); _assert_(base_input); - /*Compute max*/ thick_input->GetInputValue(&thickness,gauss); @@ -634,9 +635,5 @@ int transfermethod; - IssmDouble hmax; - IssmDouble epl_head,sediment_head; IssmDouble leakage,transfer; - IssmDouble continuum, factor; - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); @@ -648,21 +645,6 @@ break; case 1: - _assert_(sed_head_input); - _assert_(epl_head_input); - - sed_head_input->GetInputValue(&sediment_head,gauss); - epl_head_input->GetInputValue(&epl_head,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input); - - //Computing continuum function to apply to transfer term, transfer is null only if - //epl_head>sediment_head AND sediment_head>h_max - //let's try without that for a while - /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */ - /* factor=min(continuum,1.0); */ - /* transfer=leakage*factor; */ transfer=leakage; - break; default: @@ -675,9 +657,6 @@ int transfermethod; - IssmDouble hmax; - IssmDouble epl_head,sediment_head; + IssmDouble epl_head; IssmDouble leakage,transfer; - IssmDouble continuum, factor; - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); @@ -689,22 +668,8 @@ break; case 1: - _assert_(sed_head_input); _assert_(epl_head_input); - - sed_head_input->GetInputValue(&sediment_head,gauss); epl_head_input->GetInputValue(&epl_head,gauss); element->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - hmax=GetHydrologyDCInefficientHmax(element, gauss, thick_input, base_input); - - //Computing continuum function to apply to transfer term, transfer is null only if - //epl_head>sediment_head AND sediment_head>h_max - //let's try without that for a while - /* continuum=((1.0/(1.0+exp(-20.0*(sediment_head-epl_head)))))+(1.0/(1.0+exp(-20.0*(hmax-sediment_head)))); */ - /* factor=min(continuum,1.0); */ - /* transfer=epl_head*leakage*factor; */ - transfer=epl_head*leakage; - break; default: @@ -731,3 +696,4 @@ element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); } -}/*}}}*/ + delete element; +}/*}}}*/