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

Last change on this file since 26744 was 26676, checked in by vverjans, 3 years ago

NEW: added FrictionWaterPressure to stochasticforcingfields, new functions Friction::IcePressure() and Friction::SubglacialWaterPressure(), clean-up Friction::EffectivePressure()
File size: 8.4 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 //VV testing (12Nov)
15 /*
16 IssmDouble timeVV,dtVV,starttimeVV;
17 femmodel->parameters->FindParam(&timeVV,TimeEnum);
18 femmodel->parameters->FindParam(&dtVV,TimesteppingTimeStepEnum);
19 femmodel->parameters->FindParam(&starttimeVV,TimesteppingStartTimeEnum);
20 IssmDouble valMean = 0;
21 IssmDouble valSdev = 0.5;
22 int seed;
23 //seed = reCast<int,IssmDouble>((timeVV-starttimeVV)/dtVV);
24 seed = -1;
25 IssmDouble rdmVV;
26 univariateNormal_test0(&rdmVV,valMean,valSdev,seed);
27 _printf_("VV rdmVV: "<<rdmVV<<'\n');
28 */
29
30 /*Retrieve parameters*/
31 bool randomflag;
32 int M,N,numfields,my_rank;
33 int* fields = NULL;
34 int* dimensions = NULL;
35 IssmDouble* covariance = NULL;
36 femmodel->parameters->FindParam(&randomflag,StochasticForcingRandomflagEnum);
37 femmodel->parameters->FindParam(&numfields,StochasticForcingNumFieldsEnum);
38 femmodel->parameters->FindParam(&fields,&N,StochasticForcingFieldsEnum); _assert_(N==numfields);
39 femmodel->parameters->FindParam(&dimensions,&N,StochasticForcingDimensionsEnum); _assert_(N==numfields);
40 int dimtot=0;
41 for(int i=0;i<numfields;i++) dimtot = dimtot+dimensions[i];
42 femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==dimtot); _assert_(N==dimtot);
43
44 /*Compute noise terms*/
45 IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
46 my_rank=IssmComm::GetRank();
47 if(my_rank==0){
48 int fixedseed;
49 IssmDouble time,dt,starttime;
50 femmodel->parameters->FindParam(&time,TimeEnum);
51 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
52 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
53 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
54 if(randomflag) fixedseed=-1;
55 else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
56 /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
57 IssmDouble* temparray = NULL;
58 multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
59 for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
60 xDelete<IssmDouble>(temparray);
61 }
62 ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
63
64 int i=0;
65 for(int j=0;j<numfields;j++){
66 int dimenum_type,noiseenum_type;
67 IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
68 for(int k=0;k<dimensions[j];k++){
69 noisefield[k]=noiseterms[i+k];
70 }
71
72 int dimensionid;
73
74 /*Deal with the autoregressive models*/
75 if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
76 switch(fields[j]){
77 case SMBautoregressionEnum:
78 dimenum_type = SmbBasinsIdEnum;
79 noiseenum_type = SmbAutoregressionNoiseEnum;
80 break;
81 case FrontalForcingsRignotAutoregressionEnum:
82 dimenum_type = FrontalForcingsBasinIdEnum;
83 noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
84 break;
85 }
86 for(Object* &object:femmodel->elements->objects){
87 Element* element = xDynamicCast<Element*>(object);
88 int numvertices = element->GetNumberOfVertices();
89 IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
90 element->GetInputValue(&dimensionid,dimenum_type);
91 for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
92 element->AddInput(noiseenum_type,noise_element,P0Enum);
93 xDelete<IssmDouble>(noise_element);
94 }
95 }
96 else{
97 switch(fields[j]){
98 case SMBautoregressionEnum:
99 case FrontalForcingsRignotAutoregressionEnum:
100 /*Already done above*/
101 break;
102 case DefaultCalvingEnum:
103 /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
104 femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
105 for(Object* &object:femmodel->elements->objects){
106 Element* element = xDynamicCast<Element*>(object);
107 int numvertices = element->GetNumberOfVertices();
108 IssmDouble baselinecalvingrate;
109 IssmDouble calvingrate_tot[numvertices];
110 Input* baselinecalvingrate_input = NULL;
111 baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
112 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
113 Gauss* gauss = element->NewGauss();
114 for(int i=0;i<numvertices;i++){
115 gauss->GaussVertex(i);
116 baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
117 calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
118 }
119 element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
120 delete gauss;
121 }
122 break;
123 case FloatingMeltRateEnum:
124 /*Delete BasalforcingsFloatingiceMeltingRateEnum at previous time step (required if it is transient)*/
125 femmodel->inputs->DeleteInput(BasalforcingsFloatingiceMeltingRateEnum);
126 for(Object* &object:femmodel->elements->objects){
127 Element* element = xDynamicCast<Element*>(object);
128 int numvertices = element->GetNumberOfVertices();
129 IssmDouble baselinefloatingicemeltrate;
130 IssmDouble floatingicemeltrate_tot[numvertices];
131 Input* baselinefloatingicemeltrate_input = NULL;
132 baselinefloatingicemeltrate_input = element->GetInput(BaselineBasalforcingsFloatingiceMeltingRateEnum); _assert_(baselinefloatingicemeltrate_input);
133 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
134 Gauss* gauss = element->NewGauss();
135 for(int i=0;i<numvertices;i++){
136 gauss->GaussVertex(i);
137 baselinefloatingicemeltrate_input->GetInputValue(&baselinefloatingicemeltrate,gauss);
138 floatingicemeltrate_tot[i] = max(0.0,baselinefloatingicemeltrate+noisefield[dimensionid]);
139 }
140 element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&floatingicemeltrate_tot[0],P1DGEnum);
141 delete gauss;
142 }
143 break;
144 //VV working(29Nov2021)
145 case FrictionWaterPressureEnum:
146 /*Specify that WaterPressure is stochastic*/
147 femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
148 for(Object* &object:femmodel->elements->objects){
149 Element* element = xDynamicCast<Element*>(object);
150 int numvertices = element->GetNumberOfVertices();
151 IssmDouble p_water_deterministic[numvertices];
152 IssmDouble p_water[numvertices];
153 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
154 Gauss* gauss=element->NewGauss();
155 Friction* friction = new Friction(element);
156 for(int i=0;i<numvertices;i++){
157 gauss->GaussVertex(i);
158 p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
159 p_water[i] = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
160 p_water[i] = max(0.0,p_water[i]);
161 }
162 element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
163 delete gauss;
164 delete friction;
165 }
166 break;
167
168 /*
169 case FrictionEffectivePressureEnum:
170 femmodel->parameters->SetParam(true,StochasticForcingIsEffectivePressureEnum);
171 for(Object* &object:femmodel->elements->objects){
172 Element* element = xDynamicCast<Element*>(object);
173 int numvertices = element->GetNumberOfVertices();
174 IssmDouble Neff[numvertices];
175 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
176 Gauss* gauss=element->NewGauss();
177 Friction* friction = new Friction(element);
178 for(int i=0;i<numvertices;i++){
179 gauss->GaussVertex(i);
180 Neff[i] = friction->EffectivePressure(gauss);
181 Neff[i] = Neff[i]+noisefield[dimensionid];
182 }
183 element->AddInput(FrictionEffectivePressureEnum,Neff,P1DGEnum);
184 delete gauss;
185 delete friction;
186 }
187 break;
188 */
189
190 default:
191 _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
192 }
193 }
194 i=i+dimensions[j];
195 xDelete<IssmDouble>(noisefield);
196 }
197
198 /*Cleanup*/
199 xDelete<int>(fields);
200 xDelete<int>(dimensions);
201 xDelete<IssmDouble>(covariance);
202 xDelete<IssmDouble>(noiseterms);
203}/*}}}*/
Note: See TracBrowser for help on using the repository browser.