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

Last change on this file since 27035 was 27035, checked in by Mathieu Morlighem, 3 years ago

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