![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23665
- Timestamp:
- 01/28/19 02:29:02 (6 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 3 edited
-
cores/hydrology_core.cpp (modified) (9 diffs)
-
cores/transient_core.cpp (modified) (2 diffs)
-
modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/cores/hydrology_core.cpp
r23661 r23665 19 19 int solution_type; 20 20 int numoutputs = 0; 21 int smboutputs; 21 22 bool save_results; 22 23 bool modify_loads = true; 24 bool issmb; 23 25 char **requested_outputs = NULL; 26 char **requested_smb_outputs = NULL; 24 27 IssmDouble ThawedNodes; 25 28 … … 30 33 femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum); 31 34 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum); 35 femmodel->parameters->FindParam(&issmb,TransientIssmbEnum); 36 32 37 /*Using the Shreve based Model*/ 33 38 if (hydrology_model==HydrologyshreveEnum){ … … 40 45 femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); 41 46 solutionsequence_nonlinear(femmodel,modify_loads); 42 43 47 /*transfer water column thickness to old water column thickness: */ 44 48 InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); … … 50 54 /*intermediary: */ 51 55 bool isefficientlayer; 52 int hydrostep,hydroslices ;53 IssmDouble time, init_time,hydrotime,yts;56 int hydrostep,hydroslices,numaveragedinput; 57 IssmDouble time,hydrotime,yts; 54 58 IssmDouble dt,hydrodt; 55 59 /*SMB related */ 60 int smb_model; 61 62 /*recover parameters: */ 56 63 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); 57 64 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); … … 60 67 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); 61 68 69 /*recover SMB related parameters: */ 70 if(issmb){ 71 femmodel->parameters->FindParam(&smb_model,SmbEnum); 72 femmodel->parameters->FindParam(&smboutputs,SmbNumRequestedOutputsEnum); 73 if(smboutputs) femmodel->parameters->FindParam(&requested_smb_outputs,&smboutputs,SmbRequestedOutputsEnum); 74 } 75 62 76 /*first we exclude frozen nodes of the solved nodes*/ 63 77 femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); 64 78 femmodel->HydrologyIDSupdateDomainx(&ThawedNodes); 79 65 80 if(ThawedNodes>0){ 66 init_time=time-dt; //getting the time back to the start of the timestep 67 hydrotime=init_time; 81 hydrotime=time-dt; //getting the time back to the start of the timestep 68 82 hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider 69 83 hydrostep=0; 70 84 femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt)); 85 71 86 if(hydroslices>1){ 72 87 /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/ 88 numaveragedinput = 2; 89 std::vector<int> inputtostack = {EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; 90 std::vector<int> stackedinput = {EffectivePressureStackedEnum,SedimentHeadStackedEnum}; 91 std::vector<int> averagedinput = {EffectivePressureEnum,SedimentHeadEnum}; 92 73 93 if (isefficientlayer){ 74 int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; 75 int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; 76 int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; 77 femmodel->InitMeanOutputx(&stackedinput[0],4); 78 while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts 79 hydrostep+=1; 80 hydrotime+=hydrodt; 81 if(VerboseSolution()) _printf0_("hydro iteration " << hydrostep << "/" << hydroslices << " time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n"); 82 83 /*save preceding timestep*/ 84 InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); 94 numaveragedinput = 4; 95 inputtostack = {EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; 96 stackedinput = {EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; 97 averagedinput = {EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; 98 } 99 100 femmodel->InitMeanOutputx(&stackedinput[0],numaveragedinput); 101 102 while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts 103 hydrostep+=1; 104 hydrotime+=hydrodt; 105 /*Setting substep time as global time*/ 106 femmodel->parameters->SetParam(hydrotime,TimeEnum); 107 if(VerboseSolution()) _printf0_("sub iteration " << hydrostep << "/" << hydroslices << " time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n"); 108 if(issmb){ 109 if(VerboseSolution()) _printf0_(" computing mass balance\n"); 110 SmbAnalysis* analysis = new SmbAnalysis(); 111 analysis->Core(femmodel); 112 delete analysis; 113 } 114 if(VerboseSolution()) _printf0_(" computing water heads\n"); 115 /*save preceding timestep*/ 116 InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); 117 if(isefficientlayer){ 85 118 InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum); 86 119 InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum); 87 /*Proceed now to heads computations*/88 solutionsequence_hydro_nonlinear(femmodel);89 /*If we have a sub-timestep we stack the variables here*/90 femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4);91 120 } 92 femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4); 93 } 94 else{ 95 int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; 96 int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum}; 97 int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum}; 98 femmodel->InitMeanOutputx(&stackedinput[0],2); 99 while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts 100 hydrotime+=hydrodt; 101 /*save preceding timestep*/ 102 InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); 103 /*Proceed now to heads computations*/ 104 solutionsequence_hydro_nonlinear(femmodel); 105 /*If we have a sub-timestep we stack the variables here*/ 106 femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2); 107 } 108 femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2); 109 } 121 /*Proceed now to heads computations*/ 122 solutionsequence_hydro_nonlinear(femmodel); 123 /*If we have a sub-timestep we stack the variables here*/ 124 femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],numaveragedinput); 125 } 126 /*Reseting to global time*/ 127 femmodel->parameters->SetParam(time,TimeEnum); 128 femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],numaveragedinput); 110 129 } 111 130 else{ 131 if(issmb){ 132 if(VerboseSolution()) _printf0_(" computing mass balance\n"); 133 SmbAnalysis* analysis = new SmbAnalysis(); 134 analysis->Core(femmodel); 135 delete analysis; 136 } 112 137 InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); 113 138 if (isefficientlayer){ … … 116 141 } 117 142 /*Proceed now to heads computations*/ 143 if(VerboseSolution()) _printf0_(" computing water heads\n"); 118 144 solutionsequence_hydro_nonlinear(femmodel); 119 145 } 120 146 } 121 } 122 else if (hydrology_model==HydrologyshaktiEnum){ 147 else{ 148 /* If everything is frozen we still need smb */ 149 if(issmb){ 150 if(VerboseSolution()) _printf0_(" computing mass balance\n"); 151 SmbAnalysis* analysis = new SmbAnalysis(); 152 analysis->Core(femmodel); 153 delete analysis; 154 } 155 } 156 } 157 158 /*Using the SHAKTI model*/ 159 else if (hydrology_model==HydrologyshaktiEnum){ 123 160 femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum); 124 161 InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); … … 129 166 delete analysis; 130 167 } 168 169 /*Using the PISM hydrology model*/ 131 170 else if (hydrology_model==HydrologypismEnum){ 132 171 femmodel->SetCurrentConfiguration(HydrologyPismAnalysisEnum); … … 143 182 if(hydrology_model==HydrologydcEnum && ThawedNodes==0){ 144 183 if(VerboseSolution()) _printf0_(" No thawed node hydro is skiped \n");} 184 else if (hydrology_model==HydrologydcEnum && issmb){ 185 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); 186 femmodel->RequestedOutputsx(&femmodel->results,requested_smb_outputs,smboutputs); 187 } 145 188 else{ 146 189 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); … … 154 197 xDelete<char*>(requested_outputs); 155 198 } 156 199 if(issmb){ 200 if(smboutputs){ 201 for(int i=0;i<smboutputs;i++){ 202 xDelete<char>(requested_smb_outputs[i]); 203 } 204 xDelete<char*>(requested_smb_outputs); 205 } 206 } 157 207 /*End profiler*/ 158 208 femmodel->profiler->Stop(HYDROLOGYCORE); -
issm/trunk-jpl/src/c/cores/transient_core.cpp
r23649 r23665 74 74 #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_) 75 75 if(amr_frequency){ 76 76 femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); 77 77 if(amr_restart) femmodel->ReMesh(); 78 78 } … … 379 379 thermal_core(femmodel); 380 380 } 381 382 /*shifting smb position to have runoff value*/ 383 if(issmb) smb_core(femmodel); 384 385 if(ishydrology) hydrology_core(femmodel); 381 /* Using Hydrology dc coupled we need to compute smb in the hydrology inner time loop*/ 382 if(ishydrology){ 383 int hydrology_model; 384 hydrology_core(femmodel); 385 femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum); 386 if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel); 387 } 388 else{ 389 if(issmb) smb_core(femmodel); 390 } 386 391 387 392 if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) stressbalance_core(femmodel); -
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
r23644 r23665 427 427 void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/ 428 428 429 // void SurfaceMassBalancex(hd,agd,ni){430 // INPUT parameters: ni: working size of arrays431 // INPUT: surface elevation (m): hd(NA)432 // OUTPUT: mass-balance (m/yr ice): agd(NA)433 434 429 for(int i=0;i<femmodel->elements->Size();i++){ 435 430 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
Note:
See TracChangeset
for help on using the changeset viewer.
