Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18002) +++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18003) @@ -37,102 +37,72 @@ }/*}}}*/ ElementMatrix* AdjointBalancethickness2Analysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented"); Balancethickness2Analysis* analysis = new Balancethickness2Analysis(); ElementMatrix* Ke = analysis->CreateKMatrix(element); delete analysis; - /*Transpose and return Ke*/ - Ke->Transpose(); return Ke; }/*}}}*/ ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/ - _error_("not implemented"); - /*Intermediaries*/ - int domaintype; - Element* basalelement; - - /*Get basal element*/ - element->FindParam(&domaintype,DomainTypeEnum); - switch(domaintype){ - case Domain2DhorizontalEnum: - basalelement = element; - break; - case Domain3DEnum: - if(!element->IsOnBase()) return NULL; - basalelement = element->SpawnBasalElement(); - break; - default: _error_("mesh "<GetNumberOfNodes(); + int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector and vectors*/ - ElementVector* pe = basalelement->NewElementVector(SSAApproximationEnum); + ElementVector* pe = element->NewElementVector(SSAApproximationEnum); IssmDouble* basis = xNew(numnodes); IssmDouble* dbasis = xNew(2*numnodes); /*Retrieve all inputs and parameters*/ - basalelement->GetVerticesCoordinates(&xyz_list); - basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum); - basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum); - Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* thicknessobs_input = basalelement->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); - Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&num_responses,InversionNumCostFunctionsEnum); + element->FindParam(&responses,NULL,InversionCostFunctionsEnum); + Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* potential_input = element->GetInput(PotentialEnum); _assert_(potential_input); + Input* vx_input = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input); + Input* nux_input = element->GetInput(BalancethicknessNuxEnum); _assert_(nux_input); + Input* nuy_input = element->GetInput(BalancethicknessNuyEnum); _assert_(nuy_input); /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - basalelement->NodalFunctions(basis,gauss); - basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - thickness_input->GetInputValue(&thickness, gauss); - thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); - thicknessobs_input->GetInputValue(&thicknessobs, gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + nux_input->GetInputValue(&nux,gauss); + nuy_input->GetInputValue(&nuy,gauss); + potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss); + thicknessobs_input->GetInputValue(&hobs,gauss); + vbar2 = (nux*nux*vx*vx + nuy*nuy*vy*vy); + hu2 = hobs*hobs*vbar2; + + NUMx = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2); + NUMy = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2); + DEN = vbar2*vbar2; + /*Loop over all requested responses*/ for(int resp=0;respGetInputValue(&weight,gauss,responses[resp]); switch(responses[resp]){ - case ThicknessAbsMisfitEnum: - for(i=0;ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i]; + case Balancethickness2MisfitEnum: + for(i=0;ivalues[i]+=(NUMx+NUMy)/DEN *weight*Jdet*gauss->weight*basis[i]; break; - case ThicknessAbsGradientEnum: - for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight; - for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight; - break; - case ThicknessAlongGradientEnum: - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight; - break; - case ThicknessAcrossGradientEnum: - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel = sqrt(vx*vx+vy*vy); - vx = vx/(vel+1.e-9); - vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight; - break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); } @@ -144,7 +114,6 @@ xDelete(xyz_list); xDelete(basis); xDelete(dbasis); - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; delete gauss; return pe; }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18002) +++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18003) @@ -182,7 +182,7 @@ element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); - for(int i=0;ivalues[i] += - Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i]; + for(int i=0;ivalues[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i]; } /*Clean up and return*/ Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18002) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18003) @@ -11,8 +11,7 @@ void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){ /*Intermediary*/ - int i; - int counter; + int control; Element *element = NULL; Material *material = NULL; int num_control_type; @@ -32,8 +31,8 @@ iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum); - for(i=0;i(iomodel->Data(InversionControlParametersEnum)[i]); + for(int i=0;i(iomodel->Data(InversionControlParametersEnum)[i]); switch(control){ /*List of supported controls*/ case BalancethicknessThickeningRateEnum: @@ -43,6 +42,7 @@ case FrictionCoefficientEnum: case BalancethicknessNuxEnum: case BalancethicknessNuyEnum: + case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(1,control); break; @@ -55,8 +55,8 @@ } /*Update elements: */ - counter=0; - for (i=0;inumberofelements;i++){ + int counter=0; + for(int i=0;inumberofelements;i++){ if(iomodel->my_elements[i]){ element=(Element*)elements->GetObjectByOffset(counter); element->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements. @@ -65,5 +65,15 @@ } /*Free data: */ - iomodel->DeleteData(5+7,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionCostFunctionsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,ThicknessEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,DamageDEnum); + for(int i=0;i(iomodel->Data(InversionControlParametersEnum)[i]); + switch(control){ + case MaterialsRheologyBbarEnum: iomodel->DeleteData(1,MaterialsRheologyBEnum); break; + case DamageDbarEnum: iomodel->DeleteData(1,DamageDEnum); break; + default: iomodel->DeleteData(1,control); + } + } + iomodel->DeleteData(5,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionCostFunctionsEnum,InversionMinParametersEnum,InversionMaxParametersEnum); + + } Index: ../trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/control_core.cpp (revision 18002) +++ ../trunk-jpl/src/c/cores/control_core.cpp (revision 18003) @@ -144,18 +144,23 @@ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*set analysis type to compute velocity: */ - if (solution_type==SteadystateSolutionEnum || solution_type==StressbalanceSolutionEnum){ - femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + switch(solution_type){ + case SteadystateSolutionEnum: + case StressbalanceSolutionEnum: + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + break; + case BalancethicknessSolutionEnum: + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + break; + case BalancethicknessSoftSolutionEnum: + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + break; + case Balancethickness2SolutionEnum: + femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum); + break; + default: + _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet"); } - else if (solution_type==BalancethicknessSolutionEnum){ - femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); - } - else if (solution_type==BalancethicknessSoftSolutionEnum){ - femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); - } - else{ - _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet"); - } /*update parameter according to scalar: */ //false means: do not save control InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false); @@ -170,6 +175,9 @@ else if (solution_type==BalancethicknessSolutionEnum){ solutionsequence_linear(femmodel); } + else if (solution_type==Balancethickness2SolutionEnum){ + solutionsequence_linear(femmodel); + } else if (solution_type==BalancethicknessSoftSolutionEnum){ /*Don't do anything*/ } Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 18002) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 18003) @@ -435,17 +435,17 @@ case IceVolumeEnum: this->IceVolumex(responses); break; case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break; case MinVelEnum: this->MinVelx(responses); break; - case MaxVelEnum: this->MaxVelx( responses); break; + case MaxVelEnum: this->MaxVelx(responses); break; case MinVxEnum: this->MinVxx(responses); break; - case MaxVxEnum: this->MaxVxx( responses); break; - case MaxAbsVxEnum: this->MaxAbsVxx( responses); break; + case MaxVxEnum: this->MaxVxx(responses); break; + case MaxAbsVxEnum: this->MaxAbsVxx(responses); break; case MinVyEnum: this->MinVyx(responses); break; - case MaxVyEnum: this->MaxVyx( responses); break; - case MaxAbsVyEnum: this->MaxAbsVyx( responses); break; + case MaxVyEnum: this->MaxVyx(responses); break; + case MaxAbsVyEnum: this->MaxAbsVyx(responses); break; case MinVzEnum: this->MinVzx(responses); break; - case MaxVzEnum: this->MaxVzx( responses); break; - case MaxAbsVzEnum: this->MaxAbsVzx( responses); break; - case MassFluxEnum: this->MassFluxx( responses); break; + case MaxVzEnum: this->MaxVzx(responses); break; + case MaxAbsVzEnum: this->MaxAbsVzx(responses); break; + case MassFluxEnum: this->MassFluxx(responses); break; case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break; case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break; case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials,parameters); break; @@ -458,6 +458,7 @@ case RheologyBbarAbsGradientEnum: RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break; case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break; case BalancethicknessMisfitEnum: BalancethicknessMisfitx(responses); break; + case Balancethickness2MisfitEnum: Balancethickness2Misfitx(responses); break; case TotalSmbEnum: this->TotalSmbx(responses); break; case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break; case VelEnum: this->ElementResponsex(responses,VelEnum); break; @@ -1199,6 +1200,61 @@ *presponse=J; }/*}}}*/ +void FemModel::Balancethickness2Misfitx(IssmDouble* presponse){/*{{{*/ + + /*output: */ + IssmDouble J=0.; + IssmDouble J_sum; + + IssmDouble weight,thicknessobs,thickness; + IssmDouble Jdet; + IssmDouble* xyz_list = NULL; + + /*Compute Misfit: */ + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + + /*If on water, return 0: */ + if(!element->IsIceInElement()) continue; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input =element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum); + thickness_input->GetInputValue(&thickness,gauss); + thicknessobs_input->GetInputValue(&thicknessobs,gauss); + gauss->GaussPoint(ig); + + J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight; + } + + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; + } + + /*Sum all J from all cpus of the cluster:*/ + ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + J=J_sum; + + /*Assign output pointers: */ + *presponse=J; + +}/*}}}*/ void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/ /*output: */ Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18002) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18003) @@ -990,6 +990,7 @@ switch(control){ /*yts conversion*/ case BalancethicknessThickeningRateEnum: + case BalancethicknessApparentMassbalanceEnum: case VxEnum: case VyEnum: if(iomodel->Data(control)){ @@ -2802,8 +2803,11 @@ case ThicknessEnum: GradjThicknessBalancethicknessSoft(gradient,control_index); break; + case BalancethicknessApparentMassbalanceEnum: + GradjAdotBulancethickness2(gradient,control_index); + break; default: - _error_("control type not supported yet: " << control_type); + _error_("control type not supported yet: " << EnumToStringx(control_type)); } /*Now deal with ∂J/∂alpha*/ @@ -2820,6 +2824,7 @@ case ThicknessAlongGradientEnum: case ThicknessAcrossGradientEnum: case BalancethicknessMisfitEnum: + case Balancethickness2MisfitEnum: case SurfaceAbsVelMisfitEnum: case SurfaceRelVelMisfitEnum: case SurfaceLogVelMisfitEnum: @@ -3167,6 +3172,21 @@ gradient->SetValues(NUMVERTICES,vertexpidlist,gradient_g,INS_VAL); } /*}}}*/ +void Tria::GradjAdotBulancethickness2(Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + int vertexpidlist[NUMVERTICES]; + IssmDouble lambda[NUMVERTICES]; + IssmDouble gradient_g[NUMVERTICES]; + + /*Compute Gradient*/ + GradientIndexing(&vertexpidlist[0],control_index); + GetInputListOnVertices(&lambda[0],AdjointEnum); + for(int i=0;iSetValues(NUMVERTICES,vertexpidlist,gradient_g,INS_VAL); +} +/*}}}*/ void Tria::GradjVxBalancedthickness(Vector* gradient,int control_index){/*{{{*/ /*Intermediaries*/ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18002) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18003) @@ -148,6 +148,7 @@ void GradjVxBalancedthickness(Vector* gradient,int control_index); void GradjVyBalancedthickness(Vector* gradient,int control_index); void GradjThicknessBalancethicknessSoft(Vector* gradient,int control_index); + void GradjAdotBulancethickness2(Vector* gradient,int control_index); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void ControlInputGetGradient(Vector* gradient,int enum_type,int control_index); Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 18002) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 18003) @@ -76,6 +76,7 @@ void IceVolumeAboveFloatationx(IssmDouble* pV); void ElementResponsex(IssmDouble* presponse,int response_enum); void BalancethicknessMisfitx(IssmDouble* pV); + void Balancethickness2Misfitx(IssmDouble* pV); #ifdef _HAVE_DAKOTA_ void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses); #endif