Changeset 18771

Timestamp:

11/10/14 15:01:21 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added L1L2 inversions (very similar to SSA)

Location:

issm/trunk-jpl/src

Files:

• c/analyses/AdjointHorizAnalysis.cpp (modified) (8 diffs)
• c/analyses/AdjointHorizAnalysis.h (modified) (2 diffs)
• m/classes/inversion.m (modified) (1 diff)
• m/classes/inversion.py (modified) (1 diff)

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

@@ -42 +42 @@
                 case SSAApproximationEnum: 
                         return CreateKMatrixSSA(element);
+                case L1L2ApproximationEnum: 
+                        return CreateKMatrixL1L2(element);
                 case HOApproximationEnum: 
                         return CreateKMatrixHO(element);
@@ -139 +141 @@
         return Ke;
 }/*}}}*/
+ElementMatrix* AdjointHorizAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        bool incomplete_adjoint;
+
+        /*Initialize Jacobian with regular L1L2 (first part of the Gateau derivative)*/
+        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
+        ElementMatrix* Ke=analysis->CreateKMatrix(element);
+        delete analysis;
+
+        /*return*/
+        element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
+        if(!incomplete_adjoint){
+                _error_("Exact adjoint not supported yet for L1L2 model");
+        }
+        return Ke;
+}/*}}}*/
 ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
 
@@ -300 +319 @@
                 case SSAApproximationEnum: 
                         return CreatePVectorSSA(element);
+                case L1L2ApproximationEnum: 
+                        return CreatePVectorL1L2(element);
                 case HOApproximationEnum: 
                         return CreatePVectorHO(element);
@@ -897 +918 @@
         return pe;
 }/*}}}*/
+ElementVector* AdjointHorizAnalysis::CreatePVectorL1L2(Element* element){/*{{{*/
+
+        /*Same as SSA*/
+        return this->CreatePVectorSSA(element);
+}/*}}}*/
 void AdjointHorizAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");
@@ -946 +972 @@
                         switch(approximation){
                                 case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
+                                case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_index); break;
                                 case HOApproximationEnum:  GradientJDragHO( element,gradient,control_index); break;
                                 case FSApproximationEnum:  GradientJDragFS( element,gradient,control_index); break;
@@ -955 +982 @@
                         switch(approximation){
                                 case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_index); break;
+                                case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_index); break;
                                 case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_index); break;
                                 case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_index); break;
@@ -1223 +1251 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+void AdjointHorizAnalysis::GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Same as SSA*/
+        return this->GradientJDragSSA(element,gradient,control_index);
+}/*}}}*/
 void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
@@ -1446 +1479 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+void AdjointHorizAnalysis::GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+           /*Same as SSA*/
+           return this->GradientJBbarSSA(element,gradient,control_index);
+}/*}}}*/
 void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 

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

@@ -26 +26 @@
                 ElementMatrix* CreateKMatrix(Element* element);
                 ElementMatrix* CreateKMatrixSSA(Element* element);
+                ElementMatrix* CreateKMatrixL1L2(Element* element);
                 ElementMatrix* CreateKMatrixHO(Element* element);
                 ElementMatrix* CreateKMatrixFS(Element* element);
                 ElementVector* CreatePVector(Element* element);
                 ElementVector* CreatePVectorSSA(Element* element);
+                ElementVector* CreatePVectorL1L2(Element* element);
                 ElementVector* CreatePVectorHO(Element* element);
                 ElementVector* CreatePVectorFS(Element* element);
@@ -37 +39 @@
                 void GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
+                void GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
                 void GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index);

issm/trunk-jpl/src/m/classes/inversion.m

@@ -146 +146 @@
                         %Only SSA, HO and FS are supported right now
                         if solution==StressbalanceSolutionEnum()
-                                if ~(md.flowequation.isSSA || md.flowequation.isHO || md.flowequation.isFS),
+                                if ~(md.flowequation.isSSA || md.flowequation.isHO || md.flowequation.isFS || md.flowequation.isL1L2),
                                         md = checkmessage(md,['inversion can only be performed for SSA, HO or FS ice flow models']);
                                 end

issm/trunk-jpl/src/m/classes/inversion.py

@@ -131 +131 @@
                 #Only SSA, HO and FS are supported right now
                 if solution==StressbalanceSolutionEnum():
-                        if not (md.flowequation.isSSA or md.flowequation.isHO or md.flowequation.isFS):
+                        if not (md.flowequation.isSSA or md.flowequation.isHO or md.flowequation.isFS or md.flowequation.isL1L2):
                                 md.checkmessage("'inversion can only be performed for SSA, HO or FS ice flow models");