source: issm/trunk/src/c/modules/StochasticForcingx/StochasticForcingx.cpp

Last change on this file was 28013, checked in by Mathieu Morlighem, 16 months ago

merged trunk-jpl and trunk for revision 28011
File size: 11.5 KB
Line 
1/*!\file StochasticForcingx
2 * \brief: compute noise terms for the StochasticForcing fields
3 */
4
5#include "./StochasticForcingx.h"
6#include "../../classes/Loads/Friction.h"
7#include "../../shared/shared.h"
8#include "../../toolkits/toolkits.h"
9#include "../../shared/Random/random.h"
10
11void StochasticForcingx(FemModel* femmodel){/*{{{*/
12
13
14 /*Retrieve parameters*/
15 bool randomflag;
16 int M,N,numfields,numtcov,my_rank;
17 int* fields = NULL;
18 int* dimensions = NULL;
19 IssmDouble* timecovariance = NULL;
20 IssmDouble* covariance = NULL;
21 femmodel->parameters->FindParam(&randomflag,StochasticForcingRandomflagEnum);
22 femmodel->parameters->FindParam(&numfields,StochasticForcingNumFieldsEnum);
23 femmodel->parameters->FindParam(&numtcov,StochasticForcingNumTimesCovarianceEnum);
24 femmodel->parameters->FindParam(&fields,&N,StochasticForcingFieldsEnum); _assert_(N==numfields);
25 femmodel->parameters->FindParam(&dimensions,&N,StochasticForcingDimensionsEnum); _assert_(N==numfields);
26 femmodel->parameters->FindParam(&timecovariance,&N,StochasticForcingTimeCovarianceEnum); _assert_(N==numtcov);
27 int dimtot=0;
28 for(int i=0;i<numfields;i++) dimtot = dimtot+dimensions[i];
29 femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==numtcov); _assert_(N==dimtot*dimtot);
30
31 /*Check if this is a timestep for new noiseterms computation*/
32 bool isstepforstoch = false;
33 IssmDouble time,dt,starttime,tstep_stoch;
34 femmodel->parameters->FindParam(&time,TimeEnum);
35 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
36 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
37 femmodel->parameters->FindParam(&tstep_stoch,StochasticForcingTimestepEnum);
38
39 /*Check if we use HydroarmaPw*/
40 bool ispwHydro;
41 femmodel->parameters->FindParam(&ispwHydro,HydrologyIsWaterPressureArmaEnum);
42
43 #ifndef _HAVE_AD_
44 if((fmod(time,tstep_stoch)<fmod((time-dt),tstep_stoch)) || (time<=starttime+dt) || tstep_stoch==dt) isstepforstoch = true;
45 #else
46 _error_("not implemented yet");
47 #endif
48
49 /*Compute noise terms*/
50 IssmDouble* timestepcovariance = xNew<IssmDouble>(dimtot*dimtot);
51 IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
52 if(isstepforstoch){
53 /*Find covariance to be applied at current time step*/
54 int itime;
55 if(numtcov>1){
56 for(int i=0;i<numtcov;i++){
57 if(time>=timecovariance[i]) itime=i;
58 }
59 }
60 else itime=0;
61 for(int i=0;i<dimtot*dimtot;i++) timestepcovariance[i] = covariance[itime*dimtot*dimtot+i];
62 my_rank=IssmComm::GetRank();
63 if(my_rank==0){
64 int fixedseed;
65 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
66 if(randomflag) fixedseed=-1;
67 else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
68 /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
69 IssmDouble* temparray = NULL;
70 multivariateNormal(&temparray,dimtot,0.0,timestepcovariance,fixedseed);
71 for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
72 xDelete<IssmDouble>(temparray);
73 }
74 ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
75 femmodel->parameters->SetParam(noiseterms,dimtot,StochasticForcingNoisetermsEnum);
76 }
77 else{
78 IssmDouble* temparray = NULL;
79 femmodel->parameters->FindParam(&temparray,&N,StochasticForcingNoisetermsEnum); _assert_(N==dimtot);
80 for(int i=0;i<dimtot;i++) noiseterms[i] = temparray[i];
81 xDelete<IssmDouble>(temparray);
82 }
83
84 int i=0;
85 for(int j=0;j<numfields;j++){
86 int dimenum_type,noiseenum_type;
87 IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
88 for(int k=0;k<dimensions[j];k++){
89 noisefield[k]=noiseterms[i+k];
90 }
91
92 int dimensionid;
93
94 /*Deal with the ARMA models*/
95 if(fields[j]==SMBarmaEnum || fields[j]==FrontalForcingsRignotarmaEnum || fields[j]==BasalforcingsDeepwaterMeltingRatearmaEnum || fields[j]==FrontalForcingsSubglacialDischargearmaEnum || (fields[j]==FrictionWaterPressureEnum && ispwHydro)){
96 switch(fields[j]){
97 case SMBarmaEnum:
98 dimenum_type = SmbBasinsIdEnum;
99 noiseenum_type = SmbARMANoiseEnum;
100 break;
101 case FrontalForcingsRignotarmaEnum:
102 dimenum_type = FrontalForcingsBasinIdEnum;
103 noiseenum_type = ThermalforcingARMANoiseEnum;
104 break;
105 case BasalforcingsDeepwaterMeltingRatearmaEnum:
106 dimenum_type = BasalforcingsLinearBasinIdEnum;
107 noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum;
108 break;
109 case FrontalForcingsSubglacialDischargearmaEnum:
110 dimenum_type = FrontalForcingsBasinIdEnum;
111 noiseenum_type = SubglacialdischargeARMANoiseEnum;
112 break;
113 case FrictionWaterPressureEnum:
114 dimenum_type = HydrologyBasinsIdEnum;
115 noiseenum_type = FrictionWaterPressureNoiseEnum;
116 break;
117 }
118 for(Object* &object:femmodel->elements->objects){
119 Element* element = xDynamicCast<Element*>(object);
120 int numvertices = element->GetNumberOfVertices();
121 IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
122 element->GetInputValue(&dimensionid,dimenum_type);
123 for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
124 element->AddInput(noiseenum_type,noise_element,P0Enum);
125 xDelete<IssmDouble>(noise_element);
126 }
127 }
128 else{
129 switch(fields[j]){
130 case SMBarmaEnum:
131 case FrontalForcingsRignotarmaEnum:
132 case BasalforcingsDeepwaterMeltingRatearmaEnum:
133 case FrontalForcingsSubglacialDischargearmaEnum:
134 /*Already done above*/
135 break;
136 case BasalforcingsSpatialDeepwaterMeltingRateEnum:
137 /*Delete BasalforcingsSpatialDeepwaterMeltingRateEnum at previous time step (required if it is transient)*/
138 femmodel->inputs->DeleteInput(BasalforcingsSpatialDeepwaterMeltingRateEnum);
139 for(Object* &object:femmodel->elements->objects){
140 Element* element = xDynamicCast<Element*>(object);
141 int numvertices = element->GetNumberOfVertices();
142 IssmDouble baselinedeepwatermelt;
143 IssmDouble deepwatermelt_tot[numvertices];
144 Input* baselinedeepwatermelt_input = NULL;
145 baselinedeepwatermelt_input = element->GetInput(BaselineBasalforcingsSpatialDeepwaterMeltingRateEnum); _assert_(baselinedeepwatermelt_input);
146 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
147 Gauss* gauss = element->NewGauss();
148 for(int i=0;i<numvertices;i++){
149 gauss->GaussVertex(i);
150 baselinedeepwatermelt_input->GetInputValue(&baselinedeepwatermelt,gauss);
151 deepwatermelt_tot[i] = baselinedeepwatermelt+noisefield[dimensionid];
152 }
153 element->AddInput(BasalforcingsSpatialDeepwaterMeltingRateEnum,&deepwatermelt_tot[0],P1DGEnum);
154 delete gauss;
155 }
156 break;
157 case DefaultCalvingEnum:
158 /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
159 femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
160 for(Object* &object:femmodel->elements->objects){
161 Element* element = xDynamicCast<Element*>(object);
162 int numvertices = element->GetNumberOfVertices();
163 IssmDouble baselinecalvingrate;
164 IssmDouble calvingrate_tot[numvertices];
165 Input* baselinecalvingrate_input = NULL;
166 baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
167 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
168 Gauss* gauss = element->NewGauss();
169 for(int i=0;i<numvertices;i++){
170 gauss->GaussVertex(i);
171 baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
172 calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
173 }
174 element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
175 delete gauss;
176 }
177 break;
178 case FloatingMeltRateEnum:
179 /*Delete BasalforcingsFloatingiceMeltingRateEnum at previous time step (required if it is transient)*/
180 femmodel->inputs->DeleteInput(BasalforcingsFloatingiceMeltingRateEnum);
181 for(Object* &object:femmodel->elements->objects){
182 Element* element = xDynamicCast<Element*>(object);
183 int numvertices = element->GetNumberOfVertices();
184 IssmDouble baselinefloatingicemeltrate;
185 IssmDouble floatingicemeltrate_tot[numvertices];
186 Input* baselinefloatingicemeltrate_input = NULL;
187 baselinefloatingicemeltrate_input = element->GetInput(BaselineBasalforcingsFloatingiceMeltingRateEnum); _assert_(baselinefloatingicemeltrate_input);
188 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
189 Gauss* gauss = element->NewGauss();
190 for(int i=0;i<numvertices;i++){
191 gauss->GaussVertex(i);
192 baselinefloatingicemeltrate_input->GetInputValue(&baselinefloatingicemeltrate,gauss);
193 /*No check for positive melt rate because basal accretion is allowed*/
194 floatingicemeltrate_tot[i] = baselinefloatingicemeltrate+noisefield[dimensionid];
195 }
196 element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&floatingicemeltrate_tot[0],P1DGEnum);
197 delete gauss;
198 }
199 break;
200 case SMBforcingEnum:
201 /*Delete SmbMassBalanceEnum at previous time step (required if it is transient)*/
202 femmodel->inputs->DeleteInput(SmbMassBalanceEnum);
203 for(Object* &object:femmodel->elements->objects){
204 Element* element = xDynamicCast<Element*>(object);
205 int numvertices = element->GetNumberOfVertices();
206 IssmDouble baselinesmb;
207 IssmDouble smb_tot[numvertices];
208 Input* baselinesmb_input = NULL;
209 baselinesmb_input = element->GetInput(BaselineSmbMassBalanceEnum); _assert_(baselinesmb_input);
210 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
211 Gauss* gauss = element->NewGauss();
212 for(int i=0;i<numvertices;i++){
213 gauss->GaussVertex(i);
214 baselinesmb_input->GetInputValue(&baselinesmb,gauss);
215 smb_tot[i] = baselinesmb+noisefield[dimensionid];
216 }
217 element->AddInput(SmbMassBalanceEnum,&smb_tot[0],P1DGEnum);
218 delete gauss;
219 }
220 break;
221 case FrictionWaterPressureEnum:
222 /*Specify that WaterPressure is stochastic*/
223 femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
224 for(Object* &object:femmodel->elements->objects){
225 Element* element = xDynamicCast<Element*>(object);
226 int numvertices = element->GetNumberOfVertices();
227 IssmDouble p_water_deterministic[numvertices];
228 IssmDouble p_water[numvertices];
229 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
230 element->SubglacialWaterPressure(FrictionWaterPressureEnum);
231 element->GetInputListOnVertices(&p_water_deterministic[0],FrictionWaterPressureEnum);
232 for(int i=0;i<numvertices;i++) p_water[i] = p_water_deterministic[i] + noisefield[dimensionid];
233 element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
234 }
235 break;
236 default:
237 _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
238 }
239 }
240 i=i+dimensions[j];
241 xDelete<IssmDouble>(noisefield);
242 }
243
244 /*Cleanup*/
245 xDelete<int>(fields);
246 xDelete<int>(dimensions);
247 xDelete<IssmDouble>(covariance);
248 xDelete<IssmDouble>(timecovariance);
249 xDelete<IssmDouble>(timestepcovariance);
250 xDelete<IssmDouble>(noiseterms);
251}/*}}}*/
Note: See TracBrowser for help on using the repository browser.