Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21429) +++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21430) @@ -18,10 +18,10 @@ int penalty_lock; int hydro_maxiter; bool isefficientlayer; - IssmDouble sedimentlimit; IssmDouble penalty_factor; - IssmDouble leakagefactor; IssmDouble rel_tol; + IssmDouble leakagefactor; + IssmDouble sedimentlimit; /*retrieve some parameters: */ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); @@ -29,32 +29,30 @@ /*Now, do we really want DC?*/ 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(&rel_tol, "md.hydrology.rel_tol" ); - iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" ); + 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" ); - 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)); + 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)); + } }/*}}}*/ void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -108,7 +106,9 @@ /*Now, do we really want DC?*/ 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"); }/*}}}*/ @@ -130,8 +130,9 @@ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); 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); for(int i=0;inumberofvertices;i++){ @@ -189,7 +190,7 @@ bool active_element,isefficientlayer; IssmDouble D_scalar,Jdet,dt; IssmDouble sediment_transmitivity; - IssmDouble prestep_head,base_elev; + IssmDouble base_elev; IssmDouble transfer,sediment_storing; IssmDouble *xyz_list = NULL; @@ -300,7 +301,6 @@ IssmDouble Jdet; IssmDouble *xyz_list = NULL; - Input* old_wh_input = NULL; Input* active_element_input = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -320,19 +320,16 @@ Input* thick_input = basalelement->GetInput(ThicknessEnum); 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*/ if(isefficientlayer){ active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); } - /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); @@ -358,11 +355,9 @@ } } - /*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*/ active_element_input->GetInputValue(&active_element); @@ -427,20 +422,26 @@ 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(); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); @@ -478,7 +479,6 @@ kappa=kmax*pow(10.,penalty_factor); for(int i=0;iGetNode(i)); if(values[i]>h_max) { residual[i] = kappa*(values[i]-h_max); @@ -526,18 +526,22 @@ sed_head_input->GetInputValue(&prestep_head,gauss); 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; IssmDouble FullLayer_transmitivity; @@ -573,13 +577,10 @@ element->FindParam(&h_max,HydrologydcSedimentlimitEnum); 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); base_input->GetInputValue(&bed,gauss); @@ -633,11 +634,7 @@ IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer(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; - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); /*Switch between the different transfer methods cases*/ @@ -647,23 +644,8 @@ transfer=0.0; 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: _error_("no case higher than 1 for the Transfer method"); @@ -674,11 +656,8 @@ IssmDouble HydrologyDCInefficientAnalysis::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 epl_head; IssmDouble leakage,transfer; - IssmDouble continuum, factor; - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); /*Switch between the different transfer methods cases*/ @@ -688,24 +667,10 @@ transfer=0.0; 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: _error_("no case higher than 1 for the Transfer method"); @@ -730,4 +695,5 @@ } element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); } + delete element; }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21429) +++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 21430) @@ -8,6 +8,7 @@ int HydrologyDCEfficientAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ return 1; }/*}}}*/ + void HydrologyDCEfficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ int hydrology_model; @@ -32,9 +33,8 @@ 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){/*{{{*/ bool isefficientlayer; @@ -82,7 +82,9 @@ iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); 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"); }/*}}}*/ @@ -105,14 +107,22 @@ void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ /*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"); + /*Do we want an efficient layer*/ + bool isefficientlayer; + iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); + if(!isefficientlayer) return; - //create penalties for nodes: no node can have water above the max + /*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++){ if (iomodel->domaintype!=Domain3DEnum){ @@ -131,8 +141,8 @@ }/*}}}*/ 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())); xDelete(eplzigzag_counter); @@ -141,11 +151,8 @@ void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/ 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; } femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum); @@ -298,7 +305,7 @@ /*Check that all nodes are active, else return empty matrix*/ if(!active_element) { - if(domaintype!=Domain2DhorizontalEnum){ + if(domaintype!=Domain2DhorizontalEnum){ basalelement->DeleteMaterials(); delete basalelement; } @@ -341,7 +348,6 @@ Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement ->NodalFunctions(basis,gauss); @@ -395,20 +401,24 @@ 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(); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<FindParam(&transfermethod,HydrologydcTransferFlagEnum); @@ -484,23 +490,7 @@ 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=leakage*factor; */ transfer=leakage; break; default: @@ -512,13 +502,9 @@ 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 sediment_head; IssmDouble leakage,transfer; - IssmDouble continuum, factor; - HydrologyDCInefficientAnalysis* inefanalysis = NULL; - element->FindParam(&transfermethod,HydrologydcTransferFlagEnum); /*Switch between the different transfer methods cases*/ @@ -529,24 +515,9 @@ 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: _error_("no case higher than 1 for the Transfer method"); @@ -563,7 +534,6 @@ IssmDouble dt,A,B; IssmDouble EPLgrad2; IssmDouble EPL_N; - femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); @@ -618,7 +588,6 @@ element->GetInputListOnVertices(&bed[0],BaseEnum); if(!active_element){ - /*Keeping thickness to initial value if EPL is not active*/ for(int i=0;imax_thick){ thickness[i] = max_thick; @@ -694,7 +656,7 @@ int domaintype; IssmDouble h_max; IssmDouble sedheadmin; - Element* basalelement=NULL; + Element* basalelement; /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); @@ -703,7 +665,7 @@ basalelement = element; break; case Domain3DEnum: - if(!element->IsOnBase()) return; + if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<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*/ else if(old_active[i]>0.){ @@ -758,7 +721,9 @@ /*Increase of the domain is on the downstream node in term of sediment head*/ 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); + } } } } @@ -777,7 +742,7 @@ void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector* active_vec, Element* element){ /*Constants*/ int domaintype; - Element* basalelement=NULL; + Element* basalelement; /*Get basal element*/ element->FindParam(&domaintype,DomainTypeEnum); @@ -786,7 +751,7 @@ basalelement = element; break; case Domain3DEnum: - if(!element->IsOnBase()) return; + if(!element->IsOnBase()) return NULL; basalelement = element->SpawnBasalElement(); break; default: _error_("mesh "<GetNumberOfNodes(); 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);