Changeset 18059

Timestamp:

05/25/14 18:24:32 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: still working on moving gradient to analysis

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (2 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Materials/Material.h (modified) (1 diff)
• classes/Materials/Matice.cpp (modified) (1 diff)
• classes/Materials/Matice.h (modified) (1 diff)
• classes/Materials/Matpar.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

@@ -976 +976 @@
 
         /*return if floating*/
-        if(element->IsFloating())return;
+        if(element->IsFloating()) return;
 
         /*Intermediaries*/
@@ -1068 +1068 @@
 }/*}}}*/
 void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
-        _error_("not implemented yet");
+
+        /*Intermediaries*/
+        int      domaintype,dim;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DhorizontalEnum:
+                        basalelement = element;
+                        dim          = 2;
+                        break;
+                case Domain2DverticalEnum:
+                        if(!element->IsOnBase()) return;
+                        basalelement = element->SpawnBasalElement();
+                        dim          = 1;
+                        break;
+                case Domain3DEnum:
+                        if(!element->IsOnBase()) return;
+                        basalelement = element->SpawnBasalElement();
+                        dim          = 2;
+                        break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble thickness,dmudB;
+        IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 
+        IssmDouble *xyz_list= NULL;
+
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = basalelement->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
+        IssmDouble* ge            = xNew<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Retrieve all inputs we will be needing: */
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
+        Input* thickness_input = element->GetInput(ThicknessEnum);             _assert_(thickness_input);
+        Input* vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
+        Input* vy_input        = element->GetInput(VyEnum);                    _assert_(vy_input);
+        Input* adjointx_input  = element->GetInput(AdjointxEnum);              _assert_(adjointx_input);
+        Input* adjointy_input  = element->GetInput(AdjointyEnum);              _assert_(adjointy_input);
+        Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+
+                thickness_input->GetInputValue(&thickness,gauss);
+                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
+                adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
+
+                basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctionsP1(basis,gauss);
+
+                /*Build gradient vector (actually -dJ/dB): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+=-dmudB*thickness*(
+                                                (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
+                                                )*Jdet*gauss->weight*basis[i];
+                        _assert_(!xIsNan<IssmDouble>(ge[i]));
+                }
+        }
+        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(ge);
+        xDelete<int>(vertexpidlist);
+        delete gauss;
+        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -213 +213 @@
         return divergence;
 }/*}}}*/
+void       Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble dmudB;
+        IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
+        IssmDouble epsilon1d;   /* epsilon=[exx];    */
+        IssmDouble eps_eff;
+
+         if(dim==2){
+                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
+                 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+                 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
+         }
+         else{
+                 /* eps_eff^2 = 1/2 exx^2*/
+                 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
+                 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
+         }
+
+        /*Get viscosity*/
+        material->GetViscosity_B(&dmudB,eps_eff);
+
+        /*Assign output pointer*/
+        *pdmudB=dmudB;
+
+}
+/*}}}*/
 void       Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
         matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -60 +60 @@
                 void       DeepEcho();
                 void       DeleteMaterials(void);
+                void       dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
                 IssmDouble Divergence(void);
                 void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);

issm/trunk-jpl/src/c/classes/Materials/Material.h

@@ -26 +26 @@
                 virtual void       Configure(Elements* elements)=0;
                 virtual void       GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
+                virtual void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
                 virtual void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
                 virtual void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;

issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

@@ -284 +284 @@
 }
 /*}}}*/
+/*FUNCTION Matice::GetViscosity_B {{{*/
+void  Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){
+        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
+          (1-D)
+          viscosity= -----------------------
+          2 eps_eff ^[(n-1)/n]
+
+          where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
+          and n the flow law exponent.
+
+          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
+          return 10^14, initial viscosity.
+          */
+
+        /*output: */
+        IssmDouble viscosity;
+
+        /*Intermediary: */
+        IssmDouble D=0.,n;
+
+        /*Get B and n*/
+        n=GetN(); _assert_(n>0.);
+        if(this->isdamaged){
+                D=GetD();
+                _assert_(D>=0. && D<1.);
+        }
+
+        if (n==1.){
+                /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
+                viscosity=(1.-D)/2.;
+        }
+        else{
+
+                /*if no strain rate, return maximum viscosity*/
+                if(eps_eff==0.){
+                        viscosity = 1.e+14/2.;
+                }
+
+                else{
+                        viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n));
+                }
+        }
+
+        /*Checks in debugging mode*/
+        if(viscosity<=0) _error_("Negative viscosity");
+
+        /*Return: */
+        *pviscosity=viscosity;
+}
+/*}}}*/
 /*FUNCTION Matice::GetViscosityBar {{{*/
 void  Matice::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){

issm/trunk-jpl/src/c/classes/Materials/Matice.h

@@ -55 +55 @@
                 void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
                 void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
+                void       GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
                 void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
                 void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);

issm/trunk-jpl/src/c/classes/Materials/Matpar.h

@@ -76 +76 @@
                 void       Configure(Elements* elements);
                 void       GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
+                void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
                 void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
                 void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};