Changeset 19002

Timestamp:

01/09/15 08:48:50 (10 years ago)

Author:

bdef

Message:

Adding the inversion capability for hydro related friction law

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (4 diffs)
• analyses/AdjointHorizAnalysis.h (modified) (1 diff)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Inputs/TriaInput.h (modified) (1 diff)
• classes/Loads/Friction.cpp (modified) (8 diffs)
• modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (modified) (1 diff)

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

@@ -26 +26 @@
 
 /*Finite Element Analysis*/
-void           AdjointHorizAnalysis::Core(FemModel* femmodel){/*{{{*/
+void AdjointHorizAnalysis::Core(FemModel* femmodel){/*{{{*/
         _error_("not implemented");
 }/*}}}*/
@@ -1074 +1074 @@
         /*Check that control_type is supported*/
         if(control_type!=MaterialsRheologyBbarEnum && 
-                control_type!=FrictionCoefficientEnum   && 
+                control_type!=FrictionCoefficientEnum   &&
+                control_type!=FrictionAsEnum   && 
                 control_type!=DamageDbarEnum            &&
                 control_type!=MaterialsRheologyBEnum){
@@ -1101 +1102 @@
                                 case HOApproximationEnum:  GradientJDragHO( element,gradient,control_index); break;
                                 case FSApproximationEnum:  GradientJDragFS( element,gradient,control_index); break;
+                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
+                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
+                        }
+                        break;
+                case FrictionAsEnum:
+                        switch(approximation){
+                                case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_index); break;
+                                case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_index); break;
+                                case HOApproximationEnum:  GradientJDragHydroHO( element,gradient,control_index); break;
+                                case FSApproximationEnum:  GradientJDragHydroFS( element,gradient,control_index); break;
                                 case NoneApproximationEnum: /*Gradient is 0*/                    break;
                                 default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
@@ -1886 +1897 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+
+void AdjointHorizAnalysis::GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*return if floating or not on bed (gradient is 0)*/
+        if(element->IsFloating()) return;
+        if(!element->IsOnBase()) return;
+
+        /*Intermediaries*/
+        int        domaintype,dim;
+        IssmDouble Jdet,weight;
+        IssmDouble drag,dalpha2dk,normal[3];
+        IssmDouble vx,vy,vz,lambda,mu,xi;
+        IssmDouble *xyz_list_base= NULL;
+
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
+        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /* get domaintype */
+        element->FindParam(&domaintype,DomainTypeEnum);
+
+        /*Build friction element, needed later: */
+        if(domaintype!=Domain2DverticalEnum) dim=3;
+        else dim=2;
+        Friction* friction=new Friction(element,dim);
+
+        /*Retrieve all inputs we will be needing: */
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->GradientIndexing(&vertexpidlist[0],control_index);
+        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* vz_input        = NULL;
+        Input* adjointz_input  = NULL;
+        if(domaintype!=Domain2DverticalEnum){
+                vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
+                adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
+        }
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGaussBase(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                adjointx_input->GetInputValue(&lambda, gauss);
+                adjointy_input->GetInputValue(&mu, gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                if(domaintype!=Domain2DverticalEnum){
+                        adjointz_input->GetInputValue(&xi    ,gauss);
+                        vz_input->GetInputValue(&vz,gauss);
+                }
+
+                friction->GetAlphaComplement(&dalpha2dk,gauss);
+                element->NormalBase(&normal[0],xyz_list_base);
+
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctionsP1(basis,gauss);
+
+                /*Build gradient vector (actually -dJ/dk): */
+                if(domaintype!=Domain2DverticalEnum){
+                        for(int i=0;i<numvertices;i++){
+                                ge[i]+=(
+                                                        -lambda*(dalpha2dk*(vx - vz*normal[0]*normal[2]))
+                                                        -mu    *(dalpha2dk*(vy - vz*normal[1]*normal[2]))
+                                                        -xi    *(dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
+                                                 )*Jdet*gauss->weight*basis[i];
+                                _assert_(!xIsNan<IssmDouble>(ge[i]));
+                        }
+                }
+                else{
+                        for(int i=0;i<numvertices;i++){
+                                ge[i]+=(
+                                                        -lambda*dalpha2dk*vx
+                                                        -mu    *dalpha2dk*vy
+                                                 )*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_base);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(ge);
+        xDelete<int>(vertexpidlist);
+        delete gauss;
+        delete friction;
+}/*}}}*/
+void           AdjointHorizAnalysis::GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Same as SSA*/
+        return this->GradientJDragSSA(element,gradient,control_index);
+}/*}}}*/
+void           AdjointHorizAnalysis::GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*return if floating or not on bed (gradient is 0)*/
+        if(element->IsFloating()) return;
+        if(!element->IsOnBase()) return;
+
+        /*Intermediaries*/
+        IssmDouble Jdet,weight;
+        IssmDouble drag,dalpha2dk;
+        IssmDouble vx,vy,lambda,mu;
+        IssmDouble *xyz_list_base= NULL;
+
+        int      domaintype,dim;
+        element->FindParam(&domaintype,DomainTypeEnum);
+        /*Fetch number of vertices for this finite element*/
+        int numvertices = element->GetNumberOfVertices();
+
+        /*Initialize some vectors*/
+        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
+        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Build friction element, needed later: */
+        if(domaintype!=Domain2DverticalEnum) dim=3;
+        else dim=2;
+        Friction* friction=new Friction(element,dim);
+
+        /*Retrieve all inputs we will be needing: */
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->GradientIndexing(&vertexpidlist[0],control_index);
+        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
+        Input* vy_input        = NULL;
+        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
+        Input* adjointy_input  = NULL;
+        if(domaintype!=Domain2DverticalEnum){
+                vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
+                adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
+        }
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGaussBase(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                adjointx_input->GetInputValue(&lambda, gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                if(domaintype!=Domain2DverticalEnum){
+                        adjointy_input->GetInputValue(&mu, gauss);
+                        vy_input->GetInputValue(&vy,gauss);
+                }
+
+                friction->GetAlphaComplement(&dalpha2dk,gauss);
+
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctionsP1(basis,gauss);
+
+                /*Build gradient vector (actually -dJ/dD): */
+                for(int i=0;i<numvertices;i++){
+                        if(domaintype!=Domain2DverticalEnum) ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
+                        else ge[i]+=-dalpha2dk*(lambda*vx)*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_base);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(ge);
+        xDelete<int>(vertexpidlist);
+        delete gauss;
+        delete friction;
+}/*}}}*/
+
+void AdjointHorizAnalysis::GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*return if floating (gradient is 0)*/
+        if(element->IsFloating()) return;
+
+        /*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 dalpha2dk;
+        IssmDouble vx,vy,lambda,mu;
+        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            = xNewZeroInit<IssmDouble>(numvertices);
+        int*        vertexpidlist = xNew<int>(numvertices);
+
+        /*Build friction element, needed later: */
+        Friction* friction=new Friction(basalelement,dim);
+
+        /*Retrieve all inputs we will be needing: */
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
+        Input* vx_input        = basalelement->GetInput(VxEnum);          _assert_(vx_input);
+        Input* vy_input        = basalelement->GetInput(VyEnum);          _assert_(vy_input);
+        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);    _assert_(adjointx_input);
+        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);    _assert_(adjointy_input);
+
+
+
+        IssmDouble  q_exp;
+        IssmDouble  C_param;
+        IssmDouble  As;
+        IssmDouble  Neff;
+        IssmDouble  n;
+        IssmDouble  alpha;
+        IssmDouble  Chi,Gamma;
+        IssmDouble  vz,vmag;
+        IssmDouble  Uder;
+
+        /*Recover parameters: */
+        Input* qinput = basalelement->GetInput(FrictionQEnum);
+        Input* cinput = basalelement->GetInput(FrictionCEnum);
+        Input* Asinput = basalelement->GetInput(FrictionAsEnum);
+        Input* Ninput = basalelement->GetInput(FrictionEffectivePressureEnum);
+        Input* nInput =basalelement->GetInput(MaterialsRheologyNEnum);
+        
+        /* 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);
+
+                adjointx_input->GetInputValue(&lambda, gauss);
+                adjointy_input->GetInputValue(&mu, gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+
+                friction->GetAlphaComplement(&dalpha2dk,gauss);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctionsP1(basis,gauss);
+
+                /*Dealing with dalpha/du*/
+                qinput->GetInputValue(&q_exp,gauss);
+                cinput->GetInputValue(&C_param,gauss);
+                Asinput->GetInputValue(&As,gauss);
+                Ninput->GetInputValue(&Neff,gauss);
+                nInput->GetInputValue(&n,gauss);
+
+                if (q_exp==1){
+                        alpha=1;
+                }
+                else{
+                        alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
+                }
+                Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
+                Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
+                
+                Uder =Neff*C_param/(vmag*vmag*n) *
+                        (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) -
+                         n* pow(Gamma,1./n));
+                
+                /*Build gradient vector (actually -dJ/dD): */
+                for(int i=0;i<numvertices;i++){
+                        ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*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;
+        delete friction;
+        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+}/*}}}*/
+
 void           AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 

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

@@ -50 +50 @@
                 void           GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void           GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void           GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void           GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void           GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void           InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

@@ -1148 +1148 @@
                                 name==WatercolumnEnum || 
                                 name==FrictionCoefficientEnum ||
+                                name==FrictionAsEnum ||
                                 name==MaskGroundediceLevelsetEnum ||
                                 name==MaskIceLevelsetEnum ||

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

@@ -1574 +1574 @@
                                 case ThicknessEnum:
                                 case FrictionCoefficientEnum:
+                          case FrictionAsEnum:  
                                 case MaterialsRheologyBEnum:
                                         if(iomodel->Data(control)){

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

@@ -1774 +1774 @@
                                 case BalancethicknessOmegaEnum:
                                 case FrictionCoefficientEnum:
-                                case MaterialsRheologyBEnum:
+                          case FrictionAsEnum:
+                          case MaterialsRheologyBEnum:
                                         if(iomodel->Data(control)){
                                                 for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[tria_vertex_ids[j]-1];

issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

@@ -49 +49 @@
                 void GetInputValue(bool* pvalue){_error_("not implemented yet");}
                 void GetInputValue(int* pvalue){_error_("not implemented yet");}
-                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}
+                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet (Input is "<<EnumToStringx(this->enum_type)<<")");}//{_error_("not implemented yet");}
                 void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
                 void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -126 +126 @@
         IssmDouble  n;
         IssmDouble  alpha;
-        IssmDouble  Chi;
+        IssmDouble  Chi,Gamma;
         IssmDouble  vx,vy,vz,vmag;
         IssmDouble  alpha_complement;
@@ -133 +133 @@
         element->GetInputValue(&q_exp,FrictionQEnum);
         element->GetInputValue(&C_param,FrictionCEnum);
-        element->GetInputValue(&As,FrictionAsEnum);
-
+
+        element->GetInputValue(&As,gauss,FrictionAsEnum);
         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
         element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
@@ -160 +160 @@
                         _error_("not supported");
         }
+        //      vmag=100./(3600.*24.*365.);
+
         if (q_exp==1){
                 alpha=1;
@@ -166 +168 @@
                 alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
         }
-        Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
-
+        Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
+        Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
         /*Check to prevent dividing by zero if vmag==0*/
         if(vmag==0.) alpha_complement=0.;
-        else if(Neff==0.) alpha_complement=0.;
-        else alpha_complement=-(C_param*Neff/(vmag*n)) * pow((Chi/(alpha*pow(Chi,q_exp)+1)),((1-n)/n)) *(Chi/(As*(alpha*pow(Chi,q_exp)+1)))-(alpha*q_exp*pow(Chi,q_exp+1)/(As*(alpha*pow(Chi,q_exp)+1)));
+        else    if(Neff==0.) alpha_complement=0.;
+        else    alpha_complement=-(C_param*Neff/(n*vmag)) *
+                                        pow(Gamma,((1.-n)/n)) *
+                                        (Gamma/As - (alpha*q_exp*pow(Chi,q_exp-1.)* Gamma * Gamma/As));
+
         _assert_(!xIsNan<IssmDouble>(alpha_complement));
         /*Assign output pointers:*/
@@ -221 +226 @@
 
         IssmDouble  alpha;
-        IssmDouble  Chi;
+        IssmDouble  Chi,Gamma;
 
         IssmDouble  vx,vy,vz,vmag;
@@ -229 +234 @@
         element->GetInputValue(&q_exp,FrictionQEnum);
         element->GetInputValue(&C_param,FrictionCEnum);
-        element->GetInputValue(&As,FrictionAsEnum);
+        element->GetInputValue(&As,gauss,FrictionAsEnum);
 
         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
@@ -256 +261 @@
         }
 
+        //      vmag=100./(3600.*24.*365.);
         //compute alpha and Chi coefficients: */
         if (q_exp==1){
@@ -264 +270 @@
         }
         Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
-
-        /*Check to prevent dividing by zero if vmag==0*/
-        if(vmag==0.) alpha2=0.;
-        else if (Neff==0) alpha2=0.0;
-        else alpha2= Neff * C_param * pow((Chi/(1 + alpha * pow(Chi,q_exp))),1/n) * 1/vmag;
+        Gamma=(Chi/(1. + alpha * pow(Chi,q_exp)));
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(vmag==0.) alpha2=0.; 
+        else    if (Neff==0) alpha2=0.0;
+        else    alpha2=Neff * C_param * pow(Gamma,1./n) * 1/vmag;
+
         _assert_(!xIsNan<IssmDouble>(alpha2));
         /*Assign output pointers:*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

@@ -48 +48 @@
                 switch(control){
                         /*List of supported controls*/
-                        case BalancethicknessThickeningRateEnum: 
+                  case BalancethicknessThickeningRateEnum: 
                         case VxEnum:
                         case VyEnum:
                         case ThicknessEnum:
-                        case FrictionCoefficientEnum:
-                        case BalancethicknessApparentMassbalanceEnum:
+                  case FrictionCoefficientEnum:
+                  case FrictionAsEnum:
+                  case BalancethicknessApparentMassbalanceEnum:
                         case BalancethicknessOmegaEnum:
                         case MaterialsRheologyBEnum: