[19984] | 1 | #include "./SealevelriseAnalysis.h"
|
---|
[24959] | 2 | #include <math.h>
|
---|
[19984] | 3 | #include "../toolkits/toolkits.h"
|
---|
| 4 | #include "../classes/classes.h"
|
---|
[25379] | 5 | #include "../classes/Inputs/TransientInput.h"
|
---|
[19984] | 6 | #include "../shared/shared.h"
|
---|
| 7 | #include "../modules/modules.h"
|
---|
| 8 |
|
---|
| 9 | /*Model processing*/
|
---|
| 10 | void SealevelriseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
|
---|
| 11 | /*No constraints*/
|
---|
| 12 | }/*}}}*/
|
---|
| 13 | void SealevelriseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
|
---|
| 14 | /*No loads*/
|
---|
| 15 | }/*}}}*/
|
---|
[23585] | 16 | void SealevelriseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
|
---|
[19984] | 17 | ::CreateNodes(nodes,iomodel,SealevelriseAnalysisEnum,P1Enum);
|
---|
| 18 | }/*}}}*/
|
---|
| 19 | int SealevelriseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
|
---|
| 20 | return 1;
|
---|
| 21 | }/*}}}*/
|
---|
[25379] | 22 | void SealevelriseAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
|
---|
[19984] | 23 |
|
---|
[22955] | 24 | int geodetic=0;
|
---|
[24469] | 25 | int dslmodel=0;
|
---|
[22955] | 26 |
|
---|
[19984] | 27 | /*Update elements: */
|
---|
| 28 | int counter=0;
|
---|
| 29 | for(int i=0;i<iomodel->numberofelements;i++){
|
---|
| 30 | if(iomodel->my_elements[i]){
|
---|
| 31 | Element* element=(Element*)elements->GetObjectByOffset(counter);
|
---|
[25379] | 32 | element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
|
---|
[19984] | 33 | counter++;
|
---|
| 34 | }
|
---|
| 35 | }
|
---|
| 36 |
|
---|
| 37 | /*Create inputs: */
|
---|
[25379] | 38 | iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
|
---|
| 39 | iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
|
---|
[25118] | 40 | iomodel->FetchData(&geodetic,"md.solidearth.settings.computesealevelchange");
|
---|
[25379] | 41 | iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.icethicknesschange",SurfaceloadIceThicknessChangeEnum);
|
---|
| 42 | iomodel->FetchDataToInput(inputs,elements,"md.solidearth.sealevel",SealevelEnum,0);
|
---|
| 43 | iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
|
---|
| 44 | iomodel->FetchDataToInput(inputs,elements,"md.solidearth.surfaceload.waterheightchange",SurfaceloadWaterHeightChangeEnum);
|
---|
[24469] | 45 |
|
---|
| 46 | /*dynamic sea level: */
|
---|
| 47 | iomodel->FetchData(&dslmodel,"md.dsl.model");
|
---|
[24505] | 48 | if (dslmodel==1){ /*standard dsl model:{{{*/
|
---|
[24481] | 49 |
|
---|
| 50 | /*deal with global mean steric rate: */
|
---|
| 51 | IssmDouble* str=NULL;
|
---|
| 52 | IssmDouble* times = NULL;
|
---|
| 53 | int M,N;
|
---|
| 54 |
|
---|
| 55 | /*fetch str vector:*/
|
---|
| 56 | iomodel->FetchData(&str,&M,&N,"md.dsl.global_average_thermosteric_sea_level_change"); _assert_(M==2);
|
---|
| 57 |
|
---|
| 58 | //recover time vector:
|
---|
| 59 | times=xNew<IssmDouble>(N);
|
---|
| 60 | for(int t=0;t<N;t++) times[t] = str[N+t];
|
---|
| 61 |
|
---|
| 62 | /*create transient input: */
|
---|
[25379] | 63 | inputs->SetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,times,N);
|
---|
| 64 | TransientInput* transientinput = inputs->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
|
---|
[24481] | 65 |
|
---|
| 66 |
|
---|
[25539] | 67 | for(Object* & object : elements->objects){
|
---|
| 68 | Element* element=xDynamicCast<Element*>(object);
|
---|
[24481] | 69 |
|
---|
| 70 | for(int t=0;t<N;t++){
|
---|
| 71 | switch(element->ObjectEnum()){
|
---|
| 72 | case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
|
---|
| 73 | case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
|
---|
| 74 | default: _error_("Not implemented yet");
|
---|
| 75 | }
|
---|
| 76 | }
|
---|
| 77 | }
|
---|
| 78 |
|
---|
| 79 | /*cleanup:*/
|
---|
| 80 | xDelete<IssmDouble>(times);
|
---|
| 81 | iomodel->DeleteData(str,"md.dsl.global_average_thermosteric_sea_level_change");
|
---|
| 82 |
|
---|
| 83 | /*deal with dynamic sea level fields: */
|
---|
[25379] | 84 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
|
---|
| 85 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloorEnum);
|
---|
[24481] | 86 |
|
---|
[24505] | 87 | } /*}}}*/
|
---|
| 88 | else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/
|
---|
| 89 |
|
---|
| 90 | /*variables:*/
|
---|
| 91 | int nummodels;
|
---|
| 92 | IssmDouble** pstr=NULL;
|
---|
| 93 | IssmDouble* str=NULL;
|
---|
| 94 | IssmDouble* times = NULL;
|
---|
| 95 | int* pM = NULL;
|
---|
| 96 | int* pN = NULL;
|
---|
| 97 | int M,N;
|
---|
| 98 |
|
---|
[25275] | 99 | /*deal with dsl.global_average_thermosteric_sea_level_change {{{*/
|
---|
[24505] | 100 | iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change");
|
---|
| 101 |
|
---|
| 102 | /*go through the mat array and create a dataset of transient inputs:*/
|
---|
| 103 | for (int i=0;i<nummodels;i++){
|
---|
| 104 |
|
---|
| 105 | M=pM[i];
|
---|
| 106 | N=pN[i];
|
---|
| 107 | str=pstr[i];
|
---|
| 108 |
|
---|
| 109 | //recover time vector:
|
---|
| 110 | times=xNew<IssmDouble>(N);
|
---|
| 111 | for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
|
---|
| 112 |
|
---|
[25379] | 113 | TransientInput* transientinput=inputs->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
|
---|
[24505] | 114 |
|
---|
[25539] | 115 | for(Object* & object : elements->objects){
|
---|
| 116 | Element* element=xDynamicCast<Element*>(object);
|
---|
[24505] | 117 |
|
---|
| 118 | for(int t=0;t<N;t++){
|
---|
| 119 | switch(element->ObjectEnum()){
|
---|
| 120 | case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
|
---|
| 121 | case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
|
---|
| 122 | default: _error_("Not implemented yet");
|
---|
| 123 | }
|
---|
| 124 | }
|
---|
| 125 | }
|
---|
| 126 | xDelete<IssmDouble>(times);
|
---|
| 127 | }
|
---|
| 128 | /*Delete data:*/
|
---|
| 129 | for(int i=0;i<nummodels;i++){
|
---|
| 130 | IssmDouble* str=pstr[i];
|
---|
| 131 | xDelete<IssmDouble>(str);
|
---|
| 132 | }
|
---|
| 133 | xDelete<IssmDouble*>(pstr);
|
---|
| 134 | xDelete<int>(pM);
|
---|
| 135 | xDelete<int>(pN);
|
---|
| 136 | /*}}}*/
|
---|
[25379] | 137 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_surface_height_change_above_geoid",DslSeaSurfaceHeightChangeAboveGeoidEnum);
|
---|
| 138 | iomodel->FetchDataToInput(inputs,elements,"md.dsl.sea_water_pressure_change_at_sea_floor",DslSeaWaterPressureChangeAtSeaFloorEnum);
|
---|
[24505] | 139 |
|
---|
| 140 | } /*}}}*/
|
---|
[24469] | 141 | else _error_("Dsl model " << dslmodel << " not implemented yet!");
|
---|
[25060] | 142 |
|
---|
[25066] | 143 |
|
---|
[25060] | 144 |
|
---|
[22968] | 145 | /*Initialize cumdeltalthickness and sealevel rise rate input*/
|
---|
[25379] | 146 | InputUpdateFromConstantx(inputs,elements,0.,SealevelriseCumDeltathicknessEnum);
|
---|
| 147 | InputUpdateFromConstantx(inputs,elements,0.,SealevelNEsaRateEnum);
|
---|
| 148 | InputUpdateFromConstantx(inputs,elements,0.,SealevelUEsaRateEnum);
|
---|
| 149 | InputUpdateFromConstantx(inputs,elements,0.,SealevelRSLRateEnum);
|
---|
| 150 | InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticMaskEnum);
|
---|
| 151 | InputUpdateFromConstantx(inputs,elements,0.,SealevelEustaticOceanMaskEnum);
|
---|
[22968] | 152 |
|
---|
[19984] | 153 | }/*}}}*/
|
---|
| 154 | void SealevelriseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
|
---|
| 155 |
|
---|
[19986] | 156 | int nl;
|
---|
[19984] | 157 | IssmDouble* love_h=NULL;
|
---|
| 158 | IssmDouble* love_k=NULL;
|
---|
[20999] | 159 | IssmDouble* love_l=NULL;
|
---|
[25139] | 160 | IssmDouble* love_th=NULL;
|
---|
| 161 | IssmDouble* love_tk=NULL;
|
---|
| 162 | IssmDouble* love_tl=NULL;
|
---|
[24469] | 163 | int dslmodel=0;
|
---|
[23066] | 164 |
|
---|
[24921] | 165 | IssmDouble* G_rigid = NULL;
|
---|
| 166 | IssmDouble* G_rigid_local = NULL;
|
---|
[20028] | 167 | IssmDouble* G_elastic = NULL;
|
---|
[20033] | 168 | IssmDouble* G_elastic_local = NULL;
|
---|
[20999] | 169 | IssmDouble* U_elastic = NULL;
|
---|
| 170 | IssmDouble* U_elastic_local = NULL;
|
---|
| 171 | IssmDouble* H_elastic = NULL;
|
---|
| 172 | IssmDouble* H_elastic_local = NULL;
|
---|
[20033] | 173 | int M,m,lower_row,upper_row;
|
---|
| 174 | IssmDouble degacc=.01;
|
---|
[24959] | 175 | IssmDouble planetradius=0;
|
---|
| 176 | IssmDouble planetarea=0;
|
---|
[19984] | 177 |
|
---|
[20036] | 178 | int numoutputs;
|
---|
| 179 | char** requestedoutputs = NULL;
|
---|
| 180 |
|
---|
[20136] | 181 | /*transition vectors: */
|
---|
| 182 | IssmDouble **transitions = NULL;
|
---|
| 183 | int *transitions_M = NULL;
|
---|
| 184 | int *transitions_N = NULL;
|
---|
| 185 | int ntransitions;
|
---|
| 186 |
|
---|
[19986] | 187 | /*some constant parameters: */
|
---|
[24531] | 188 | parameters->AddObject(iomodel->CopyConstantObject("md.dsl.model",DslModelEnum));
|
---|
[25118] | 189 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
|
---|
| 190 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
|
---|
| 191 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
|
---|
| 192 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum));
|
---|
| 193 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
|
---|
| 194 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
|
---|
| 195 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
|
---|
| 196 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
|
---|
| 197 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
|
---|
| 198 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
|
---|
| 199 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
|
---|
| 200 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
|
---|
| 201 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
|
---|
| 202 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
|
---|
[19986] | 203 |
|
---|
[24959] | 204 | /*compute planet area and plug into parameters:*/
|
---|
[25118] | 205 | iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
|
---|
[24959] | 206 | planetarea=4*PI*planetradius*planetradius;
|
---|
[25118] | 207 | parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
|
---|
[24959] | 208 |
|
---|
[24505] | 209 | /*Deal with dsl multi-model ensembles: {{{*/
|
---|
| 210 | iomodel->FetchData(&dslmodel,"md.dsl.model");
|
---|
| 211 | parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
|
---|
| 212 | if(dslmodel==2){
|
---|
[25273] | 213 | IssmDouble modelid;
|
---|
[24505] | 214 | int nummodels;
|
---|
[25273] | 215 |
|
---|
| 216 | /*create double param, not int param, because Dakota will be updating it as
|
---|
| 217 | * a double potentially: */
|
---|
| 218 | iomodel->FetchData(&modelid,"md.dsl.modelid");
|
---|
[25274] | 219 | parameters->AddObject(new DoubleParam(DslModelidEnum,modelid));
|
---|
[24505] | 220 | parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
|
---|
| 221 | iomodel->FetchData(&nummodels,"md.dsl.nummodels");
|
---|
| 222 |
|
---|
| 223 | /*quick checks: */
|
---|
| 224 | if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!");
|
---|
| 225 | if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!");
|
---|
| 226 | } /*}}}*/
|
---|
[24481] | 227 | /*Deal with elasticity {{{*/
|
---|
[24921] | 228 | /*love numbers: */
|
---|
[25146] | 229 | iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
|
---|
| 230 | iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
|
---|
| 231 | iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
|
---|
| 232 | iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
|
---|
| 233 | iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
|
---|
| 234 | iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
|
---|
| 235 | parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
|
---|
[19984] | 236 |
|
---|
[25139] | 237 | parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
|
---|
| 238 | parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
|
---|
| 239 | parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
|
---|
| 240 | parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
|
---|
| 241 | parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
|
---|
| 242 | parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
|
---|
| 243 |
|
---|
[24921] | 244 | /*compute elastic green function for a range of angles*/
|
---|
[25118] | 245 | iomodel->FetchData(°acc,"md.solidearth.settings.degacc");
|
---|
[24921] | 246 | M=reCast<int,IssmDouble>(180./degacc+1.);
|
---|
[20022] | 247 |
|
---|
[24921] | 248 | // AD performance is sensitive to calls to ensurecontiguous.
|
---|
| 249 | // // Providing "t" will cause ensurecontiguous to be called.
|
---|
| 250 | #ifdef _HAVE_AD_
|
---|
| 251 | G_rigid=xNew<IssmDouble>(M,"t");
|
---|
| 252 | G_elastic=xNew<IssmDouble>(M,"t");
|
---|
| 253 | U_elastic=xNew<IssmDouble>(M,"t");
|
---|
| 254 | H_elastic=xNew<IssmDouble>(M,"t");
|
---|
| 255 | #else
|
---|
| 256 | G_rigid=xNew<IssmDouble>(M);
|
---|
| 257 | G_elastic=xNew<IssmDouble>(M);
|
---|
| 258 | U_elastic=xNew<IssmDouble>(M);
|
---|
| 259 | H_elastic=xNew<IssmDouble>(M);
|
---|
| 260 | #endif
|
---|
[21742] | 261 |
|
---|
[24921] | 262 | /*compute combined legendre + love number (elastic green function:*/
|
---|
| 263 | m=DetermineLocalSize(M,IssmComm::GetComm());
|
---|
| 264 | GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
|
---|
| 265 | // AD performance is sensitive to calls to ensurecontiguous.
|
---|
| 266 | // // Providing "t" will cause ensurecontiguous to be called.
|
---|
| 267 | #ifdef _HAVE_AD_
|
---|
| 268 | G_elastic_local=xNew<IssmDouble>(m,"t");
|
---|
| 269 | G_rigid_local=xNew<IssmDouble>(m,"t");
|
---|
| 270 | U_elastic_local=xNew<IssmDouble>(m,"t");
|
---|
| 271 | H_elastic_local=xNew<IssmDouble>(m,"t");
|
---|
| 272 | #else
|
---|
| 273 | G_elastic_local=xNew<IssmDouble>(m);
|
---|
| 274 | G_rigid_local=xNew<IssmDouble>(m);
|
---|
| 275 | U_elastic_local=xNew<IssmDouble>(m);
|
---|
| 276 | H_elastic_local=xNew<IssmDouble>(m);
|
---|
| 277 | #endif
|
---|
[21742] | 278 |
|
---|
[24921] | 279 | for(int i=lower_row;i<upper_row;i++){
|
---|
| 280 | IssmDouble alpha,x;
|
---|
| 281 | alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
|
---|
[20033] | 282 |
|
---|
[24921] | 283 | G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
|
---|
| 284 | G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
|
---|
| 285 | U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
|
---|
| 286 | H_elastic_local[i-lower_row]= 0;
|
---|
| 287 | IssmDouble Pn = 0.;
|
---|
| 288 | IssmDouble Pn1 = 0.;
|
---|
| 289 | IssmDouble Pn2 = 0.;
|
---|
| 290 | IssmDouble Pn_p = 0.;
|
---|
| 291 | IssmDouble Pn_p1 = 0.;
|
---|
| 292 | IssmDouble Pn_p2 = 0.;
|
---|
[20029] | 293 |
|
---|
[24921] | 294 | for (int n=0;n<nl;n++) {
|
---|
| 295 | IssmDouble deltalove_G;
|
---|
| 296 | IssmDouble deltalove_U;
|
---|
[21908] | 297 |
|
---|
[24921] | 298 | deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
|
---|
| 299 | deltalove_U = (love_h[n]-love_h[nl-1]);
|
---|
[20033] | 300 |
|
---|
[24921] | 301 | /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
|
---|
| 302 | if(n==0){
|
---|
| 303 | Pn=1;
|
---|
| 304 | Pn_p=0;
|
---|
| 305 | }
|
---|
| 306 | else if(n==1){
|
---|
| 307 | Pn = cos(alpha);
|
---|
| 308 | Pn_p = 1;
|
---|
| 309 | }
|
---|
| 310 | else{
|
---|
| 311 | Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
|
---|
| 312 | Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
|
---|
| 313 | }
|
---|
| 314 | Pn2=Pn1; Pn1=Pn;
|
---|
| 315 | Pn_p2=Pn_p1; Pn_p1=Pn_p;
|
---|
[23066] | 316 |
|
---|
[24921] | 317 | G_elastic_local[i-lower_row] += deltalove_G*Pn; // gravitational potential
|
---|
| 318 | U_elastic_local[i-lower_row] += deltalove_U*Pn; // vertical (up) displacement
|
---|
| 319 | H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p; // horizontal displacements
|
---|
[20033] | 320 | }
|
---|
[24921] | 321 | }
|
---|
[20033] | 322 |
|
---|
[24921] | 323 | /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
|
---|
| 324 | int* recvcounts=xNew<int>(IssmComm::GetSize());
|
---|
| 325 | int* displs=xNew<int>(IssmComm::GetSize());
|
---|
[20033] | 326 |
|
---|
[24921] | 327 | //recvcounts:
|
---|
| 328 | ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
|
---|
[20033] | 329 |
|
---|
[24921] | 330 | /*displs: */
|
---|
| 331 | ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
|
---|
[20033] | 332 |
|
---|
[24921] | 333 | /*All gather:*/
|
---|
| 334 | ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
---|
| 335 | ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
---|
| 336 | ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
---|
| 337 | ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
|
---|
| 338 | /*free ressources: */
|
---|
| 339 | xDelete<int>(recvcounts);
|
---|
| 340 | xDelete<int>(displs);
|
---|
[20033] | 341 |
|
---|
[24921] | 342 | /*}}}*/
|
---|
[20033] | 343 |
|
---|
[24921] | 344 | /*Avoid singularity at 0: */
|
---|
| 345 | G_rigid[0]=G_rigid[1];
|
---|
| 346 | parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
|
---|
| 347 | G_elastic[0]=G_elastic[1];
|
---|
| 348 | parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
|
---|
| 349 | U_elastic[0]=U_elastic[1];
|
---|
| 350 | parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
|
---|
| 351 | H_elastic[0]=H_elastic[1];
|
---|
| 352 | parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
|
---|
[20028] | 353 |
|
---|
[24921] | 354 | /*free ressources: */
|
---|
| 355 | xDelete<IssmDouble>(love_h);
|
---|
| 356 | xDelete<IssmDouble>(love_k);
|
---|
| 357 | xDelete<IssmDouble>(love_l);
|
---|
[25139] | 358 | xDelete<IssmDouble>(love_th);
|
---|
| 359 | xDelete<IssmDouble>(love_tk);
|
---|
| 360 | xDelete<IssmDouble>(love_tl);
|
---|
[24921] | 361 | xDelete<IssmDouble>(G_rigid);
|
---|
| 362 | xDelete<IssmDouble>(G_rigid_local);
|
---|
| 363 | xDelete<IssmDouble>(G_elastic);
|
---|
| 364 | xDelete<IssmDouble>(G_elastic_local);
|
---|
| 365 | xDelete<IssmDouble>(U_elastic);
|
---|
| 366 | xDelete<IssmDouble>(U_elastic_local);
|
---|
| 367 | xDelete<IssmDouble>(H_elastic);
|
---|
| 368 | xDelete<IssmDouble>(H_elastic_local);
|
---|
[24977] | 369 |
|
---|
| 370 | /*Indicate we have not yet ran the Geometry Core module: */
|
---|
| 371 | parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
|
---|
| 372 |
|
---|
[24921] | 373 | /*}}}*/
|
---|
[24481] | 374 | /*Transitions:{{{ */
|
---|
[25118] | 375 | iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");
|
---|
[20136] | 376 | if(transitions){
|
---|
| 377 | parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
|
---|
| 378 |
|
---|
| 379 | for(int i=0;i<ntransitions;i++){
|
---|
| 380 | IssmDouble* transition=transitions[i];
|
---|
| 381 | xDelete<IssmDouble>(transition);
|
---|
| 382 | }
|
---|
| 383 | xDelete<IssmDouble*>(transitions);
|
---|
[20138] | 384 | xDelete<int>(transitions_M);
|
---|
| 385 | xDelete<int>(transitions_N);
|
---|
[24481] | 386 | } /*}}}*/
|
---|
| 387 | /*Requested outputs {{{*/
|
---|
[25118] | 388 | iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.solidearth.requested_outputs");
|
---|
[20036] | 389 | if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
|
---|
[25118] | 390 | iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
|
---|
[24481] | 391 | /*}}}*/
|
---|
[20029] | 392 |
|
---|
[19984] | 393 | }/*}}}*/
|
---|
| 394 |
|
---|
| 395 | /*Finite Element Analysis*/
|
---|
| 396 | void SealevelriseAnalysis::Core(FemModel* femmodel){/*{{{*/
|
---|
| 397 | _error_("not implemented");
|
---|
| 398 | }/*}}}*/
|
---|
| 399 | ElementVector* SealevelriseAnalysis::CreateDVector(Element* element){/*{{{*/
|
---|
| 400 | /*Default, return NULL*/
|
---|
| 401 | return NULL;
|
---|
| 402 | }/*}}}*/
|
---|
| 403 | ElementMatrix* SealevelriseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
|
---|
| 404 | _error_("Not implemented");
|
---|
| 405 | }/*}}}*/
|
---|
| 406 | ElementMatrix* SealevelriseAnalysis::CreateKMatrix(Element* element){/*{{{*/
|
---|
| 407 | _error_("not implemented yet");
|
---|
| 408 | }/*}}}*/
|
---|
| 409 | ElementVector* SealevelriseAnalysis::CreatePVector(Element* element){/*{{{*/
|
---|
| 410 | _error_("not implemented yet");
|
---|
| 411 | }/*}}}*/
|
---|
| 412 | void SealevelriseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
|
---|
| 413 | _error_("not implemented yet");
|
---|
| 414 | }/*}}}*/
|
---|
[25317] | 415 | void SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
|
---|
[19984] | 416 | _error_("Not implemented yet");
|
---|
| 417 | }/*}}}*/
|
---|
| 418 | void SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
|
---|
[22955] | 419 | _error_("not implemeneted yet!");
|
---|
[20153] | 420 |
|
---|
[19984] | 421 | }/*}}}*/
|
---|
| 422 | void SealevelriseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
|
---|
| 423 | /*Default, do nothing*/
|
---|
| 424 | return;
|
---|
| 425 | }/*}}}*/
|
---|