Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 5250) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 5251) @@ -226,4 +226,5 @@ TemperatureOldEnum, ThicknessEnum, + ThicknessObsEnum, TypeEnum, VelEnum, @@ -346,10 +347,10 @@ /*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 Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 5250) +++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp (revision 5251) @@ -200,4 +200,5 @@ case TemperatureOldEnum : return "TemperatureOld"; case ThicknessEnum : return "Thickness"; + case ThicknessObsEnum : return "ThicknessObs"; case TypeEnum : return "Type"; case VelEnum : return "Vel"; Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 5250) +++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp (revision 5251) @@ -198,4 +198,5 @@ 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; Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 5250) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 5251) @@ -32,8 +32,14 @@ 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,3 +67,6 @@ xfree((void**)&iomodel->melting_rate); xfree((void**)&iomodel->accumulation_rate); + xfree((void**)&iomodel->thickness_obs); + xfree((void**)&iomodel->weights); + xfree((void**)&iomodel->control_parameter); } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5250) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5251) @@ -389,4 +389,7 @@ InputUpdateFromSolutionOneDof(solution,ThicknessEnum); } + else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){ + InputUpdateFromSolutionAdjointBalancedthickness( solution); + } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ InputUpdateFromSolutionOneDof(solution,VelEnum); @@ -657,9 +660,5 @@ /*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,5 +677,5 @@ 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,4 +735,7 @@ else ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + } + else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){ + CreatePVectorAdjointBalancedthickness(pg); } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ @@ -869,4 +871,7 @@ 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,4 +1210,87 @@ } /*}}}*/ +/*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; igGetParameterValue(&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;iDuplicateInput(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,4 +2259,8 @@ 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,4 +4330,94 @@ } /*}}}*/ +/*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; igGetParameterValue(&thickness, &gauss_l1l2l3[0]); + thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]); + weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]); + + for (i=0;iinputs->AddInput(new TriaVertexInput(AdjointEnum,lambda)); + + /*Free ressources:*/ + xfree((void**)&doflist); + +} +/*}}}*/ /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/ void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5250) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5251) @@ -82,4 +82,5 @@ 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,4 +131,5 @@ void CreatePVectorAdjointHoriz(Vec pg); void CreatePVectorAdjointStokes(Vec pg); + void CreatePVectorAdjointBalancedthickness(Vec pg); void CreatePVectorDiagnosticHutter(Vec pg); void CreatePVectorPrognostic_CG(Vec pg); @@ -149,4 +151,5 @@ 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); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 5250) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 5251) @@ -56,4 +56,5 @@ xfree((void**)&this->gridonstokes); xfree((void**)&this->borderstokes); + xfree((void**)&this->thickness_obs); xfree((void**)&this->thickness); xfree((void**)&this->surface); @@ -267,4 +268,5 @@ this->gridonstokes=NULL; this->borderstokes=NULL; + this->thickness_obs=NULL; this->thickness=NULL; this->surface=NULL; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 5250) +++ /issm/trunk/src/c/objects/IoModel.h (revision 5251) @@ -59,4 +59,5 @@ /*observations: */ + double* thickness_obs; double* vx_obs; double* vy_obs; Index: /issm/trunk/src/c/solutions/adjoint_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/adjoint_core.cpp (revision 5250) +++ /issm/trunk/src/c/solutions/adjoint_core.cpp (revision 5251) @@ -26,13 +26,27 @@ /*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); Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 5250) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 5251) @@ -28,14 +28,14 @@ /*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,10 +49,18 @@ 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,6 +68,16 @@ /*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.*/