Changeset 5251

Timestamp:

08/13/10 15:25:10 (15 years ago)

Author:

Mathieu Morlighem

Message:

prepared CM balancedthickness (To be continued)

Location:

issm/trunk/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• EnumDefinitions/EnumToString.cpp (modified) (1 diff)
• EnumDefinitions/StringToEnum.cpp (modified) (1 diff)
• modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (modified) (2 diffs)
• objects/Elements/Tria.cpp (modified) (10 diffs)
• objects/Elements/Tria.h (modified) (3 diffs)
• objects/IoModel.cpp (modified) (2 diffs)
• objects/IoModel.h (modified) (1 diff)
• solutions/adjoint_core.cpp (modified) (1 diff)
• solutions/objectivefunctionC.cpp (modified) (3 diffs)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -226 +226 @@
         TemperatureOldEnum,
         ThicknessEnum,
+        ThicknessObsEnum,
         TypeEnum,
         VelEnum,
@@ -346 +347 @@
 
 /*Functions on enums: */
-int EnumIsElement(int en);
-int EnumIsLoad(int en);
-int EnumIsMaterial(int en);
-char* EnumToString(int enum_type);
-int StringToEnum(char* string);
-char* EnumToModelField(int en);
+int   EnumIsElement(int en);
+int   EnumIsLoad(int en);
+int   EnumIsMaterial(int en);
+char *EnumToString(int enum_type);
+int   StringToEnum(char *string);
+char *EnumToModelField(int  en);
 
 #endif

issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

@@ -200 +200 @@
                 case TemperatureOldEnum : return "TemperatureOld";
                 case ThicknessEnum : return "Thickness";
+                case ThicknessObsEnum : return "ThicknessObs";
                 case TypeEnum : return "Type";
                 case VelEnum : return "Vel";

issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

@@ -198 +198 @@
         else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
         else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
+        else if (strcmp(name,"ThicknessObs")==0) return ThicknessObsEnum;
         else if (strcmp(name,"Type")==0) return TypeEnum;
         else if (strcmp(name,"Vel")==0) return VelEnum;

issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp

@@ -32 +32 @@
         IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
         IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
-
         if (iomodel->dim==3){
                 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
                 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
+        }
+        if(iomodel->control_analysis){
+                IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs");
+                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
+                if(strcmp(iomodel->control_type,"dhdt")==0){
+                        IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
+                }
         }
 
@@ -61 +67 @@
         xfree((void**)&iomodel->melting_rate);
         xfree((void**)&iomodel->accumulation_rate);
+        xfree((void**)&iomodel->thickness_obs);
+        xfree((void**)&iomodel->weights);
+        xfree((void**)&iomodel->control_parameter);
 }

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -389 +389 @@
                 InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
         }
+        else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){
+                InputUpdateFromSolutionAdjointBalancedthickness( solution);
+        }
         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
                 InputUpdateFromSolutionOneDof(solution,VelEnum);
@@ -657 +660 @@
 
         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==DiagnosticHorizAnalysisEnum){
-                CreateKMatrixDiagnosticMacAyeal( Kgg);
-        }
-        else if (analysis_type==AdjointHorizAnalysisEnum){
-                /*Same as diagnostic*/
+        if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){
                 CreateKMatrixDiagnosticMacAyeal( Kgg);
         }
@@ -678 +677 @@
                  ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType()));
         }
-        else if (analysis_type==BalancedthicknessAnalysisEnum){
+        else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){
                 if (GetElementType()==P1Enum)
                  CreateKMatrixBalancedthickness_CG( Kgg);
@@ -736 +735 @@
                 else
                  ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType()));
+        }
+        else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){
+                CreatePVectorAdjointBalancedthickness(pg);
         }
         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
@@ -869 +871 @@
         else if (control_type==RheologyB2dEnum){
                 GradjB(gradient);
+        }
+        else if (control_type==DhDtEnum){
+                GradjDhDt(gradient);
         }
         else ISSMERROR("%s%i","control type not supported yet: ",control_type);
@@ -1205 +1210 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GradjDhDt{{{1*/
+void  Tria::GradjDhDt(Vec gradient){
+
+        int i;
+
+        /* node data: */
+        const int    numgrids=3;
+        const int    NDOF1=1;
+        const int    numdof=NDOF1*numgrids;
+        double       xyz_list[numgrids][3];
+        int          doflist1[numgrids];
+
+        /* grid data: */
+        double DhDt[numgrids];
+
+        /* gaussian points: */
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /*element vector at the gaussian points: */
+        double  grade_g[numgrids]={0.0};
+        double  grade_g_gaussian[numgrids];
+
+        /* Jacobian: */
+        double Jdet;
+
+        /*nodal functions: */
+        double l1l2l3[3];
+        double  lambda;
+
+        /*Inputs*/
+        Input* adjoint_input=NULL;
+
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+        GetDofList1(&doflist1[0]);
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
+
+        /*Retrieve all inputs we will be needing: */
+        adjoint_input=inputs->GetInput(AdjointEnum);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /*Get thickness: */
+                adjoint_input->GetParameterValue(&lambda, gauss_l1l2l3);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+                /* Get nodal functions value at gaussian point:*/
+                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+                /*DhDtuild gradje_g_gaussian vector (actually -dJ/dDhDt): */
+                for (i=0;i<numgrids;i++){
+                        //standard gradient dJ/dki
+                        grade_g_gaussian[i]=lambda*Jdet*gauss_weight*l1l2l3[i]; 
+                }
+
+                /*Add grade_g_gaussian to grade_g: */
+                for( i=0; i<numgrids;i++) grade_g[i]+=grade_g_gaussian[i];
+        }
+
+        /*Add grade_g to global vector gradient: */
+        VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
+
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+}
+/*}}}*/
 /*FUNCTION Tria::InputControlUpdate{{{1*/
 void  Tria::InputControlUpdate(double scalar,bool save_parameter){
@@ -1253 +1341 @@
                 /*Finally: save input if requested*/
                 if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
+
+        }
+        else if(control_type==DhDtEnum){
+
+                /*First, get revert to previous parameter value (kept in ControlParameter input)*/
+                this->inputs->DuplicateInput(ControlParameterEnum,DhDtEnum);
+
+                /*Now, update using the gradient new = old + scalar * gradient*/
+                this->inputs->AXPY(DhDtEnum,scalar,GradientEnum);
+
+                /*Constrain input*/
+                input=(Input*)this->inputs->GetInput(DhDtEnum); ISSMASSERT(input);
+                input->Constrain(cm_min,cm_max);
+
+                /*Finally: save input if requested*/
+                if (save_parameter) inputs->DuplicateInput(DhDtEnum,ControlParameterEnum);
 
         }
@@ -2155 +2259 @@
                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->control_parameter[tria_vertex_ids[i]-1];
                 this->inputs->AddInput(new TriaVertexInput(ControlParameterEnum,nodeinputs));
+        }
+        if (iomodel->thickness_obs) {
+                for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_obs[tria_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new TriaVertexInput(ThicknessObsEnum,nodeinputs));
         }
         if (iomodel->melting_rate) {
@@ -4222 +4330 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreatePVectorAdjointBalancedthickness{{{1*/
+void Tria::CreatePVectorAdjointBalancedthickness(Vec p_g){
+
+        int i;
+
+        /* node data: */
+        const int    numgrids=3;
+        const int    numdof=1*numgrids;
+        double       xyz_list[numgrids][3];
+        int*         doflist=NULL;
+
+        /* gaussian points: */
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /* parameters: */
+        double  thickness,thicknessobs,weight;
+
+        /*element vector : */
+        double  pe_g[numdof]={0.0};
+        double  pe_g_gaussian[numdof];
+
+        /* Jacobian: */
+        double Jdet;
+
+        /*nodal functions: */
+        double l1l2l3[3];
+
+        /*inputs: */
+        Input* thickness_input=NULL;
+        Input* thicknessobs_input=NULL;
+        Input* weights_input=NULL;
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist);
+
+        /*Retrieve all inputs we will be needing: */
+        thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
+        thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
+        weights_input     =inputs->GetInput(ThicknessObsEnum);ISSMASSERT(weights_input);
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+
+                /* Get nodal functions value at gaussian point:*/
+                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
+
+                /*Get parameters at gauss point*/
+                thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
+                thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]);
+                weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]);
+
+                for (i=0;i<numgrids;i++){
+                        pe_g_gaussian[i]=(thicknessobs-thickness)*weight*Jdet*gauss_weight*l1l2l3[i]; 
+                }
+
+                /*Add pe_g_gaussian vector to pe_g: */
+                for( i=0; i<numdof; i++){
+                        pe_g[i]+=pe_g_gaussian[i];
+                }
+        }
+
+        /*Add pe_g to global vector p_g: */
+        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+        /*Free ressources:*/
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+        xfree((void**)&doflist);
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
 void Tria::CreatePVectorAdjointHoriz(Vec p_g){
@@ -5833 +6031 @@
 }
 /*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancedthickness {{{1*/
+void  Tria::InputUpdateFromSolutionAdjointBalancedthickness(double* solution){
+
+        int i;
+
+        const int    numvertices=3;
+        const int    numdofpervertex=1;
+        const int    numdof=numdofpervertex*numvertices;
+        int*         doflist=NULL;
+        double       values[numdof];
+        double       lambda[numvertices];
+
+        /*Get dof list: */
+        GetDofList(&doflist);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof;i++){
+                values[i]=solution[doflist[i]];
+        }
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numvertices;i++){
+                lambda[i]=values[i*numdofpervertex+0];
+        }
+
+        /*Add vx and vy as inputs to the tria element: */
+        this->inputs->AddInput(new TriaVertexInput(AdjointEnum,lambda));
+
+        /*Free ressources:*/
+        xfree((void**)&doflist);
+
+}
+/*}}}*/
 /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
 void  Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){

issm/trunk/src/c/objects/Elements/Tria.h

@@ -82 +82 @@
                 void   GradjB(Vec gradient);
                 void   GradjDrag(Vec gradient);
+                void   GradjDhDt(Vec gradient);
                 void   InputControlUpdate(double scalar,bool save_parameter);
                 bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
@@ -130 +131 @@
                 void      CreatePVectorAdjointHoriz(Vec pg);
                 void      CreatePVectorAdjointStokes(Vec pg);
+                void      CreatePVectorAdjointBalancedthickness(Vec pg);
                 void      CreatePVectorDiagnosticHutter(Vec pg);
                 void      CreatePVectorPrognostic_CG(Vec pg);
@@ -149 +151 @@
                 void    GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
                 void      GradjDragStokes(Vec gradient);
+                void      InputUpdateFromSolutionAdjointBalancedthickness( double* solution);
                 void      InputUpdateFromSolutionAdjointHoriz( double* solution);
                 void      InputUpdateFromSolutionDiagnosticHoriz( double* solution);

issm/trunk/src/c/objects/IoModel.cpp

@@ -56 +56 @@
         xfree((void**)&this->gridonstokes);
         xfree((void**)&this->borderstokes);
+        xfree((void**)&this->thickness_obs);
         xfree((void**)&this->thickness);
         xfree((void**)&this->surface);
@@ -267 +268 @@
         this->gridonstokes=NULL;
         this->borderstokes=NULL;
+        this->thickness_obs=NULL;
         this->thickness=NULL;
         this->surface=NULL;

issm/trunk/src/c/objects/IoModel.h

@@ -59 +59 @@
 
                 /*observations: */
+                double*  thickness_obs;
                 double*  vx_obs;
                 double*  vy_obs;

issm/trunk/src/c/solutions/adjoint_core.cpp

@@ -26 +26 @@
 
         /*set analysis type to compute velocity: */
-        _printf_("%s\n","      computing velocities");
-        if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
-        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-        solver_diagnostic_nonlinear(femmodel,conserve_loads); 
+        if (solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum){
+                _printf_("%s\n","      computing velocities");
+                if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
+                else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
+                solver_diagnostic_nonlinear(femmodel,conserve_loads); 
+        }
+        else if (solution_type==BalancedthicknessSolutionEnum){
+                femmodel->SetCurrentConfiguration(BalancedthicknessAnalysisEnum);
+                solver_linear(femmodel); 
+        }
+        else{
+                ISSMERROR("Solution %s not implemented yet",EnumToString(solution_type));
+        }
 
         /*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
         _printf_("%s\n","      computing adjoint");
-        if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
-        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
+        if (solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum){
+                if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
+                else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
+        }
+        else if (solution_type==BalancedthicknessSolutionEnum){
+                femmodel->SetCurrentConfiguration(BalancedthicknessAnalysisEnum,AdjointBalancedthicknessAnalysisEnum);
+        }
         solver_adjoint_linear(femmodel);
         

issm/trunk/src/c/solutions/objectivefunctionC.cpp

@@ -28 +28 @@
         
         /*parameters: */
-        FemModel* femmodel=NULL;
-        int     n;
-        double* optscal=NULL;
-        double* fit=NULL;
-        int     control_type;
-        int     analysis_type;
-        bool    control_steady;
-        bool  isstokes=false;
-        bool    conserve_loads=true;
-        bool    process_units=false;
+        FemModel *femmodel       = NULL;
+        int       n;
+        double   *optscal        = NULL;
+        double   *fit            = NULL;
+        int       control_type;
+        int       solution_type;
+        int       analysis_type;
+        bool      isstokes       = false;
+        bool      conserve_loads = true;
+        bool      process_units  = false;
         
         /*Recover finite element model: */
@@ -49 +49 @@
         femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
         femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
-        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
         femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 
         /*set analysis type to compute velocity: */
-        if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
-        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
+        if (solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum){
+                if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
+                else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
+        }
+        else if (solution_type==BalancedthicknessSolutionEnum){
+                femmodel->SetCurrentConfiguration(BalancedthicknessAnalysisEnum);
+        }
+        else{
+                ISSMERROR("Solution %s not implemented yet",EnumToString(solution_type));
+        }
 
         /*update parameter according to scalar: */ //false means: do not copy updated parameter onto ControlParameter input
@@ -60 +68 @@
 
         /*Run diagnostic with updated inputs: */
-        if(!control_steady) solver_diagnostic_nonlinear(femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
-        else                diagnostic_core(femmodel);  //We need a 3D velocity!! (vz is required for the next thermal run)
+        if (solution_type==SteadystateSolutionEnum){
+                diagnostic_core(femmodel);      //We need a 3D velocity!! (vz is required for the next thermal run)
+        }
+        else if (solution_type==DiagnosticSolutionEnum){
+                solver_diagnostic_nonlinear(femmodel,conserve_loads); 
+        }
+        else if (solution_type==BalancedthicknessSolutionEnum){
+                solver_linear(femmodel); 
+        }
+        else{
+                ISSMERROR("Solution %s not implemented yet",EnumToString(solution_type));
+        }
 
         /*Compute misfit for this velocity field.*/