source: issm/oecreview/Archive/25834-26739/ISSM-26118-26119.diff@ 26740

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

CHG: added 25834-26739
File size: 7.0 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/cores/cores.h
2===================================================================
3--- ../trunk-jpl/src/c/cores/cores.h (revision 26118)
4+++ ../trunk-jpl/src/c/cores/cores.h (revision 26119)
5@@ -62,9 +62,9 @@
6 void sealevelchange_geometry(FemModel* femmodel);
7 #endif
8 void grd_core(FemModel* femmodel);
9+void grd_core_optim(FemModel* femmodel);
10 void solidearthexternal_core(FemModel* femmodel);
11 void dynstr_core(FemModel* femmodel);
12-SealevelMasks* sealevel_masks(FemModel* femmodel);
13 Vector<IssmDouble>* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
14 Vector<IssmDouble>* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_barystatic,IssmDouble oceanarea);
15 void sealevelchange_core_deformation(Vector<IssmDouble>** pN_radial, Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);
16Index: ../trunk-jpl/src/c/cores/sealevelchange_core.cpp
17===================================================================
18--- ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26118)
19+++ ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26119)
20@@ -16,6 +16,7 @@
21 void TransferSealevel(FemModel* femmodel,int forcingenum);
22 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
23 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
24+SealevelMasks* sealevel_masks(FemModel* femmodel);
25 /*}}}*/
26
27 /*main cores:*/
28@@ -389,19 +390,18 @@
29 }
30 }; /*}}}*/
31
32-void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
33+void grd_core_optim(FemModel* femmodel) { /*{{{*/
34
35 /*variables:{{{*/
36 int nel;
37 BarystaticContributions* barycontrib=NULL;
38+ SealevelMasks* masks=NULL;
39 GenericParam<BarystaticContributions*>* barycontribparam=NULL;
40- IssmDouble barystatic;
41
42 Vector<IssmDouble>* loads=NULL;
43 IssmDouble* allloads=NULL;
44 Vector<IssmDouble>* sealevelloads=NULL;
45 Vector<IssmDouble>* oldsealevelloads=NULL;
46- IssmDouble sealevelloadsaverage;
47 IssmDouble* allsealevelloads=NULL;
48 Vector<IssmDouble>* oceanareas=NULL;
49 IssmDouble oceanarea;
50@@ -412,13 +412,33 @@
51 IssmDouble eps_abs;
52 int step;
53 IssmDouble time;
54-
55- IssmDouble cumbslc;
56- IssmDouble cumbslcice;
57- IssmDouble cumbslchydro;
58+
59+ int modelid,earthid;
60+ int horiz;
61+ int count,frequency,iscoupling;
62+ int grd=0;
63 /*}}}*/
64
65- /*retrieve parameters: */
66+ /*Verbose: */
67+ if(VerboseSolution()) _printf0_(" computing GRD patterns\n");
68+
69+ /*retrieve parameters:{{{*/
70+ femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
71+ femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
72+ femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
73+ /*}}}*/
74+
75+ /*only run if grd was requested, if we are the earth, and we have reached
76+ * the necessary number of time steps dictated by :*/
77+ if(!grd) return;
78+ if(count!=frequency)return;
79+ femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
80+ if(iscoupling){
81+ femmodel->parameters->FindParam(&modelid,ModelIdEnum);
82+ femmodel->parameters->FindParam(&earthid,EarthIdEnum);
83+ if(modelid!=earthid)return;
84+ }
85+ /*retrieve parameters: {{{*/
86 femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
87 barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
88 barycontrib=barycontribparam->GetParameterValue();
89@@ -425,12 +445,17 @@
90 femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
91 femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
92 femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
93-
94+ femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
95+ /*}}}*/
96+
97 /*initialize matrices and vectors:*/
98 femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
99 loads=new Vector<IssmDouble>(nel);
100 sealevelloads=new Vector<IssmDouble>(nel);
101 oceanareas=new Vector<IssmDouble>(nel);
102+
103+ /*call masks core: */
104+ masks=sealevel_masks(femmodel);
105
106 /*buildup loads: */
107 for(Object* & object : femmodel->elements->objects){
108@@ -438,7 +463,7 @@
109 element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
110 }
111
112- //Communicate loads from local to global:
113+ //broadcast loads
114 loads->Assemble(); allloads=loads->ToMPISerial();
115
116 /*convolve loads:*/
117@@ -446,17 +471,16 @@
118 Element* element = xDynamicCast<Element*>(object);
119 element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
120 }
121-
122+ sealevelloads->Assemble();
123+
124 //Get ocean area:
125 oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
126 if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
127
128- //Get sea level loads ocean average:
129- sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
130+ //substract ocean average and barystatic contributionfrom sea level loads:
131+ sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
132
133- //substract ocean average and barystatic contributionfrom sea level loads:
134- barystatic=barycontrib->Total()/oceanarea/rho_water;
135- sealevelloads->Shift(-sealevelloadsaverage+barystatic);
136+ //broadcast sea level loads
137 allsealevelloads=sealevelloads->ToMPISerial();
138
139 bool converged=false;
140@@ -470,12 +494,12 @@
141 element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
142 }
143 sealevelloads->Assemble();
144+
145
146 //substract ocean average and barystatic contribution
147- sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
148- sealevelloads->Shift(-sealevelloadsaverage+barystatic);
149-
150- //broadcast loads
151+ sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
152+
153+ //broadcast sea level loads
154 allsealevelloads=sealevelloads->ToMPISerial();
155
156 //convergence?
157@@ -483,10 +507,25 @@
158 if (converged)break;
159 }
160
161+ /*convolve loads and sea level loads to get the deformation:*/
162+ for(Object* & object : femmodel->elements->objects){
163+ Element* element = xDynamicCast<Element*>(object);
164+ //element->SealevelchangeDeformationConvolution(allsealevelloads, allloads,masks);
165+ }
166+
167+ /*Update bedrock motion and geoid:*/
168+ /*inputs->AXPY(SealevelEnum,SealevelGRD,1);
169+ inputs->AXPY(BedEnum,BedGRD,1);
170+ if(horiz){
171+ inputs->AXPY(BedrockeastEnum,BedEastGRD,1);
172+ inputs->AXPY(BedrocknorthEnum,BedNorthGRD,1);
173+ }*/
174+
175 /*cumulate barystatic contributions and save to results: */
176 barycontrib->Cumulate(femmodel->parameters);
177 barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
178
179+
180 }
181 /*}}}*/
182
Note: See TracBrowser for help on using the repository browser.