![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 17024
- Timestamp:
- 12/16/13 15:45:23 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 edited
-
analyses/HydrologyDCEfficientAnalysis.cpp (modified) (3 diffs)
-
solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
r17021 r17024 128 128 129 129 /*Check that all nodes are active, else return empty matrix*/ 130 if(!basalelement->AllActive()) return NULL; 131 130 if(!basalelement->AllActive()) { 131 if(meshtype!=Mesh2DhorizontalEnum){ 132 basalelement->DeleteMaterials(); 133 delete basalelement; 134 } 135 return NULL; 136 } 132 137 /* Intermediaries */ 133 138 IssmDouble D_scalar,Jdet,dt; … … 186 191 xDelete<IssmDouble>(B); 187 192 delete gauss; 193 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 188 194 return Ke; 189 195 … … 209 215 210 216 /*Check that all nodes are active, else return empty matrix*/ 211 if(!basalelement->AllActive()) return NULL; 212 217 if(!basalelement->AllActive()){ 218 if(meshtype!=Mesh2DhorizontalEnum){ 219 basalelement->DeleteMaterials(); 220 delete basalelement; 221 } 222 return NULL; 223 } 213 224 /*Intermediaries */ 214 225 IssmDouble dt,scalar,water_head,connectivity; -
issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
r17021 r17024 23 23 Vector<IssmDouble>* ug_epl_main_iter=NULL; 24 24 25 26 25 Vector<IssmDouble>* ys=NULL; 27 26 Vector<IssmDouble>* dug=NULL; 28 29 //testing stuff30 27 Vector<IssmDouble>* duf=NULL; 31 28 32 29 Matrix<IssmDouble>* Kff=NULL; 33 30 Matrix<IssmDouble>* Kfs=NULL; 31 34 32 Vector<IssmDouble>* pf=NULL; 35 33 Vector<IssmDouble>* df=NULL; … … 51 49 femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum); 52 50 femmodel->parameters->FindParam(&time,TimeEnum); 51 53 52 /*FIXME, hardcoded, put on an enum*/ 54 53 hydro_maxiter=150; 54 55 55 hydrocount=1; 56 56 hydroconverged=false; 57 /*We don't need the outer loop if only one layer is used*/ 58 if(!isefficientlayer) hydroconverged=true; 57 59 58 60 /*Retrieve inputs as the initial state for the non linear iteration*/ … … 63 65 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); 64 66 GetSolutionFromInputsx(&ug_epl,femmodel); 65 66 /*Initialize the transfer input*/67 68 /*Initialize the element mask*/ 67 69 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); 68 70 analysis->ElementizeEplMask(femmodel); … … 72 74 femmodel->HydrologyTransferx(); 73 75 74 /* Iteration on the two layers*/76 /*The real computation starts here, outermost loop is on the two layer system*/ 75 77 for(;;){ 76 78 sedcount=1; 77 79 eplcount=1; 78 //save pointer to old velocity 79 ug_sed_main_iter=ug_sed->Duplicate(); 80 ug_sed->Copy(ug_sed_main_iter); 81 80 81 /*If there is two layers we need an outer loop value to compute convergence*/ 82 82 if(isefficientlayer){ 83 ug_sed_main_iter=ug_sed->Duplicate(); 84 ug_sed->Copy(ug_sed_main_iter); 83 85 ug_epl_main_iter=ug_epl->Duplicate(); 84 86 ug_epl->Copy(ug_epl_main_iter); 85 87 } 86 88 /*Loop on sediment layer to deal with transfer and head value*/ 87 89 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 88 90 InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); … … 90 92 femmodel->UpdateConstraintsx(); 91 93 femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); 92 93 /*Reset constraint on the ZigZag Lock */94 95 /*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/ 94 96 ResetConstraintsx(femmodel); 95 96 /*Iteration on the sediment layer*/97 97 sedconverged=false; 98 /* {{{ *//*Treating the sediment*/ 98 99 for(;;){ 99 /* if(isefficientlayer){ */100 /* /\*Updating Elemental Mask*\/ */101 /* HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */102 /* analysis->ElementizeEplMask(femmodel); */103 /* delete analysis; */104 /* femmodel->HydrologyTransferx(); */105 /* } */106 100 uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter); 107 101 uf_sed->Copy(uf_sed_sub_iter); 108 SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); 109 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); 110 Reduceloadx(pf,Kfs,ys); delete Kfs; 111 delete uf_sed; 112 Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); 113 delete Kff; delete pf; delete df; 114 delete ug_sed; 115 Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; 116 InputUpdateFromSolutionx(femmodel,ug_sed); 117 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); 118 119 if (!sedconverged){ 120 if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); 121 if(num_unstable_constraints==0) sedconverged = true; 122 if (sedcount>=hydro_maxiter){ 123 //sedconverged=true; 124 _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); 125 } 126 } 127 128 sedcount++; 129 102 /* {{{ *//*Loop on the sediment layer to deal with the penalization*/ 103 for(;;){ 104 /* {{{ *//*Core of the computation*/ 105 SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel); 106 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); 107 Reduceloadx(pf,Kfs,ys); delete Kfs; 108 delete uf_sed; 109 Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters); 110 delete Kff; delete pf; delete df; 111 delete ug_sed; 112 Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys; 113 InputUpdateFromSolutionx(femmodel,ug_sed); 114 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); 115 /* }}} */ 116 if (!sedconverged){ 117 if(VerboseConvergence()) _printf0_(" # Sediment unstable constraints = " << num_unstable_constraints << "\n"); 118 if(num_unstable_constraints==0) sedconverged = true; 119 if (sedcount>=hydro_maxiter){ 120 _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); 121 } 122 } 123 /*Add an iteration and get out of the loop if the penalisation is converged*/ 124 sedcount++; 125 if(sedconverged)break; 126 } 127 /* }}} *//*End of the sediment penalization loop*/ 128 sedconverged=false; 129 /*Update Elemental Mask and compute transfer*/ 130 if(isefficientlayer){ 131 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); 132 analysis->ElementizeEplMask(femmodel); 133 delete analysis; 134 femmodel->HydrologyTransferx(); 135 } 136 /*Checking convegence on the value of the sediment head*/ 137 duf=uf_sed_sub_iter->Duplicate();_assert_(duf); 138 uf_sed_sub_iter->Copy(duf); 139 duf->AYPX(uf_sed,-1.0); 140 ndu_sed=duf->Norm(NORM_TWO); 141 delete duf; 142 nu_sed=uf_sed_sub_iter->Norm(NORM_TWO); 143 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!"); 144 if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the layer is empty*/ 145 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 146 if((ndu_sed/nu_sed)<eps_hyd){ 147 if(VerboseConvergence()) _printf0_(" # Inner sediment convergence achieve \n"); 148 sedconverged=true; 149 } 150 delete uf_sed_sub_iter; 130 151 if(sedconverged){ 131 sedconverged=false;132 /*Checking convegence on the value of the sediment head*/133 duf=uf_sed_sub_iter->Duplicate();_assert_(duf);134 uf_sed_sub_iter->Copy(duf);135 duf->AYPX(uf_sed,-1.0);136 ndu_sed=duf->Norm(NORM_TWO);137 delete duf;138 nu_sed=uf_sed_sub_iter->Norm(NORM_TWO);139 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");140 if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the EPL is used but empty*/141 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n");142 if((ndu_sed/nu_sed)<eps_hyd){143 if(VerboseConvergence()) _printf0_(" # Inner sediment convergence achieve \n");144 sedconverged=true;145 }146 }147 delete uf_sed_sub_iter;148 149 if(sedconverged){150 if(isefficientlayer){151 /*Updating Elemental Mask*/152 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();153 analysis->ElementizeEplMask(femmodel);154 delete analysis;155 femmodel->HydrologyTransferx();156 }157 152 femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); 158 153 InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum); … … 162 157 } 163 158 } 164 165 /*Second layer*/ 159 /* }}} *//*End of the global sediment loop*/ 160 161 /* {{{ *//*Now dealing with the EPL in the same way*/ 166 162 if(isefficientlayer){ 167 163 168 164 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); 165 /*updating mask*/ 166 femmodel->HydrologyEPLupdateDomainx(); 169 167 InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); 170 168 InputUpdateFromConstantx(femmodel,false,ConvergedEnum); 171 169 femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); 170 171 eplconverged = false; 172 172 173 174 175 /*Iteration on the EPL layer*/176 eplconverged = false;177 173 for(;;){ 178 179 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)180 femmodel->HydrologyEPLupdateDomainx();181 /* Start by retrieving the EPL head slopes and compute EPL Thickness*/174 ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); 175 ug_epl->Copy(ug_epl_sub_iter); 176 177 /* {{{ *//*Start by retrieving the EPL head slopes and compute EPL Thickness*/ 182 178 if(VerboseSolution()) _printf0_("computing EPL Head slope...\n"); 183 179 femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum); … … 188 184 solutionsequence_linear(femmodel); 189 185 femmodel->HydrologyEPLThicknessx(); 186 //updating mask after the computation of the epl thickness (Allow to close too thin EPL) 187 femmodel->HydrologyEPLupdateDomainx(); 188 /* }}} */ 189 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); 190 190 191 femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); 192 193 /*Updating Elemental Mask*/ 191 /* {{{ *//*Loop on the epl layer to deal with the penalization*/ 192 for(;;){ 193 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); 194 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); 195 Reduceloadx(pf,Kfs,ys); delete Kfs; 196 delete uf_epl; 197 Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); 198 delete Kff; delete pf; delete df; 199 delete uf_epl_sub_iter; 200 uf_epl_sub_iter=uf_epl->Duplicate(); 201 uf_epl->Copy(uf_epl_sub_iter); 202 delete ug_epl; 203 Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; 204 InputUpdateFromSolutionx(femmodel,ug_epl); 205 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); 206 207 if (!eplconverged){ 208 if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n"); 209 if(num_unstable_constraints==0) eplconverged = true; 210 if (eplcount>=hydro_maxiter){ 211 _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); 212 } 213 } 214 /*Add an iteration and get out of the loop if the penalisation is converged*/ 215 eplcount++; 216 if(eplconverged) break; 217 } 218 /* }}} */ /*End of the EPL penalization loop*/ 219 220 eplconverged=false; 221 /*Update Elemental Mask and compute transfer*/ 194 222 HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); 195 223 analysis->ElementizeEplMask(femmodel); 196 224 delete analysis; 197 225 femmodel->HydrologyTransferx(); 198 199 ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter); 200 ug_epl->Copy(ug_epl_sub_iter); 201 202 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); 203 CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); 204 Reduceloadx(pf,Kfs,ys); delete Kfs; 205 delete uf_epl; 206 Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters); 207 delete Kff; delete pf; delete df; 208 delete uf_epl_sub_iter; 209 uf_epl_sub_iter=uf_epl->Duplicate(); 210 uf_epl->Copy(uf_epl_sub_iter); 211 delete ug_epl; 212 Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; 213 InputUpdateFromSolutionx(femmodel,ug_epl); 214 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); 215 216 217 if (!eplconverged){ 218 if(VerboseConvergence()) _printf0_(" # EPL unstable constraints = " << num_unstable_constraints << "\n"); 219 if(num_unstable_constraints==0) eplconverged = true; 220 if (eplcount>=hydro_maxiter){ 221 //eplconverged =true; 222 _error_(" maximum number of EPL iterations (" << hydro_maxiter << ") exceeded"); 223 } 224 } 225 eplcount++; 226 226 227 dug=ug_epl_sub_iter->Duplicate();_assert_(dug); 228 ug_epl_sub_iter->Copy(dug); 229 dug->AYPX(ug_epl,-1.0); 230 ndu_epl=dug->Norm(NORM_TWO); 231 delete dug; 232 nu_epl=ug_epl_sub_iter->Norm(NORM_TWO); 233 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!"); 234 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ 235 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 236 if((ndu_epl/nu_epl)<eps_hyd) eplconverged=true; 237 238 delete ug_epl_sub_iter; 227 239 if(eplconverged){ 228 eplconverged=false;229 230 dug=ug_epl_sub_iter->Duplicate();_assert_(dug);231 ug_epl_sub_iter->Copy(dug);232 dug->AYPX(ug_epl,-1.0);233 ndu_epl=dug->Norm(NORM_TWO);234 delete dug;235 nu_epl=ug_epl_sub_iter->Norm(NORM_TWO);236 237 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");238 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/239 if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n");240 if((ndu_epl/nu_epl)<eps_hyd)eplconverged=true;241 }242 delete ug_epl_sub_iter;243 244 if(eplconverged){245 246 240 InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum); 247 241 InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum); … … 251 245 } 252 246 } 253 254 /*System convergence check*/ 247 /* }}} */ /*End of the global EPL loop*/ 248 249 /* {{{ */ /*Now dealing with the convergence of the whole system*/ 255 250 if(!hydroconverged){ 256 251 //compute norm(du)/norm(u) … … 264 259 if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!"); 265 260 if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the Sediment is used but empty*/ 261 dug=ug_epl_main_iter->Duplicate();_assert_(dug); 262 ug_epl_main_iter->Copy(dug); 263 dug->AYPX(ug_epl,-1.0); 264 ndu_epl=dug->Norm(NORM_TWO); 265 delete dug; 266 nu_epl=ug_epl_main_iter->Norm(NORM_TWO); 267 delete ug_epl_main_iter; 268 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!"); 269 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ 266 270 if (!xIsNan<IssmDouble>(eps_hyd)){ 267 if (!isefficientlayer){ 268 if ((ndu_sed/nu_sed)<eps_hyd){ 269 if(VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n"); 270 hydroconverged=true; 271 } 272 else{ 273 if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " > " << eps_hyd*100 << " %\n"); 274 hydroconverged=false; 275 } 276 } 277 else{ 278 dug=ug_epl_main_iter->Duplicate();_assert_(dug); 279 ug_epl_main_iter->Copy(dug); 280 dug->AYPX(ug_epl,-1.0); 281 ndu_epl=dug->Norm(NORM_TWO); 282 delete dug; 283 nu_epl=ug_epl_main_iter->Norm(NORM_TWO); 284 delete ug_epl_main_iter; 285 if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!"); 286 if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/ 287 if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd*10)){ 288 if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n"); 289 hydroconverged=true; 290 } 291 else{ 292 if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 293 if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 294 hydroconverged=false; 295 } 271 if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd*10)){ 272 if (VerboseConvergence()) _printf0_(setw(50) << left << " Converged after, " << hydrocount << " iterations \n"); 273 hydroconverged=true; 274 } 275 else{ 276 if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 277 if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << " aiming lower than " << eps_hyd*100 << " %\n"); 278 hydroconverged=false; 296 279 } 297 280 } 298 281 else _printf0_(setw(50) << left << " Convergence criterion:" << ndu_sed/nu_sed*100 << " %\n"); 299 282 if (hydrocount>=hydro_maxiter){ 300 /*Hacking to get the results of non converged runs*/ 301 //hydroconverged = true; 302 _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); 283 _error_(" maximum number for hydrological global iterations (" << hydro_maxiter << ") exceeded"); 303 284 } 304 285 } 305 286 hydrocount++; 306 287 if(hydroconverged)break; 307 }308 288 } 289 /* }}} */ 309 290 InputUpdateFromSolutionx(femmodel,ug_sed); 310 291 if(isefficientlayer)InputUpdateFromSolutionx(femmodel,ug_epl); 311 312 292 /*Free ressources: */ 313 293 delete ug_epl; … … 316 296 delete uf_epl; 317 297 delete uf_epl_sub_iter; 318 298 319 299 }
Note:
See TracChangeset
for help on using the changeset viewer.
