source: issm/oecreview/Archive/17984-18295/ISSM-18058-18059.diff@ 18296

Last change on this file since 18296 was 18296, checked in by Mathieu Morlighem, 11 years ago

Added 17984-18295
File size: 10.4 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18058)
4+++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18059)
5@@ -975,7 +975,7 @@
6 void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
7
8 /*return if floating*/
9- if(element->IsFloating())return;
10+ if(element->IsFloating()) return;
11
12 /*Intermediaries*/
13 int domaintype,dim;
14@@ -1067,7 +1067,89 @@
15 _error_("not implemented yet");
16 }/*}}}*/
17 void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
18- _error_("not implemented yet");
19+
20+ /*Intermediaries*/
21+ int domaintype,dim;
22+ Element* basalelement;
23+
24+ /*Get basal element*/
25+ element->FindParam(&domaintype,DomainTypeEnum);
26+ switch(domaintype){
27+ case Domain2DhorizontalEnum:
28+ basalelement = element;
29+ dim = 2;
30+ break;
31+ case Domain2DverticalEnum:
32+ if(!element->IsOnBase()) return;
33+ basalelement = element->SpawnBasalElement();
34+ dim = 1;
35+ break;
36+ case Domain3DEnum:
37+ if(!element->IsOnBase()) return;
38+ basalelement = element->SpawnBasalElement();
39+ dim = 2;
40+ break;
41+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
42+ }
43+
44+ /*Intermediaries*/
45+ IssmDouble Jdet,weight;
46+ IssmDouble thickness,dmudB;
47+ IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
48+ IssmDouble *xyz_list= NULL;
49+
50+ /*Fetch number of vertices for this finite element*/
51+ int numvertices = basalelement->GetNumberOfVertices();
52+
53+ /*Initialize some vectors*/
54+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
55+ IssmDouble* ge = xNew<IssmDouble>(numvertices);
56+ int* vertexpidlist = xNew<int>(numvertices);
57+
58+ /*Retrieve all inputs we will be needing: */
59+ basalelement->GetVerticesCoordinates(&xyz_list);
60+ basalelement->GradientIndexing(&vertexpidlist[0],control_index);
61+ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
62+ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
63+ Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
64+ Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input);
65+ Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input);
66+ Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
67+
68+ /* Start looping on the number of gaussian points: */
69+ Gauss* gauss=basalelement->NewGauss(4);
70+ for(int ig=gauss->begin();ig<gauss->end();ig++){
71+ gauss->GaussPoint(ig);
72+
73+
74+ thickness_input->GetInputValue(&thickness,gauss);
75+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
76+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
77+ adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
78+ adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
79+
80+ basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
81+
82+ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
83+ basalelement->NodalFunctionsP1(basis,gauss);
84+
85+ /*Build gradient vector (actually -dJ/dB): */
86+ for(int i=0;i<numvertices;i++){
87+ ge[i]+=-dmudB*thickness*(
88+ (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
89+ )*Jdet*gauss->weight*basis[i];
90+ _assert_(!xIsNan<IssmDouble>(ge[i]));
91+ }
92+ }
93+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
94+
95+ /*Clean up and return*/
96+ xDelete<IssmDouble>(xyz_list);
97+ xDelete<IssmDouble>(basis);
98+ xDelete<IssmDouble>(ge);
99+ xDelete<int>(vertexpidlist);
100+ delete gauss;
101+ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
102 }/*}}}*/
103 void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
104 _error_("not implemented yet");
105Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
106===================================================================
107--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18058)
108+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18059)
109@@ -283,6 +283,56 @@
110 *pviscosity=viscosity;
111 }
112 /*}}}*/
113+/*FUNCTION Matice::GetViscosity_B {{{*/
114+void Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){
115+ /*From a string tensor and a material object, return viscosity, using Glen's flow law.
116+ (1-D)
117+ viscosity= -----------------------
118+ 2 eps_eff ^[(n-1)/n]
119+
120+ where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
121+ and n the flow law exponent.
122+
123+ If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
124+ return 10^14, initial viscosity.
125+ */
126+
127+ /*output: */
128+ IssmDouble viscosity;
129+
130+ /*Intermediary: */
131+ IssmDouble D=0.,n;
132+
133+ /*Get B and n*/
134+ n=GetN(); _assert_(n>0.);
135+ if(this->isdamaged){
136+ D=GetD();
137+ _assert_(D>=0. && D<1.);
138+ }
139+
140+ if (n==1.){
141+ /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
142+ viscosity=(1.-D)/2.;
143+ }
144+ else{
145+
146+ /*if no strain rate, return maximum viscosity*/
147+ if(eps_eff==0.){
148+ viscosity = 1.e+14/2.;
149+ }
150+
151+ else{
152+ viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n));
153+ }
154+ }
155+
156+ /*Checks in debugging mode*/
157+ if(viscosity<=0) _error_("Negative viscosity");
158+
159+ /*Return: */
160+ *pviscosity=viscosity;
161+}
162+/*}}}*/
163 /*FUNCTION Matice::GetViscosityBar {{{*/
164 void Matice::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){
165 /*From a string tensor and a material object, return viscosity, using Glen's flow law.
166Index: ../trunk-jpl/src/c/classes/Materials/Material.h
167===================================================================
168--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18058)
169+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18059)
170@@ -25,6 +25,7 @@
171 virtual Material* copy2(Element* element)=0;
172 virtual void Configure(Elements* elements)=0;
173 virtual void GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
174+ virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
175 virtual void GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
176 virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
177 virtual void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
178Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
179===================================================================
180--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18058)
181+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18059)
182@@ -54,6 +54,7 @@
183 Material* copy2(Element* element);
184 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
185 void GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
186+ void GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
187 void GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
188 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
189 void GetViscosityDComplement(IssmDouble*, IssmDouble*);
190Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
191===================================================================
192--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18058)
193+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18059)
194@@ -75,6 +75,7 @@
195 Material* copy2(Element* element){_error_("not implemented");};
196 void Configure(Elements* elements);
197 void GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
198+ void GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
199 void GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
200 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
201 void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
202Index: ../trunk-jpl/src/c/classes/Elements/Element.h
203===================================================================
204--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18058)
205+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18059)
206@@ -59,6 +59,7 @@
207 void Echo();
208 void DeepEcho();
209 void DeleteMaterials(void);
210+ void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
211 IssmDouble Divergence(void);
212 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
213 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
214Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
215===================================================================
216--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18058)
217+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18059)
218@@ -212,6 +212,33 @@
219 delete gauss;
220 return divergence;
221 }/*}}}*/
222+void Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
223+
224+ /*Intermediaries*/
225+ IssmDouble dmudB;
226+ IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */
227+ IssmDouble epsilon1d; /* epsilon=[exx]; */
228+ IssmDouble eps_eff;
229+
230+ if(dim==2){
231+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
232+ this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
233+ eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
234+ }
235+ else{
236+ /* eps_eff^2 = 1/2 exx^2*/
237+ this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
238+ eps_eff = sqrt(epsilon1d*epsilon1d/2.);
239+ }
240+
241+ /*Get viscosity*/
242+ material->GetViscosity_B(&dmudB,eps_eff);
243+
244+ /*Assign output pointer*/
245+ *pdmudB=dmudB;
246+
247+}
248+/*}}}*/
249 void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
250 matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
251 }/*}}}*/
Note: See TracBrowser for help on using the repository browser.