Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23867) @@ -553,4 +553,7 @@ /*Nothing in P vector*/ break; + case RheologyBInitialguessMisfitEnum: + /*Nothing in P vector*/ + break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); @@ -788,4 +791,7 @@ /*Nothing in P vector*/ break; + case RheologyBInitialguessMisfitEnum: + /*Nothing in P vector*/ + break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); @@ -1036,4 +1042,7 @@ /*Nothing in P vector*/ break; + case RheologyBInitialguessMisfitEnum: + /*Nothing in P vector*/ + break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); @@ -1097,4 +1106,5 @@ case RheologyBbarAbsGradientEnum: GradientJBbarGradient(element,gradient,control_index); break; case RheologyBAbsGradientEnum: GradientJBGradient(element,gradient,control_index); break; + case RheologyBInitialguessMisfitEnum: GradientJBinitial(element,gradient,control_index); break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); } @@ -1398,4 +1408,68 @@ delete gauss; }/*}}}*/ +void AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + int domaintype; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + break; + case Domain2DverticalEnum: + break; + case Domain3DEnum: + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* rheology_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheology_input); + Input* rheology0_input = element->GetInput(RheologyBInitialguessEnum); _assert_(rheology0_input); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(basis,gauss); + weights_input->GetInputValue(&weight,gauss,RheologyBInitialguessMisfitEnum); + + /*Build alpha_complement_list: */ + rheology_input->GetInputValue(&B,gauss); + rheology0_input->GetInputValue(&B0,gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for(int i=0;iweight*basis[i]*(B-B0); + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; +}/*}}}*/ void AdjointHorizAnalysis::GradientJBHO(Element* element,Vector* gradient,int control_index){/*{{{*/ /*Intermediaries*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 23866) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 23867) @@ -38,4 +38,5 @@ void GradientJBbarFS(Element* element,Vector* gradient,int control_index); void GradientJBbarGradient(Element* element,Vector* gradient,int control_index); + void GradientJBinitial(Element* element,Vector* gradient,int control_index); void GradientJBbarL1L2(Element* element,Vector* gradient,int control_index); void GradientJBbarHO(Element* element,Vector* gradient,int control_index); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23867) @@ -2223,4 +2223,5 @@ case RheologyBbarAbsGradientEnum: RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; case RheologyBAbsGradientEnum: RheologyBAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case RheologyBInitialguessMisfitEnum: RheologyBInitialguessMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; case BalancethicknessMisfitEnum: BalancethicknessMisfitx(&double_result); break; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23867) @@ -60,4 +60,5 @@ if( cost_function==ThicknessAbsMisfitEnum) iomodel->FetchDataToInput(elements,"md.inversion.thickness_obs",InversionThicknessObsEnum); else if(cost_function==SurfaceAbsMisfitEnum) iomodel->FetchDataToInput(elements,"md.inversion.surface_obs",InversionSurfaceObsEnum); + else if(cost_function==RheologyBInitialguessMisfitEnum) iomodel->FetchDataToInput(elements,"md.materials.rheology_B",RheologyBInitialguessEnum); else if(cost_function==SurfaceAbsVelMisfitEnum || cost_function==SurfaceRelVelMisfitEnum Index: /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.cpp (revision 23867) @@ -78,2 +78,75 @@ return Jelem; } + +void RheologyBInitialguessMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ + + /*output: */ + IssmDouble J=0.; + IssmDouble J_sum; + + /*Compute Misfit: */ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + J+=RheologyBInitialguessMisfit(element); + } + + /*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: */ + *pJ=J; +} + +IssmDouble RheologyBInitialguessMisfit(Element* element){ + + int domaintype,numcomponents; + IssmDouble Jelem=0.; + IssmDouble misfit,Jdet; + IssmDouble B,B0,weight; + IssmDouble* xyz_list = NULL; + + /*If on water, return 0: */ + if(!element->IsIceInElement()) return 0.; + + /*Get problem dimension*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DverticalEnum: numcomponents = 1; break; + case Domain3DEnum: numcomponents = 2; break; + case Domain2DhorizontalEnum: numcomponents = 2; break; + default: _error_("not supported yet"); + } + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input=element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* rheologyb_input=element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* rheologyb0_input=element->GetInput(RheologyBInitialguessEnum); _assert_(rheologyb0_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,RheologyBInitialguessMisfitEnum); + rheologyb_input->GetInputValue(&B,gauss); + rheologyb0_input->GetInputValue(&B0,gauss); + + /*Tikhonov regularization: J = 1/2 (B-B0)^2 */ + Jelem+=weight*1./2.*(B-B0)*(B-B0)*Jdet*gauss->weight; + } + + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; + return Jelem; +} Index: /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.h (revision 23866) +++ /issm/trunk-jpl/src/c/modules/RheologyBAbsGradientx/RheologyBAbsGradientx.h (revision 23867) @@ -12,3 +12,6 @@ IssmDouble RheologyBAbsGradient(Element* element); +void RheologyBInitialguessMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); +IssmDouble RheologyBInitialguessMisfit(Element* element); + #endif Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23866) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23867) @@ -601,4 +601,6 @@ RheologyBAbsGradientEnum, RheologyBbarAbsGradientEnum, + RheologyBInitialguessMisfitEnum, + RheologyBInitialguessEnum, SealevelEnum, SealevelriseCumDeltathicknessEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23867) @@ -607,4 +607,6 @@ case RheologyBAbsGradientEnum : return "RheologyBAbsGradient"; case RheologyBbarAbsGradientEnum : return "RheologyBbarAbsGradient"; + case RheologyBInitialguessMisfitEnum : return "RheologyBInitialguessMisfit"; + case RheologyBInitialguessEnum : return "RheologyBInitialguess"; case SealevelEnum : return "Sealevel"; case SealevelriseCumDeltathicknessEnum : return "SealevelriseCumDeltathickness"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23866) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23867) @@ -619,4 +619,6 @@ else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; + else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else if (strcmp(name,"SealevelriseCumDeltathickness")==0) return SealevelriseCumDeltathicknessEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum; else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; - else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; - else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; + if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; + else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum; + else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum; else if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; - else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; - else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; + if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; + else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; + else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; - else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; - else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum; + if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; + else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; + else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum; else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum; else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; @@ -996,10 +998,10 @@ else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum; - else if (strcmp(name,"Fset")==0) return FsetEnum; - else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"FSSolver")==0) return FSSolverEnum; + if (strcmp(name,"Fset")==0) return FsetEnum; + else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; + else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; - else if (strcmp(name,"None")==0) return NoneEnum; - else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + if (strcmp(name,"None")==0) return NoneEnum; + else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; + else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"XY")==0) return XYEnum; else if (strcmp(name,"XYZ")==0) return XYZEnum; - else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum; - else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; + if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum; + else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum; + else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum; else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;