source: issm/oecreview/Archive/25834-26739/ISSM-26109-26110.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: 60.9 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h
2===================================================================
3--- ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 26109)
4+++ ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 26110)
5@@ -51,6 +51,7 @@
6 virtual void PointwiseDivide(IssmAbsVec* x,IssmAbsVec* y)=0;
7 virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0;
8 virtual void Pow(doubletype scale_factor)=0;
9+ virtual void Sum(doubletype* pvalue)=0;
10 };
11
12 #endif //#ifndef _ISSM_ABS_VEC_H_
13Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
14===================================================================
15--- ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 26109)
16+++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 26110)
17@@ -593,6 +593,10 @@
18
19 }
20 /*}}}*/
21+ void Sum(doubletype* pvalue){/*{{{*/
22+ _error_("not support yet!");
23+ }
24+ /*}}}*/
25 void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){/*{{{*/
26
27 /*intermediary: */
28Index: ../trunk-jpl/src/c/toolkits/issm/IssmVec.h
29===================================================================
30--- ../trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 26109)
31+++ ../trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 26110)
32@@ -219,6 +219,10 @@
33 vector->Pow(scale_factor);
34 }
35 /*}}}*/
36+ void Sum(doubletype* pvalue){/*{{{*/
37+ vector->Sum(pvalue);
38+ }
39+ /*}}}*/
40 };
41
42 #endif //#ifndef _ISSMVEC_H_
43Index: ../trunk-jpl/src/c/toolkits/objects/Vector.h
44===================================================================
45--- ../trunk-jpl/src/c/toolkits/objects/Vector.h (revision 26109)
46+++ ../trunk-jpl/src/c/toolkits/objects/Vector.h (revision 26110)
47@@ -405,5 +405,16 @@
48 else this->ivector->Pow(scale_factor);
49 }
50 /*}}}*/
51-};
52+void Sum(doubletype* pvalue){ /*{{{*/
53+ _assert_(this);/*{{{*/
54+
55+ if(type==PetscVecType){
56+ #ifdef _HAVE_PETSC_
57+ this->pvector->Sum(pvalue);
58+ #endif
59+ }
60+ else this->ivector->Sum(pvalue);
61+}
62+/*}}}*/
63+}; /*}}}*/
64 #endif //#ifndef _VECTOR_H_
65Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h
66===================================================================
67--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 26109)
68+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 26110)
69@@ -57,6 +57,7 @@
70 IssmDouble Max(void);
71 void Scale(IssmDouble scale_factor);
72 void Pow(IssmDouble scale_factor);
73+ void Sum(IssmDouble* pvalue);
74 void PointwiseDivide(PetscVec* x,PetscVec* y);
75 void PointwiseMult(PetscVec* x,PetscVec* y);
76 IssmDouble Dot(PetscVec* vector);
77Index: ../trunk-jpl/src/c/Makefile.am
78===================================================================
79--- ../trunk-jpl/src/c/Makefile.am (revision 26109)
80+++ ../trunk-jpl/src/c/Makefile.am (revision 26110)
81@@ -96,6 +96,7 @@
82 ./classes/Vertex.cpp \
83 ./classes/Hook.cpp \
84 ./classes/Radar.cpp \
85+ ./classes/BarystaticContributions.cpp \
86 ./classes/ExternalResults/Results.cpp \
87 ./classes/Elements/Element.cpp \
88 ./classes/Elements/Elements.cpp \
89Index: ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
90===================================================================
91--- ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26109)
92+++ ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26110)
93@@ -182,7 +182,9 @@
94 parameters->AddObject(new DoubleMatParam(CumBslcOceanPartitionEnum,bslcocean_partition,npartocean,1));
95 xDelete<IssmDouble>(partitionocean);
96 }
97-
98+ /*New optimized code:*/
99+ BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
100+ parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));
101
102 /*Deal with external multi-model ensembles: {{{*/
103 if(isexternal){
104Index: ../trunk-jpl/src/c/classes/BarystaticContributions.cpp
105===================================================================
106--- ../trunk-jpl/src/c/classes/BarystaticContributions.cpp (nonexistent)
107+++ ../trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26110)
108@@ -0,0 +1,149 @@
109+/*!\file BarystaticContributions.c
110+ * \brief: implementation of the BarystaticContributions object
111+ */
112+
113+/*Include files: {{{*/
114+#ifdef HAVE_CONFIG_H
115+ #include <config.h>
116+#else
117+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
118+#endif
119+
120+#include "./BarystaticContributions.h"
121+#include "../toolkits/toolkits.h"
122+#include "./classes.h"
123+/*}}}*/
124+
125+/*Constructors and destructors:*/
126+BarystaticContributions::BarystaticContributions(IoModel* iomodel ){ /*{{{*/
127+
128+ int nel;
129+
130+ iomodel->FetchData(&nice,"md.solidearth.npartice");
131+ if(nice){
132+ iomodel->FetchData(&pice,&nel,NULL,"md.solidearth.partitionice");
133+ ice=new Vector<IssmDouble>(nice);
134+ cumice=new Vector<IssmDouble>(nice); cumice->Set(0); cumice->Assemble();
135+ }
136+
137+ iomodel->FetchData(&nhydro,"md.solidearth.nparthydro");
138+ if(nhydro){
139+ iomodel->FetchData(&phydro,&nel,NULL,"md.solidearth.partitionhydro");
140+ hydro=new Vector<IssmDouble>(nhydro);
141+ cumhydro=new Vector<IssmDouble>(nhydro); cumhydro->Set(0); cumhydro->Assemble();
142+ }
143+ iomodel->FetchData(&nocean,"md.solidearth.npartocean");
144+ if(nocean){
145+ iomodel->FetchData(&pocean,&nel,NULL,"md.solidearth.partitionocean");
146+ ocean=new Vector<IssmDouble>(nocean);
147+ cumocean=new Vector<IssmDouble>(nocean); cumocean->Set(0); cumocean->Assemble();
148+ }
149+
150+} /*}}}*/
151+BarystaticContributions::~BarystaticContributions(){ /*{{{*/
152+ delete ice; delete cumice;
153+ delete hydro; delete cumhydro;
154+ delete ocean; delete cumocean;
155+ if(nice)xDelete<IssmDouble>(pice);
156+ if(nhydro)xDelete<IssmDouble>(phydro);
157+ if(nocean)xDelete<IssmDouble>(pocean);
158+}; /*}}}*/
159+
160+/*Support routines:*/
161+IssmDouble BarystaticContributions::Total(){ /*{{{*/
162+
163+ IssmDouble sumice,sumhydro,sumocean;
164+
165+ ice->Assemble();
166+ hydro->Assemble();
167+ ocean->Assemble();
168+
169+ ice->Sum(&sumice);
170+ hydro->Sum(&sumhydro);
171+ ocean->Sum(&sumocean);
172+
173+ return sumice+sumhydro+sumocean;
174+
175+} /*}}}*/
176+IssmDouble BarystaticContributions::CumTotal(){ /*{{{*/
177+
178+ IssmDouble sumice,sumhydro,sumocean;
179+
180+ cumice->Assemble();
181+ cumhydro->Assemble();
182+ cumocean->Assemble();
183+
184+ cumice->Sum(&sumice);
185+ cumhydro->Sum(&sumhydro);
186+ cumocean->Sum(&sumocean);
187+
188+
189+ return sumice+sumhydro+sumocean;
190+
191+} /*}}}*/
192+void BarystaticContributions::Cumulate(Parameters* parameters){ /*{{{*/
193+
194+ cumice->AXPY(ice,1);
195+ cumocean->AXPY(ocean,1);
196+ cumhydro->AXPY(hydro,1);
197+
198+
199+} /*}}}*/
200+void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/
201+
202+ int step;
203+ IssmDouble time;
204+ IssmDouble rho_water;
205+
206+ IssmDouble* cumice_serial=NULL;
207+ IssmDouble* cumhydro_serial=NULL;
208+ IssmDouble* cumocean_serial=NULL;
209+
210+ IssmDouble sumice,sumhydro,sumocean;
211+
212+ parameters->FindParam(&step,StepEnum);
213+ parameters->FindParam(&time,TimeEnum);
214+ parameters->FindParam(&rho_water,TimeEnum);
215+
216+ ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean);
217+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time));
218+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time));
219+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time));
220+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time));
221+
222+ cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean);
223+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcEnum,this->CumTotal()/oceanarea/rho_water,step,time));
224+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcIceEnum,sumice/oceanarea/rho_water,step,time));
225+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcHydroEnum,sumhydro/oceanarea/rho_water,step,time));
226+ results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time));
227+
228+ cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;i<nice;i++)cumice_serial[i]=cumice_serial[i]/oceanarea/rho_water;
229+ cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;i<nhydro;i++)cumhydro_serial[i]=cumhydro_serial[i]/oceanarea/rho_water;
230+ cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;i<nocean;i++)cumocean_serial[i]=cumocean_serial[i]/oceanarea/rho_water;
231+
232+ results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time));
233+ results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time));
234+ results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time));
235+
236+ if(IssmComm::GetRank()==0){
237+ xDelete<IssmDouble>(cumice_serial);
238+ xDelete<IssmDouble>(cumhydro_serial);
239+ xDelete<IssmDouble>(cumocean_serial);
240+ }
241+ return;
242+
243+} /*}}}*/
244+void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
245+
246+ int id;
247+
248+ id=reCast<int>(pice[eid]);
249+ ice->SetValue(id,icevalue,ADD_VAL);
250+
251+ id=reCast<int>(phydro[eid]);
252+ hydro->SetValue(id,hydrovalue,ADD_VAL);
253+
254+ id=reCast<int>(pocean[eid]);
255+ ocean->SetValue(id,oceanvalue,ADD_VAL);
256+
257+} /*}}}*/
258Index: ../trunk-jpl/src/c/classes/BarystaticContributions.h
259===================================================================
260--- ../trunk-jpl/src/c/classes/BarystaticContributions.h (nonexistent)
261+++ ../trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26110)
262@@ -0,0 +1,46 @@
263+/*!\file BarystaticContributions.h
264+ * \brief: header file for barystatic contribution object
265+ */
266+
267+#ifndef _BARYSTATICCONTRIBUTIONS_H_
268+#define _BARYSTATICCONTRIBUTIONS_H_
269+
270+/*Headers:*/
271+class IoModel;
272+class Parameters;
273+class Results;
274+template <class doubletype> class Vector;
275+#include "../shared/shared.h"
276+
277+class BarystaticContributions {
278+
279+ public:
280+
281+ Vector<IssmDouble>* ice; //contributions to every ice partition
282+ Vector<IssmDouble>* cumice; //cumulated contributions to every ice partition
283+ int nice; //number of ice partitions
284+ IssmDouble* pice; //ice partition
285+
286+ Vector<IssmDouble>* hydro; //contributions to every hydro partition
287+ Vector<IssmDouble>* cumhydro; //cumulated contributions to every hydro partition
288+ int nhydro; //number of hydro partitions
289+ IssmDouble* phydro; //hydro partition
290+
291+ Vector<IssmDouble>* ocean; //contributions to every ocean partition
292+ Vector<IssmDouble>* cumocean; //cumulated contributions to every ocean partition
293+ int nocean; //number of ocean partitions
294+ IssmDouble* pocean; //ocean partition
295+
296+ /*BarystaticContributions constructors, destructors :*/
297+ BarystaticContributions(IoModel* iomodel );
298+ ~BarystaticContributions();
299+
300+ /*routines:*/
301+ IssmDouble Total();
302+ IssmDouble CumTotal();
303+ void Cumulate(Parameters* parameters);
304+ void Save(Results* results, Parameters* parameters, IssmDouble oceanarea);
305+ void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue);
306+
307+};
308+#endif /* _BARYSTATICCONTRIBUTIONS_H_ */
309Index: ../trunk-jpl/src/c/classes/Elements/Element.h
310===================================================================
311--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26109)
312+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26110)
313@@ -37,6 +37,7 @@
314 template <class doubletype> class Vector;
315 class ElementMatrix;
316 class ElementVector;
317+class BarystaticContributions;
318 /*}}}*/
319
320 class Element: public Object{
321@@ -388,6 +389,12 @@
322 virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0;
323 virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0;
324 virtual void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
325+
326+ virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
327+ virtual void SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
328+ virtual void SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0;
329+ virtual void SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0;
330+ virtual void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
331 #endif
332
333 };
334Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
335===================================================================
336--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 26109)
337+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 26110)
338@@ -226,6 +226,12 @@
339 void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
340 void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
341 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
342+
343+ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
344+ void SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
345+ void SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
346+ void SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
347+ void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
348 #endif
349
350 /*}}}*/
351Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
352===================================================================
353--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 26109)
354+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 26110)
355@@ -181,6 +181,13 @@
356 void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
357 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
358 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
359+
360+ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
361+ void SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
362+ void SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
363+ void SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
364+ void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
365+
366 #endif
367
368 #ifdef _HAVE_DAKOTA_
369Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
370===================================================================
371--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26109)
372+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26110)
373@@ -188,6 +188,12 @@
374 void DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
375 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
376 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
377+ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
378+ void SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
379+ void SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
380+ void SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
381+ void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
382+
383 #endif
384
385 #ifdef _HAVE_DAKOTA_
386Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
387===================================================================
388--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26109)
389+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26110)
390@@ -5431,6 +5431,7 @@
391 /*}}}*/
392 #endif
393 #ifdef _HAVE_SEALEVELCHANGE_
394+//old code
395 IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
396
397 if(masks->isoceanin[this->lid]){
398@@ -5450,7 +5451,7 @@
399
400 }
401 /*}}}*/
402-void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
403+void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
404 /*early return if we are not on an ice cap OR ocean:*/
405 if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
406 dI_list[0] = 0.0; // this is important!!!
407@@ -5544,7 +5545,7 @@
408
409 return;
410 }/*}}}*/
411-void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
412+void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
413
414 masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
415 masks->isoceanin[this->lid]=this->IsOceanInElement();
416@@ -5557,10 +5558,9 @@
417 /*are we not fully grounded: */
418 if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
419 else masks->notfullygrounded[this->lid]=false;
420-
421 }
422 /*}}}*/
423-void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
424+void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
425
426 /*diverse:*/
427 int gsize;
428@@ -5950,7 +5950,7 @@
429 return bslchydro;
430 }
431 /*}}}*/
432-void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
433+void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
434
435 /*diverse:*/
436 int gsize;
437@@ -6000,7 +6000,7 @@
438 return;
439 }
440 /*}}}*/
441-void Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
442+void Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
443
444 /*diverse:*/
445 int gsize,dummy;
446@@ -6042,7 +6042,7 @@
447 return;
448 }
449 /*}}}*/
450-void Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
451+void Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
452
453 /*diverse:*/
454 int gsize;
455@@ -6131,7 +6131,7 @@
456 return;
457 }
458 /*}}}*/
459-void Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
460+void Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
461
462 IssmDouble xyz_list[NUMVERTICES][3];
463
464@@ -6237,8 +6237,8 @@
465 }
466 /*}}}*/
467
468-//new logic
469-void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
470+//new code
471+void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
472
473 /*diverse:*/
474 int nel;
475@@ -6427,7 +6427,7 @@
476 return;
477 }
478 /*}}}*/
479-void Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
480+void Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
481
482 /*diverse:*/
483 IssmDouble area;
484@@ -6435,7 +6435,6 @@
485 IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
486 IssmDouble I,W,BP; //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4)
487 bool notfullygrounded=false;
488- bool scaleoceanarea= false;
489 bool computerigid= false;
490 int glfraction=1;
491 int npartice,nparthydro,npartocean;
492@@ -6468,9 +6467,6 @@
493 #ifdef _ISSM_DEBUG_
494 constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
495 #endif
496- /*barystatic_contribution[0]+=0;
497- barystatic_contribution[1]+=0;
498- barystatic_contribution[2]+=0;*/
499 return;
500 }
501 }
502@@ -6482,9 +6478,6 @@
503 #ifdef _ISSM_DEBUG_
504 this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
505 #endif
506- /*barystatic_contribution[0]+=0;
507- barystatic_contribution[1]+=0;
508- barystatic_contribution[2]+=0;*/
509 return;
510 }
511 }
512@@ -6493,9 +6486,6 @@
513 * hydrology:*/
514 if(!isice && !ishydro){
515 if(!masks->isoceanin[this->lid]){
516- /*barystatic_contribution[0]+=0;
517- barystatic_contribution[1]+=0;
518- barystatic_contribution[2]+=0;*/
519 return;
520 }
521
522@@ -6508,7 +6498,6 @@
523 rho_ice=FindParam(MaterialsRhoIceEnum);
524 rho_water=FindParam(MaterialsRhoSeawaterEnum);
525 this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
526- this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
527 this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
528 this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
529 this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
530@@ -6517,9 +6506,9 @@
531 this->Element::GetInputValue(&area,AreaEnum);
532
533 /*Deal with ice loads if we are on grounded ice:*/
534- if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/
535+ if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){
536
537- /*Compute fraction of the element that is grounded: */
538+ /*Compute fraction of the element that is grounded: {{{*/
539 if(notfullygrounded){
540 IssmDouble xyz_list[NUMVERTICES][3];
541 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
542@@ -6531,6 +6520,7 @@
543 #endif
544 }
545 else phi_ice=1.0;
546+ /*}}}*/
547
548 /*Inform mask: */
549 constant+=100; //1 for ice.
550@@ -6568,15 +6558,13 @@
551 }
552 /*}}}*/
553
554- /*Compute barystatic contribution:*/
555- _assert_(oceanarea>0.);
556- if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
557- bslcice = rho_ice*area*I/(oceanarea*rho_water);
558+ /*Compute barystatic contribution in kg:*/
559+ bslcice = rho_ice*area*I;
560 _assert_(!xIsNan<IssmDouble>(bslcice));
561
562 /*Transfer thickness change into kg/m^2:*/
563 I=I*rho_ice*phi_ice;
564- } /*}}}*/
565+ }
566
567 /*Deal with water loads if we are on ground:*/
568 if(!masks->isfullyfloating[this->lid]){
569@@ -6586,10 +6574,8 @@
570 if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!");
571 deltathickness_input->GetInputAverage(&W);
572
573- /*Compute barystatic component:*/
574- _assert_(oceanarea>0.);
575- if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
576- bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water);
577+ /*Compute barystatic component in kg:*/
578+ bslchydro = rho_freshwater*area*phi_water*W;
579 _assert_(!xIsNan<IssmDouble>(bslchydro));
580
581 /*convert from m to kg/m^2:*/
582@@ -6604,8 +6590,8 @@
583 if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!");
584 bottompressure_change_input->GetInputAverage(&BP);
585
586- /*Compute barystatic component:*/
587- bslcbp = rho_water*area*BP/(oceanarea*rho_water);
588+ /*Compute barystatic component in kg:*/
589+ bslcbp = rho_water*area*BP;
590
591 /*convert from m to kg/m^2:*/
592 BP=BP*rho_water;
593@@ -6612,28 +6598,19 @@
594 }
595
596 /*Plug all loads into total load vector:*/
597- localloads[this->lid]+=I+W+BP;
598+ loads->SetValue(this->sid,I+W+BP,INS_VAL);
599
600- /*Plug bslcice into barystatic contribution vector:*/
601- if(barystatic_contribution_onpartition){
602- int idi=reCast<int>(partition[this->Sid()])+1;
603- int idj=0;
604- idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL);
605- idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL);
606- idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL);
607- }
608+ /*Keep track of barystatic contributions:*/
609+ barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
610
611- barystatic_contribution[0]+=bslcice;
612- barystatic_contribution[1]+=bslchydro;
613- barystatic_contribution[2]+=bslcbp;
614-
615 }
616 /*}}}*/
617-IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/
618+void Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
619
620 /*sal green function:*/
621 IssmDouble* G=NULL;
622- IssmDouble Sealevel[NUMVERTICES]={0,0,0};
623+ IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
624+ IssmDouble oceanaverage,oceanarea=0;
625
626 bool sal = false;
627 int size;
628@@ -6646,15 +6623,115 @@
629
630 for(int i=0;i<NUMVERTICES;i++) {
631 for (int e=0;e<nel;e++){
632- Sealevel[i]+=G[i*nel+e]*loads[e];
633+ SealevelGRD[i]+=G[i*nel+e]*loads[e];
634 }
635 }
636+ }
637
638- this->AddInput(SealevelRSLEnum,Sealevel,P1Enum);
639+ /*compute ocean average over element:*/
640+ OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
641+
642+ /*add ocean average in the global sealevelloads vector:*/
643+ sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
644+ oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
645+
646+ return;
647+} /*}}}*/
648+void Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
649+
650+ bool converged=false;
651+ IssmDouble SealevelGRD[3]={0,0,0};
652+ IssmDouble oceanaverage,oceanarea=0;
653+ int nel;
654+ bool sal = false;
655+ IssmDouble* G=NULL;
656+ int size;
657+
658+ this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
659+ this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
660+
661+ if(sal){
662+ this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
663+
664+ for(int i=0;i<NUMVERTICES;i++) {
665+ for (int e=0;e<nel;e++){
666+ SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
667+ }
668+ }
669 }
670+ OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
671+ newsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
672
673 } /*}}}*/
674+void Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
675
676+ IssmDouble phi=1.0;
677+ bool iscoastline=false;
678+ IssmDouble area;
679+ IssmDouble Sg_avg=0; //output
680+
681+ /*Do we have an ocean?:*/
682+ if(!masks->isoceanin[this->lid]){
683+ *poceanarea=0;
684+ *poceanaverage=0;
685+ }
686+
687+ /*Do we have a coastline?:*/
688+ if(!masks->isfullyfloating[this->lid])iscoastline=true;
689+
690+ if(iscoastline){
691+ IssmDouble xyz_list[NUMVERTICES][3];
692+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
693+ phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]);
694+ }
695+
696+ /*Get area of element:*/
697+ this->Element::GetInputValue(&area,AreaEnum);
698+
699+ /*Average over ocean if there is no coastline, area of the ocean
700+ *is the areaa of the element:*/
701+ if(!iscoastline){
702+
703+ /*Average Sg over vertices:*/
704+ for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[i]/NUMVERTICES;
705+
706+ *poceanaverage=Sg_avg;
707+ *poceanarea=area;
708+ return;
709+ }
710+
711+ /*Average over the ocean only if there is a coastline. Area of the
712+ * ocean will be the fraction times the area of the element:*/
713+ area=phi*area;
714+
715+ IssmDouble total_weight=0;
716+ bool mainlyfloating = true;
717+ int point1;
718+ IssmDouble fraction1,fraction2;
719+
720+ /*Recover portion of element that is grounded*/
721+ this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
722+ //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part
723+ Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2);
724+
725+ /* Start looping on the number of gaussian points and average over these gaussian points: */
726+ total_weight=0;
727+ Sg_avg=0;
728+ while(gauss->next()){
729+ IssmDouble Sg_gauss=0;
730+ TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum);
731+ Sg_avg+=Sg_gauss*gauss->weight;
732+ total_weight+=gauss->weight;
733+ }
734+ Sg_avg=Sg_avg/total_weight;
735+ delete gauss;
736+
737+ *poceanaverage=Sg_avg;
738+ *poceanarea=area;
739+ return;
740+
741+}
742+/*}}}*/
743 #endif
744
745 #ifdef _HAVE_DAKOTA_
746Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
747===================================================================
748--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 26109)
749+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 26110)
750@@ -161,20 +161,22 @@
751 void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz);
752 #endif
753 #ifdef _HAVE_SEALEVELCHANGE_
754+ void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
755+ void SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
756+ IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
757+ IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
758+ void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
759+ void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
760+ void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
761+ void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
762 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
763- void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
764- void SetSealevelMasks(SealevelMasks* masks);
765- void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
766- void SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
767-
768- IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
769- IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
770- void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
771- void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
772- void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
773- void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
774- void SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea);
775- IssmDouble SealevelchangeConvolution(IssmDouble* loads);
776+ void SetSealevelMasks(SealevelMasks* masks);
777+
778+ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
779+ void SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
780+ void SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks);
781+ void SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks);
782+ void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
783 #endif
784 /*}}}*/
785 /*Tria specific routines:{{{*/
786Index: ../trunk-jpl/src/c/classes/FemModel.cpp
787===================================================================
788--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 26109)
789+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 26110)
790@@ -4748,164 +4748,6 @@
791 /*}}}*/
792 #endif
793 #ifdef _HAVE_SEALEVELCHANGE_
794-void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/
795-
796- /*serialized vectors:*/
797- IssmDouble bslcice = 0.;
798- IssmDouble bslcice_cpu = 0.;
799- IssmDouble bslchydro = 0.;
800- IssmDouble bslchydro_cpu = 0.;
801- IssmDouble area = 0.;
802- IssmDouble oceanarea = 0.;
803- IssmDouble oceanarea_cpu = 0.;
804- int bp_compute_fingerprints= 0;
805- bool isoceantransport=false;
806-
807- Vector<IssmDouble>* bslcice_partition=NULL;
808- IssmDouble* bslcice_partition_serial=NULL;
809- IssmDouble* partitionice=NULL;
810- int npartice,nel;
811-
812- Vector<IssmDouble>* bslchydro_partition=NULL;
813- IssmDouble* bslchydro_partition_serial=NULL;
814- IssmDouble* partitionhydro=NULL;
815- bool istws=0;
816- int nparthydro;
817-
818- int npartocean;
819- Vector<IssmDouble>* bslcocean_partition=NULL;
820- IssmDouble* bslcocean_partition_serial=NULL;
821- IssmDouble* partitionocean=NULL;
822-
823- /*Initialize temporary vector that will be used to sum barystatic components
824- * on all local elements, prior to assembly:*/
825- int gsize = this->nodes->NumberOfDofs(GsetEnum);
826- IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
827- int* indices=xNew<int>(gsize);
828- for(int i=0;i<gsize;i++) indices[i]=i;
829-
830- /*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
831- int i = -1;
832- for(Object* & object : this->elements->objects){
833- i +=1;
834- Element* element = xDynamicCast<Element*>(object);
835- element->GetInputValue(&area,AreaEnum);
836- if (masks->isoceanin[i]) oceanarea_cpu += area;
837- }
838- ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
839- ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
840- _assert_(oceanarea>0.);
841-
842- /*Initialize partition vectors to retrieve barystatic contributions: */
843- this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
844- if(npartice){
845- this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum);
846- bslcice_partition= new Vector<IssmDouble>(npartice);
847- }
848-
849- this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
850- if(nparthydro){
851- this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
852- bslchydro_partition= new Vector<IssmDouble>(nparthydro);
853- }
854-
855- this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
856- if(npartocean){
857- this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
858- bslchydro_partition= new Vector<IssmDouble>(npartocean);
859- }
860- /*For later:
861- npartbarystatic=npartice;
862- if(nparthydro>npartbarystatic)npartbarystatic=nparthydro;
863- if(npartocean>npartbarystatic)npartbarystatic=npartocean;
864- bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3);
865-
866- bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0;
867- for(Object* & object : this->elements->objects){
868- Element* element = xDynamicCast<Element*>(object);
869- element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea);
870- }
871- MPI Bcast localloads -> loads
872-
873- for(Object* & object : this->elements->objects){
874- Element* element = xDynamicCast<Element*>(object);
875- element->SealevelchangeConvolution(loads);
876- }
877- */
878-
879-
880-
881-
882-
883- /*Call the barystatic sea level change core for ice : */
884- bslcice_cpu=0;
885- for(Object* & object : this->elements->objects){
886- Element* element = xDynamicCast<Element*>(object);
887- bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea);
888- }
889-
890- /*Call the barystatic sea level change core for hydro: */
891- bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0.
892- this->parameters->FindParam(&istws,TransientIshydrologyEnum);
893- if(istws){
894- for(int i=0;i<elements->Size();i++){
895- Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
896- bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea);
897- }
898- }
899-
900- /*Call the barystatic sea level change core for bottom pressures: */
901- this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
902- this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum);
903- if(bp_compute_fingerprints && isoceantransport){
904- for(int i=0;i<elements->Size();i++){
905- Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
906- element->SealevelchangeBarystaticBottomPressure(RSLgi,masks);
907- }
908- }
909-
910- /*Plug values once and assemble: */
911- pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
912- pRSLgi->Assemble();
913-
914- /*Sum all barystatic components from all cpus:*/
915- ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
916- ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
917- _assert_(!xIsNan<IssmDouble>(bslcice));
918-
919- ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
920- ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
921- _assert_(!xIsNan<IssmDouble>(bslchydro));
922-
923- /*Take care of partition vectors:*/
924- if(bslcice_partition){
925- bslcice_partition->Assemble();
926- bslcice_partition_serial=bslcice_partition->ToMPISerial();
927- }
928- if(bslchydro_partition){
929- bslchydro_partition->Assemble();
930- bslchydro_partition_serial=bslchydro_partition->ToMPISerial();
931- }
932-
933-
934- /*Free ressources:*/
935- xDelete<int>(indices);
936- xDelete<IssmDouble>(RSLgi);
937- if(bslchydro_partition)delete bslchydro_partition;
938- if(bslcice_partition)delete bslcice_partition;
939- if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
940- if(partitionice)xDelete<IssmDouble>(partitionice);
941-
942- /*Assign output pointers:*/
943- *poceanarea = oceanarea;
944- *pbslcice = bslcice;
945- *pbslchydro = bslchydro;
946- *pbslc=bslchydro+bslcice;
947- *pbslcice_partition=bslcice_partition_serial;
948- *pbslchydro_partition=bslchydro_partition_serial;
949-
950-}
951-/*}}}*/
952 void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, SealevelMasks* masks, bool verboseconvolution){/*{{{*/
953
954 /*serialized vectors:*/
955Index: ../trunk-jpl/src/c/classes/classes.h
956===================================================================
957--- ../trunk-jpl/src/c/classes/classes.h (revision 26109)
958+++ ../trunk-jpl/src/c/classes/classes.h (revision 26110)
959@@ -18,6 +18,7 @@
960 #include "./Massfluxatgate.h"
961 #include "./Misfit.h"
962 #include "./SealevelMasks.h"
963+#include "./BarystaticContributions.h"
964 #include "./Nodalvalue.h"
965 #include "./Numberedcostfunction.h"
966 #include "./Cfsurfacesquare.h"
967Index: ../trunk-jpl/src/c/cores/cores.h
968===================================================================
969--- ../trunk-jpl/src/c/cores/cores.h (revision 26109)
970+++ ../trunk-jpl/src/c/cores/cores.h (revision 26110)
971@@ -80,10 +80,7 @@
972 void WriteLockFile(char* filename);
973 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
974 void PrintBanner(void);
975-void TransferForcing(FemModel* femmodel,int forcingenum);
976-void TransferSealevel(FemModel* femmodel,int forcingenum);
977 void EarthMassTransport(FemModel* femmodel);
978-void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
979
980 //solution configuration
981 void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);
982Index: ../trunk-jpl/src/c/cores/sealevelchange_core.cpp
983===================================================================
984--- ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26109)
985+++ ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26110)
986@@ -11,6 +11,12 @@
987 #include "../shared/shared.h"
988 #include "../modules/modules.h"
989 #include "../solutionsequences/solutionsequences.h"
990+/*support routines local definitions:{{{*/
991+void TransferForcing(FemModel* femmodel,int forcingenum);
992+void TransferSealevel(FemModel* femmodel,int forcingenum);
993+void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
994+IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
995+/*}}}*/
996
997 /*main cores:*/
998 void sealevelchange_core(FemModel* femmodel){ /*{{{*/
999@@ -383,6 +389,107 @@
1000 }
1001 }; /*}}}*/
1002
1003+void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
1004+
1005+ /*variables:{{{*/
1006+ int nel;
1007+ BarystaticContributions* barycontrib=NULL;
1008+ GenericParam<BarystaticContributions*>* barycontribparam=NULL;
1009+ IssmDouble barystatic;
1010+
1011+ Vector<IssmDouble>* loads=NULL;
1012+ IssmDouble* allloads=NULL;
1013+ Vector<IssmDouble>* sealevelloads=NULL;
1014+ Vector<IssmDouble>* oldsealevelloads=NULL;
1015+ IssmDouble sealevelloadsaverage;
1016+ IssmDouble* allsealevelloads=NULL;
1017+ Vector<IssmDouble>* oceanareas=NULL;
1018+ IssmDouble oceanarea;
1019+ bool scaleoceanarea=false;
1020+ IssmDouble rho_water;
1021+
1022+ IssmDouble eps_rel;
1023+ IssmDouble eps_abs;
1024+ int step;
1025+ IssmDouble time;
1026+
1027+ IssmDouble cumbslc;
1028+ IssmDouble cumbslcice;
1029+ IssmDouble cumbslchydro;
1030+ /*}}}*/
1031+
1032+ /*retrieve parameters: */
1033+ femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
1034+ barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
1035+ barycontrib=barycontribparam->GetParameterValue();
1036+ femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
1037+ femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
1038+ femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
1039+
1040+ /*initialize matrices and vectors:*/
1041+ femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
1042+ loads=new Vector<IssmDouble>(nel);
1043+ sealevelloads=new Vector<IssmDouble>(nel);
1044+ oceanareas=new Vector<IssmDouble>(nel);
1045+
1046+ /*buildup loads: */
1047+ for(Object* & object : femmodel->elements->objects){
1048+ Element* element = xDynamicCast<Element*>(object);
1049+ element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
1050+ }
1051+
1052+ //Communicate loads from local to global:
1053+ loads->Assemble(); allloads=loads->ToMPISerial();
1054+
1055+ /*convolve loads:*/
1056+ for(Object* & object : femmodel->elements->objects){
1057+ Element* element = xDynamicCast<Element*>(object);
1058+ element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
1059+ }
1060+
1061+ //Get ocean area:
1062+ oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
1063+ if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
1064+
1065+ //Get sea level loads ocean average:
1066+ sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
1067+
1068+ //substract ocean average and barystatic contributionfrom sea level loads:
1069+ barystatic=barycontrib->Total()/oceanarea/rho_water;
1070+ sealevelloads->Shift(-sealevelloadsaverage+barystatic);
1071+ allsealevelloads=sealevelloads->ToMPISerial();
1072+
1073+ bool converged=false;
1074+ for(;;){
1075+
1076+ oldsealevelloads=sealevelloads->Duplicate();
1077+
1078+ /*convolve load and sealevel loads on oceans:*/
1079+ for(Object* & object : femmodel->elements->objects){
1080+ Element* element = xDynamicCast<Element*>(object);
1081+ element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
1082+ }
1083+ sealevelloads->Assemble();
1084+
1085+ //substract ocean average and barystatic contribution
1086+ sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
1087+ sealevelloads->Shift(-sealevelloadsaverage+barystatic);
1088+
1089+ //broadcast loads
1090+ allsealevelloads=sealevelloads->ToMPISerial();
1091+
1092+ //convergence?
1093+ slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
1094+ if (converged)break;
1095+ }
1096+
1097+ /*cumulate barystatic contributions and save to results: */
1098+ barycontrib->Cumulate(femmodel->parameters);
1099+ barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
1100+
1101+}
1102+/*}}}*/
1103+
1104 //Geometry:
1105 SealevelMasks* sealevel_masks(FemModel* femmodel) { /*{{{*/
1106
1107@@ -1061,3 +1168,11 @@
1108 *pconverged=converged;
1109
1110 } /*}}}*/
1111+IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
1112+ IssmDouble sealevelloadsaverage;
1113+ Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
1114+ sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
1115+ sealevelloadsvolume->Sum(&sealevelloadsaverage);
1116+ delete sealevelloadsvolume;
1117+ return sealevelloadsaverage/oceanarea;
1118+} /*}}}*/
1119Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim
1120===================================================================
1121--- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26109)
1122+++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26110)
1123@@ -746,6 +746,7 @@
1124 syn keyword cConstant BslcEnum
1125 syn keyword cConstant BslcIceEnum
1126 syn keyword cConstant BslcHydroEnum
1127+syn keyword cConstant BslcOceanEnum
1128 syn keyword cConstant BslcRateEnum
1129 syn keyword cConstant GmtslcEnum
1130 syn keyword cConstant SealevelRSLBarystaticEnum
1131@@ -1070,6 +1071,7 @@
1132 syn keyword cConstant BalancevelocitySolutionEnum
1133 syn keyword cConstant BasalforcingsIsmip6Enum
1134 syn keyword cConstant BasalforcingsPicoEnum
1135+syn keyword cConstant BarystaticContributionsEnum
1136 syn keyword cConstant BeckmannGoosseFloatingMeltRateEnum
1137 syn keyword cConstant BedSlopeSolutionEnum
1138 syn keyword cConstant BoolExternalResultEnum
1139@@ -1432,6 +1434,7 @@
1140 syn keyword cType AmrBamg
1141 syn keyword cType AmrNeopz
1142 syn keyword cType ArrayInput
1143+syn keyword cType BarystaticContributions
1144 syn keyword cType BoolInput
1145 syn keyword cType BoolParam
1146 syn keyword cType Cfdragcoeffabsgrad
1147Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
1148===================================================================
1149--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26109)
1150+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26110)
1151@@ -742,6 +742,7 @@
1152 BslcEnum,
1153 BslcIceEnum,
1154 BslcHydroEnum,
1155+ BslcOceanEnum,
1156 BslcRateEnum,
1157 GmtslcEnum,
1158 SealevelRSLBarystaticEnum,
1159@@ -1069,6 +1070,7 @@
1160 BalancevelocitySolutionEnum,
1161 BasalforcingsIsmip6Enum,
1162 BasalforcingsPicoEnum,
1163+ BarystaticContributionsEnum,
1164 BeckmannGoosseFloatingMeltRateEnum,
1165 BedSlopeSolutionEnum,
1166 BoolExternalResultEnum,
1167Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
1168===================================================================
1169--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26109)
1170+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26110)
1171@@ -748,6 +748,7 @@
1172 case BslcEnum : return "Bslc";
1173 case BslcIceEnum : return "BslcIce";
1174 case BslcHydroEnum : return "BslcHydro";
1175+ case BslcOceanEnum : return "BslcOcean";
1176 case BslcRateEnum : return "BslcRate";
1177 case GmtslcEnum : return "Gmtslc";
1178 case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic";
1179@@ -1072,6 +1073,7 @@
1180 case BalancevelocitySolutionEnum : return "BalancevelocitySolution";
1181 case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6";
1182 case BasalforcingsPicoEnum : return "BasalforcingsPico";
1183+ case BarystaticContributionsEnum : return "BarystaticContributions";
1184 case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate";
1185 case BedSlopeSolutionEnum : return "BedSlopeSolution";
1186 case BoolExternalResultEnum : return "BoolExternalResult";
1187Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
1188===================================================================
1189--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26109)
1190+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26110)
1191@@ -766,6 +766,7 @@
1192 else if (strcmp(name,"Bslc")==0) return BslcEnum;
1193 else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
1194 else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
1195+ else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
1196 else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
1197 else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
1198 else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
1199@@ -873,11 +874,11 @@
1200 else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
1201 else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
1202 else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
1203- else if (strcmp(name,"SmbT")==0) return SmbTEnum;
1204 else stage=8;
1205 }
1206 if(stage==8){
1207- if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
1208+ if (strcmp(name,"SmbT")==0) return SmbTEnum;
1209+ else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
1210 else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
1211 else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
1212 else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
1213@@ -996,11 +997,11 @@
1214 else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
1215 else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
1216 else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
1217- else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
1218 else stage=9;
1219 }
1220 if(stage==9){
1221- if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
1222+ if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
1223+ else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
1224 else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
1225 else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
1226 else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
1227@@ -1096,6 +1097,7 @@
1228 else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
1229 else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
1230 else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
1231+ else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum;
1232 else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
1233 else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
1234 else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
1235@@ -1118,12 +1120,12 @@
1236 else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
1237 else if (strcmp(name,"Closed")==0) return ClosedEnum;
1238 else if (strcmp(name,"Colinear")==0) return ColinearEnum;
1239- else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
1240- else if (strcmp(name,"Contact")==0) return ContactEnum;
1241 else stage=10;
1242 }
1243 if(stage==10){
1244- if (strcmp(name,"Contour")==0) return ContourEnum;
1245+ if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
1246+ else if (strcmp(name,"Contact")==0) return ContactEnum;
1247+ else if (strcmp(name,"Contour")==0) return ContourEnum;
1248 else if (strcmp(name,"Contours")==0) return ContoursEnum;
1249 else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
1250 else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum;
1251@@ -1241,12 +1243,12 @@
1252 else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
1253 else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
1254 else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
1255- else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
1256- else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
1257 else stage=11;
1258 }
1259 if(stage==11){
1260- if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
1261+ if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
1262+ else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
1263+ else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
1264 else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
1265 else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
1266 else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
1267@@ -1364,12 +1366,12 @@
1268 else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
1269 else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
1270 else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
1271- else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
1272- else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
1273 else stage=12;
1274 }
1275 if(stage==12){
1276- if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
1277+ if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
1278+ else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
1279+ else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
1280 else if (strcmp(name,"Scaled")==0) return ScaledEnum;
1281 else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
1282 else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
1283Index: ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h
1284===================================================================
1285--- ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 26109)
1286+++ ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 26110)
1287@@ -323,6 +323,14 @@
1288
1289 }
1290 /*}}}*/
1291+ void Sum(doubletype* pvalue){/*{{{*/
1292+
1293+ doubletype value=0;
1294+ int i;
1295+ for(i=0;i<this->M;i++)value+=this->vector[i];
1296+ *pvalue=value;
1297+ }
1298+ /*}}}*/
1299
1300 };
1301 #endif //#ifndef _ISSM_SEQ_VEC_H_
1302Index: ../trunk-jpl/src/c/toolkits/objects/Matrix.h
1303===================================================================
1304--- ../trunk-jpl/src/c/toolkits/objects/Matrix.h (revision 26109)
1305+++ ../trunk-jpl/src/c/toolkits/objects/Matrix.h (revision 26110)
1306@@ -281,6 +281,22 @@
1307 return output;
1308 }
1309 /*}}}*/
1310+ doubletype* ToMPISerial(void){/*{{{*/
1311+
1312+ doubletype* output=NULL;
1313+
1314+ if(type==PetscMatType){
1315+ #ifdef _HAVE_PETSC_
1316+ output=this->pmatrix->ToMPISerial();
1317+ #endif
1318+ }
1319+ else{
1320+ _error_("not implemented yet!");
1321+ }
1322+
1323+ return output;
1324+ }
1325+ /*}}}*/
1326 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/
1327
1328 if(type==PetscMatType){
1329@@ -306,10 +322,8 @@
1330
1331 }
1332 /*}}}*/
1333- /*
1334- * sets all values to 0 but keeps the structure of a sparse matrix
1335- */
1336 void SetZero(void) {/*{{{*/
1337+ // sets all values to 0 but keeps the structure of a sparse matrix
1338 if(type==PetscMatType){
1339 #ifdef _HAVE_PETSC_
1340 this->pmatrix->SetZero();
1341Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp
1342===================================================================
1343--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 26109)
1344+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 26110)
1345@@ -248,6 +248,13 @@
1346
1347 }
1348 /*}}}*/
1349+void PetscVec::Sum(IssmDouble* pvalue){/*{{{*/
1350+
1351+ _assert_(this->vector);
1352+ VecSum(this->vector,pvalue);
1353+
1354+}
1355+/*}}}*/
1356 IssmDouble PetscVec::Dot(PetscVec* input){/*{{{*/
1357
1358 IssmDouble dot;
1359Index: ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h
1360===================================================================
1361--- ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 26109)
1362+++ ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 26110)
1363@@ -31,6 +31,7 @@
1364 void PetscOptionsDetermineSolverType(int* psolver_type);
1365 void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
1366 void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
1367+void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
1368 Vec SerialToVec(double* vector,int vector_size);
1369 InsertMode ISSMToPetscInsertMode(InsMode mode);
1370 NormType ISSMToPetscNormMode(NormMode mode);
1371Index: ../trunk-jpl/test/NightlyRun/test2004.m
1372===================================================================
1373--- ../trunk-jpl/test/NightlyRun/test2004.m (revision 26109)
1374+++ ../trunk-jpl/test/NightlyRun/test2004.m (revision 26110)
1375@@ -133,6 +133,10 @@
1376
1377 disp(' reading bedrock');
1378 md.geometry.bed=-ones(md.mesh.numberofvertices,1);
1379+ md.geometry.base=md.geometry.bed;
1380+ md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
1381+ md.geometry.surface=md.geometry.bed+md.geometry.thickness;
1382+
1383 end % }}}
1384 %Slc: {{{
1385 if bas.iscontinentany('antarctica'),
1386@@ -166,11 +170,11 @@
1387 md.masstransport.spcthickness=mean(delHAIS(md.mesh.elements),2)/100;
1388 end
1389
1390- md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
1391+ md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
1392
1393- md.dsl.global_average_thermosteric_sea_level_change=[0;0];
1394- md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
1395- md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
1396+ md.dsl.global_average_thermosteric_sea_level=[0;0];
1397+ md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
1398+ md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
1399
1400 end %}}}
1401 % material properties: {{{
1402@@ -312,6 +316,9 @@
1403 %geometry: {{{
1404 di=md.materials.rho_ice/md.materials.rho_water;
1405 md.geometry.bed=-ones(md.mesh.numberofvertices,1);
1406+ md.geometry.base=md.geometry.bed;
1407+ md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
1408+ md.geometry.surface=md.geometry.bed+md.geometry.thickness;
1409 % }}}
1410 %materials: {{{
1411 md.materials=materials('hydro');
1412@@ -345,6 +352,9 @@
1413 sl.transfer('mask.ice_levelset');
1414 sl.transfer('mask.ocean_levelset');
1415 sl.transfer('geometry.bed');
1416+sl.transfer('geometry.surface');
1417+sl.transfer('geometry.thickness');
1418+sl.transfer('geometry.base');
1419 sl.transfer('mesh.lat');
1420 sl.transfer('mesh.long');
1421 sl.transfer('masstransport.spcthickness');
1422@@ -399,6 +409,16 @@
1423 md.transient.ismasstransport=1;
1424 md.transient.isslc=1;
1425
1426+%Initializations:
1427+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
1428+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
1429+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
1430+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
1431+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
1432+md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1);
1433+md.initialization.dsl=zeros(md.mesh.numberofvertices,1);
1434+md.initialization.str=0;
1435+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
1436
1437 %max number of iterations reverted back to 10 (i.e. the original default value)
1438 md.solidearth.settings.maxiter=10;
Note: See TracBrowser for help on using the repository browser.