Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18003) @@ -38,71 +38,61 @@ 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*/ @@ -111,26 +101,6 @@ switch(responses[resp]){ - case ThicknessAbsMisfitEnum: - for(i=0;ivalues[i]+=(thicknessobs-thickness)*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; + case Balancethickness2MisfitEnum: + for(i=0;ivalues[i]+=(NUMx+NUMy)/DEN *weight*Jdet*gauss->weight*basis[i]; break; default: @@ -145,5 +115,4 @@ xDelete(basis); xDelete(dbasis); - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; delete gauss; return pe; Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18003) @@ -183,5 +183,5 @@ 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]; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18003) @@ -991,4 +991,5 @@ /*yts conversion*/ case BalancethicknessThickeningRateEnum: + case BalancethicknessApparentMassbalanceEnum: case VxEnum: case VyEnum: @@ -2803,6 +2804,9 @@ 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)); } @@ -2821,4 +2825,5 @@ case ThicknessAcrossGradientEnum: case BalancethicknessMisfitEnum: + case Balancethickness2MisfitEnum: case SurfaceAbsVelMisfitEnum: case SurfaceRelVelMisfitEnum: @@ -3164,4 +3169,19 @@ GetInputListOnVertices(&lambda[0],AdjointEnum); for(int i=0;iSetValues(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); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18002) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18003) @@ -149,4 +149,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18003) @@ -436,15 +436,15 @@ 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; @@ -459,4 +459,5 @@ 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; @@ -1200,4 +1201,59 @@ }/*}}}*/ +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){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18002) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18003) @@ -77,4 +77,5 @@ 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); Index: /issm/trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 18003) @@ -145,15 +145,20 @@ /*set analysis type to compute velocity: */ - if (solution_type==SteadystateSolutionEnum || solution_type==StressbalanceSolutionEnum){ - femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); - } - 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"); + 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"); } @@ -171,4 +176,7 @@ solutionsequence_linear(femmodel); } + else if (solution_type==Balancethickness2SolutionEnum){ + solutionsequence_linear(femmodel); + } else if (solution_type==BalancethicknessSoftSolutionEnum){ /*Don't do anything*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18002) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18003) @@ -12,6 +12,5 @@ /*Intermediary*/ - int i; - int counter; + int control; Element *element = NULL; Material *material = NULL; @@ -33,6 +32,6 @@ 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*/ @@ -44,4 +43,5 @@ case BalancethicknessNuxEnum: case BalancethicknessNuyEnum: + case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(1,control); break; @@ -56,6 +56,6 @@ /*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); @@ -66,4 +66,14 @@ /*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); + + }