Changeset 17971

Timestamp:

05/09/14 11:30:18 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: inversions for FS blowbands

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (4 diffs)
• classes/Elements/Seg.cpp (modified) (1 diff)
• classes/Elements/Seg.h (modified) (2 diffs)
• classes/Elements/Tria.cpp (modified) (3 diffs)
• cores/adjointstressbalance_core.cpp (modified) (1 diff)

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

@@ -332 +332 @@
         /*Prepare coordinate system list*/
         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
 
@@ -513 +514 @@
         delete gauss;
         return pe;
-        _error_("S");
 }/*}}}*/
 ElementVector* AdjointHorizAnalysis::CreatePVectorHO(Element* element){/*{{{*/
@@ -951 +951 @@
 }/*}}}*/
 void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
-        int          i;
+        int          i,dim;
         int*         vdoflist=NULL;
         int*         pdoflist=NULL;
         IssmDouble   FSreconditioning;
 
+        element->FindParam(&dim,DomainDimensionEnum);
+
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->NumberofNodesVelocity();
         int pnumnodes = element->NumberofNodesPressure();
-        int vnumdof   = vnumnodes*3;
+        int vnumdof   = vnumnodes*dim;
         int pnumdof   = pnumnodes*1;
 
@@ -970 +972 @@
 
         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-        for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
 

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

@@ -116 +116 @@
 
 }/*}}}*/
+void       Seg::GradientIndexing(int* indexing,int control_index){/*{{{*/
+
+        /*Get some parameters*/
+        int num_controls;
+        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+
+        /*get gradient indices*/
+        for(int i=0;i<NUMVERTICES;i++){
+                indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
+        }
+
+}
+/*}}}*/
+void       Seg::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        int        i;
+        int        analysis_type;
+        int        vertexpidlist[NUMVERTICES];
+        int        connectivity[NUMVERTICES];
+        IssmDouble vx,lambda,alpha_complement,drag,Jdet;
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble dk[NDOF2]; 
+        IssmDouble grade_g[NUMVERTICES]={0.0};
+        IssmDouble grade_g_gaussian[NUMVERTICES];
+        IssmDouble basis[3];
+        Friction*  friction=NULL;
+        GaussSeg  *gauss=NULL;
+
+        if(IsFloating())return;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        GradientIndexing(&vertexpidlist[0],control_index);
+        this->GetVerticesConnectivityList(&connectivity[0]);
+
+        /*Build frictoin element, needed later: */
+        friction=new Friction(this,1);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                   _assert_(adjointx_input);
+        Input* vx_input=inputs->GetInput(VxEnum);                               _assert_(vx_input);
+        Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussSeg(4);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctions(basis, gauss);
+
+                /*Build alpha_complement_list: */
+                friction->GetAlphaComplement(&alpha_complement,gauss);
+
+                dragcoefficient_input->GetInputValue(&drag, gauss);
+                adjointx_input->GetInputValue(&lambda, gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for (i=0;i<NUMVERTICES;i++){
+                        grade_g_gaussian[i]=-2*drag*alpha_complement*(lambda*vx)*Jdet*gauss->weight*basis[i];
+                }
+
+                /*Add gradje_g_gaussian vector to gradje_g: */
+                for(i=0;i<NUMVERTICES;i++){
+                        _assert_(!xIsNan<IssmDouble>(grade_g[i]));
+                        grade_g[i]+=grade_g_gaussian[i];
+                }
+        }
+        gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+}
+/*}}}*/
 bool       Seg::IsIcefront(void){/*{{{*/
 

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

@@ -169 +169 @@
 
                 IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};
-                void       GradientIndexing(int* indexing,int control_index){_error_("not implemented yet");};
+                void       GradientIndexing(int* indexing,int control_index);
                 void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
                 void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
@@ -176 +176 @@
                 void       GradjDSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
                 void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
-                void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
+                void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
                 void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
                 void       GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};

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

@@ -2863 +2863 @@
         /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
 
+        int   approximation;
+        Seg*  seg=NULL;
+
         /*If on water, grad = 0: */
         if(!IsIceInElement()) return;
@@ -2869 +2872 @@
         switch(control_type){
                 case FrictionCoefficientEnum:
+                        inputs->GetInputValue(&approximation,ApproximationEnum);
+                        switch(approximation){
+                                case SSAApproximationEnum:
+                                        GradjDragSSA(gradient,control_index);
+                                        break;
+                                case FSApproximationEnum:
+                                        GradjDragFS(gradient,control_index);
+                                        break;
+                                case NoneApproximationEnum:
+                                        /*Gradient is 0*/
+                                        break;
+                                default:
+                                        _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
+                        }
+                        break;
                         GradjDragSSA(gradient,control_index);
                         break;
@@ -3163 +3181 @@
         delete gauss;
         delete friction;
+}
+/*}}}*/
+void       Tria::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+        /*Gradient is 0 if on shelf or not on bed*/
+        if(IsFloating() || !IsOnBase()) return;
+
+        int index1,index2;
+        this->EdgeOnBaseIndices(&index1,&index2);
+        Seg* seg = SpawnSeg(index1,index2);
+        seg->GradjDragFS(gradient,control_index);
+        seg->DeleteMaterials(); delete seg;
 }
 /*}}}*/

issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp

@@ -38 +38 @@
                 if(VerboseSolution()) _printf0_("   saving results\n");
                 if (isFS){
-                        int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum};
-                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
+                        //int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum};
+                        //femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
+                        int outputs[3] = {AdjointxEnum,AdjointyEnum,AdjointpEnum};
+                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
                 }
                 else{