Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18565) @@ -45,5 +45,63 @@ }/*}}}*/ ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/ - _error_("Not implemented yet"); + + /*Intermediaries */ + int num_responses,i; + IssmDouble vx,vy,vel,Jdet; + IssmDouble surface,surfaceobs,weight; + int *responses = NULL; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and vectors*/ + ElementVector* pe = element->NewElementVector(SSAApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(2*numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&num_responses,InversionNumCostFunctionsEnum); + element->FindParam(&responses,NULL,InversionCostFunctionsEnum); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* surfaceobs_input = element->GetInput(InversionSurfaceObsEnum); _assert_(surfaceobs_input); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_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->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + surface_input->GetInputValue(&surface, gauss); + surfaceobs_input->GetInputValue(&surfaceobs, gauss); + + /*Loop over all requested responses*/ + for(int resp=0;respGetInputValue(&weight,gauss,responses[resp]); + + switch(responses[resp]){ + case SurfaceAbsMisfitEnum: + for(i=0;ivalues[i]+=(surfaceobs-surface)*weight*Jdet*gauss->weight*basis[i]; + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + + /*Clean up and return*/ + xDelete(responses); + xDelete(xyz_list); + xDelete(basis); + xDelete(dbasis); + delete gauss; + return pe; }/*}}}*/ @@ -71,5 +129,5 @@ /*Deal with first part (partial derivative a J with respect to k)*/ for(resp=0;resp* gradient,int control_index){/*{{{*/ +void AdjointBalancethickness2Analysis::GradientJOmega(Element* element,Vector* gradient,int control_index){/*{{{*/ /*Intermediaries*/ - IssmDouble lambda,Jdet; + IssmDouble dlambda[2],ds[2],D0,Jdet; IssmDouble *xyz_list= NULL; @@ -103,5 +162,7 @@ element->GetVerticesCoordinates(&xyz_list); element->GradientIndexing(&vertexpidlist[0],control_index); - Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); + Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); + Input* s_input = element->GetInput(SurfaceEnum); _assert_(s_input); + Input* D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input); Gauss* gauss=element->NewGauss(2); @@ -111,9 +172,12 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctionsP1(basis,gauss); - adjoint_input->GetInputValue(&lambda,gauss); + + D0_input->GetInputValue(&D0,gauss); + adjoint_input->GetInputDerivativeValue(&dlambda[0],xyz_list,gauss); + s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); /*Build gradient vector (actually -dJ/da): */ for(int i=0;iweight*basis[i]*lambda; + ge[i]+= Jdet*gauss->weight*basis[i]*(ds[0]*dlambda[0] + ds[1]*dlambda[1]); _assert_(!xIsNan(ge[i])); } Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 18564) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 18565) @@ -28,5 +28,5 @@ void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); - void GradientJAdot(Element* element,Vector* gradient,int control_index); + void GradientJOmega(Element* element,Vector* gradient,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18565) @@ -15,4 +15,14 @@ /*Finite element type*/ int finiteelement = P1Enum; + + /*Load variables in element*/ + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BaseEnum); + iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); + iomodel->FetchDataToInput(elements,BalancethicknessOmegaEnum); /*Update elements: */ @@ -22,17 +32,9 @@ Element* element=(Element*)elements->GetObjectByOffset(counter); element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; } } - iomodel->FetchDataToInput(elements,ThicknessEnum); - iomodel->FetchDataToInput(elements,SurfaceEnum); - iomodel->FetchDataToInput(elements,BaseEnum); - iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); - iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); - iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); - iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum); - iomodel->FetchDataToInput(elements,BalancethicknessCmuEnum); }/*}}}*/ void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ @@ -66,5 +68,5 @@ /*Intermediaries */ - IssmDouble Jdet,D; + IssmDouble Jdet,D0,omega; IssmDouble* xyz_list = NULL; @@ -76,10 +78,12 @@ IssmDouble* dbasis = xNew(2*numnodes); - /*Create input D*/ - this->CreateDiffusionCoefficient(element); - /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input); + Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); + Input* D0_input = element->GetInput(BalancethicknessD0Enum); + if(!D0_input){ + this->CreateD0(element); + D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input); + } /* Start looping on the number of gaussian points: */ @@ -89,9 +93,10 @@ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); element->JacobianDeterminant(&Jdet,xyz_list,gauss); - D_input->GetInputValue(&D,gauss); + D0_input->GetInputValue(&D0,gauss); + omega_input->GetInputValue(&omega,gauss); for(int i=0;ivalues[i*numnodes+j] += D*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + Ke->values[i*numnodes+j] += D0*omega*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); } } @@ -107,5 +112,5 @@ /*Intermediaries */ - IssmDouble dhdt,mb,ms,D,db[2],Jdet; + IssmDouble dhdt,mb,ms,Jdet; IssmDouble* xyz_list = NULL; @@ -116,5 +121,4 @@ ElementVector* pe = element->NewElementVector(); IssmDouble* basis = xNew(numnodes); - IssmDouble* dbasis = xNew(2*numnodes); /*Retrieve all inputs and parameters*/ @@ -123,6 +127,4 @@ Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input); Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); - Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input); - Input* bed_input = element->GetInput(BaseEnum); _assert_(bed_input); /* Start looping on the number of gaussian points: */ @@ -133,17 +135,11 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); ms_input->GetInputValue(&ms,gauss); mb_input->GetInputValue(&mb,gauss); dhdt_input->GetInputValue(&dhdt,gauss); - D_input->GetInputValue(&D,gauss); - bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss); - //db[0]=0.; - //db[1]=0.; for(int i=0;ivalues[i]+=Jdet*gauss->weight*( (ms-mb-dhdt)*basis[i] - //-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i] ); } @@ -220,5 +216,5 @@ /*Specifics*/ -void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/ +void Balancethickness2Analysis::CreateD0(Element* element){/*{{{*/ /*Intermediaries */ @@ -231,14 +227,10 @@ /*Fetch number of vertices and allocate output*/ int numnodes = element->GetNumberOfNodes(); - IssmDouble* D = xNew(numnodes); - IssmDouble* Dgradb = xNew(numnodes); + IssmDouble* D0 = xNew(numnodes); /*retrieve what we need: */ element->GetVerticesCoordinates(&xyz_list); - Input* bed_input = element->GetInput(BaseEnum); _assert_(bed_input); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - Input* k_input = element->GetInput(FrictionCoefficientEnum);_assert_(k_input); - Input* Cmu_input = element->GetInput(BalancethicknessCmuEnum);_assert_(Cmu_input); IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); @@ -249,24 +241,17 @@ B_input->GetInputValue(&B,gauss); - k_input->GetInputValue(&k,gauss); - //thickness_input->GetInputValue(&h,gauss); surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); - surface_input->GetInputValue(&s,gauss); - bed_input->GetInputValue(&b,gauss); - h = s-b; mu0 = pow(2.,(1-3*n)/(2.*n)) * B; Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n); - //Cmu = pow(mu0,n)/Hstar/(k*k); - Cmu_input->GetInputValue(&Cmu,gauss); - - D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h + Cmu); + + D0[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)/(n+2); } /*Add input*/ - element->AddInput(BalancethicknessDiffusionCoefficientEnum,D,element->GetElementType()); + element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType()); /*Clean up and return*/ - xDelete(D); + xDelete(D0); delete gauss; }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 18564) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h (revision 18565) @@ -32,5 +32,5 @@ /*Specifics*/ - void CreateDiffusionCoefficient(Element* element); + void CreateD0(Element* element); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18565) @@ -1034,4 +1034,5 @@ /*No yts conversion*/ case ThicknessEnum: + case BalancethicknessOmegaEnum: case FrictionCoefficientEnum: if(iomodel->Data(control)){ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18565) @@ -641,5 +641,4 @@ 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; @@ -719,5 +718,5 @@ case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; case BalancethicknessMisfitEnum: BalancethicknessMisfitx(&double_result); break; - case Balancethickness2MisfitEnum: Balancethickness2Misfitx(&double_result); break; + case SurfaceAbsMisfitEnum: SurfaceAbsMisfitx(&double_result); break; /*Vector */ @@ -1444,5 +1443,5 @@ }/*}}}*/ -void FemModel::Balancethickness2Misfitx(IssmDouble* presponse){/*{{{*/ +void FemModel::SurfaceAbsMisfitx(IssmDouble* presponse){/*{{{*/ /*output: */ @@ -1450,6 +1449,5 @@ IssmDouble J_sum; - IssmDouble weight,thicknessobs,thickness,potential,dpotential[2]; - IssmDouble vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nu; + IssmDouble surface,surfaceobs,weight; IssmDouble Jdet; IssmDouble* xyz_list = NULL; @@ -1461,5 +1459,30 @@ /*If on water, return 0: */ if(!element->IsIceInElement()) continue; - _error_("Not implemented"); + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + Input* weights_input =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* surfaceobs_input=element->GetInput(InversionSurfaceObsEnum); _assert_(surfaceobs_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,SurfaceAbsMisfitEnum); + surface_input->GetInputValue(&surface,gauss); + surfaceobs_input->GetInputValue(&surfaceobs,gauss); + + /*Compute SurfaceAbsMisfitEnum*/ + J+=0.5*(surface-surfaceobs)*(surface-surfaceobs)*weight*Jdet*gauss->weight; + } } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18564) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18565) @@ -81,5 +81,4 @@ 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); @@ -93,4 +92,5 @@ void CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn); void ThicknessAbsGradientx( IssmDouble* pJ); + void SurfaceAbsMisfitx( IssmDouble* pJ); #ifdef _HAVE_GIA_ void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18565) @@ -22,12 +22,11 @@ if(VerboseSolution()) _printf0_("call computational core:\n"); - //solutionsequence_linear(femmodel); - solutionsequence_nonlinear(femmodel,false); + solutionsequence_linear(femmodel); + //solutionsequence_nonlinear(femmodel,false); if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - const int numoutputs = 5; - //int outputs[numoutputs] = {ThicknessEnum}; - int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum,BalancethicknessDiffusionCoefficientEnum}; + const int numoutputs = 4; + int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum}; femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 18565) @@ -29,4 +29,5 @@ iomodel->FetchDataToInput(elements,InversionVyObsEnum,0.); iomodel->FetchDataToInput(elements,InversionThicknessObsEnum,0.); + iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.); iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum); @@ -42,4 +43,5 @@ case FrictionCoefficientEnum: case BalancethicknessApparentMassbalanceEnum: + case BalancethicknessOmegaEnum: iomodel->FetchData(1,control); break; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18564) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18565) @@ -165,4 +165,5 @@ InversionStepThresholdEnum, InversionThicknessObsEnum, + InversionSurfaceObsEnum, InversionVxObsEnum, InversionVyObsEnum, @@ -307,4 +308,6 @@ BalancethicknessDiffusionCoefficientEnum, BalancethicknessCmuEnum, + BalancethicknessOmegaEnum, + BalancethicknessD0Enum, /*}}}*/ /*Surfaceforcings{{{*/ @@ -527,4 +530,5 @@ TemperaturePicardEnum, ThicknessAbsMisfitEnum, + SurfaceAbsMisfitEnum, VelEnum, VelocityEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18565) @@ -173,4 +173,5 @@ case InversionStepThresholdEnum : return "InversionStepThreshold"; case InversionThicknessObsEnum : return "InversionThicknessObs"; + case InversionSurfaceObsEnum : return "InversionSurfaceObs"; case InversionVxObsEnum : return "InversionVxObs"; case InversionVyObsEnum : return "InversionVyObs"; @@ -315,4 +316,6 @@ case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient"; case BalancethicknessCmuEnum : return "BalancethicknessCmu"; + case BalancethicknessOmegaEnum : return "BalancethicknessOmega"; + case BalancethicknessD0Enum : return "BalancethicknessD0"; case SurfaceforcingsEnum : return "Surfaceforcings"; case SMBEnum : return "SMB"; @@ -518,4 +521,5 @@ case TemperaturePicardEnum : return "TemperaturePicard"; case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit"; + case SurfaceAbsMisfitEnum : return "SurfaceAbsMisfit"; case VelEnum : return "Vel"; case VelocityEnum : return "Velocity"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18564) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18565) @@ -176,4 +176,5 @@ else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; + else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; - else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum; + if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; + else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum; else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum; else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum; @@ -321,4 +322,6 @@ else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum; else if (strcmp(name,"BalancethicknessCmu")==0) return BalancethicknessCmuEnum; + else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum; + else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum; else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum; else if (strcmp(name,"SMB")==0) return SMBEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; - else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; - else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; - else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; + if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; + else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; + else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; + else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else if (strcmp(name,"Boundary")==0) return BoundaryEnum; - else if (strcmp(name,"Converged")==0) return ConvergedEnum; - else if (strcmp(name,"Fill")==0) return FillEnum; - else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Friction")==0) return FrictionEnum; + if (strcmp(name,"Converged")==0) return ConvergedEnum; + else if (strcmp(name,"Fill")==0) return FillEnum; + else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; + else if (strcmp(name,"Friction")==0) return FrictionEnum; else if (strcmp(name,"Internal")==0) return InternalEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; @@ -530,4 +533,5 @@ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; + else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; @@ -625,12 +629,12 @@ else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum; - else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum; else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum; else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum; else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum; + else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum; else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum; else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum; @@ -748,12 +752,12 @@ else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum; else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum; - else if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum; else if (strcmp(name,"Matseaice")==0) return MatseaiceEnum; else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum; else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum; + else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum; else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum; else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum; Index: /issm/trunk-jpl/src/m/classes/balancethickness.m =================================================================== --- /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 18564) +++ /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 18565) @@ -10,5 +10,5 @@ stabilization = 0; - Cmu = NaN; + omega = NaN; end methods @@ -35,5 +35,5 @@ md = checkfield(md,'fieldname','balancethickness.stabilization','size',[1 1],'values',[0 1 2 3]); - %md = checkfield(md,'fieldname','balancethickness.Cmu','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0); + %md = checkfield(md,'fieldname','balancethickness.omega','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0); end % }}} function disp(obj) % {{{ @@ -53,5 +53,5 @@ WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer'); - WriteData(fid,'object',obj,'fieldname','Cmu','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','omega','format','DoubleMat','mattype',1); end % }}} end Index: /issm/trunk-jpl/src/m/classes/m1qn3inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/m1qn3inversion.m (revision 18564) +++ /issm/trunk-jpl/src/m/classes/m1qn3inversion.m (revision 18565) @@ -22,4 +22,5 @@ vel_obs = NaN thickness_obs = NaN + surface_obs = NaN end @@ -72,5 +73,5 @@ md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',... {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar',... - 'Vx' 'Vy' 'Thickness' 'BalancethicknessNu' 'BalancethicknessApparentMassbalance'}); + 'Vx' 'Vy' 'Thickness' 'BalancethicknessOmega' 'BalancethicknessApparentMassbalance'}); md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0); md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0); @@ -84,4 +85,5 @@ if solution==BalancethicknessSolutionEnum() md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); + md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1); elseif solution==BalancethicknessSoftSolutionEnum() md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1); @@ -110,4 +112,5 @@ fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]'); fielddisplay(obj,'thickness_obs','observed thickness [m]'); + fielddisplay(obj,'surface_obs','observed surface elevation [m]'); disp('Available cost functions:'); disp(' 101: SurfaceAbsVelMisfit'); @@ -145,4 +148,5 @@ end WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype); + WriteData(fid,'object',obj,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype); %process control parameters @@ -170,4 +174,5 @@ pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum(); pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum(); + pos=find(obj.cost_functions==601); data(pos)=SurfaceAbsMisfitEnum(); WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3); WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer'); Index: /issm/trunk-jpl/src/m/enum/BalancethicknessD0Enum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessD0Enum.m (revision 18565) +++ /issm/trunk-jpl/src/m/enum/BalancethicknessD0Enum.m (revision 18565) @@ -0,0 +1,11 @@ +function macro=BalancethicknessD0Enum() +%BALANCETHICKNESSD0ENUM - Enum of BalancethicknessD0 +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=BalancethicknessD0Enum() + +macro=StringToEnum('BalancethicknessD0'); Index: /issm/trunk-jpl/src/m/enum/BalancethicknessOmegaEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/BalancethicknessOmegaEnum.m (revision 18565) +++ /issm/trunk-jpl/src/m/enum/BalancethicknessOmegaEnum.m (revision 18565) @@ -0,0 +1,11 @@ +function macro=BalancethicknessOmegaEnum() +%BALANCETHICKNESSOMEGAENUM - Enum of BalancethicknessOmega +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=BalancethicknessOmegaEnum() + +macro=StringToEnum('BalancethicknessOmega'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18564) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18565) @@ -165,4 +165,5 @@ def InversionStepThresholdEnum(): return StringToEnum("InversionStepThreshold")[0] def InversionThicknessObsEnum(): return StringToEnum("InversionThicknessObs")[0] +def InversionSurfaceObsEnum(): return StringToEnum("InversionSurfaceObs")[0] def InversionVxObsEnum(): return StringToEnum("InversionVxObs")[0] def InversionVyObsEnum(): return StringToEnum("InversionVyObs")[0] @@ -307,4 +308,6 @@ def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0] def BalancethicknessCmuEnum(): return StringToEnum("BalancethicknessCmu")[0] +def BalancethicknessOmegaEnum(): return StringToEnum("BalancethicknessOmega")[0] +def BalancethicknessD0Enum(): return StringToEnum("BalancethicknessD0")[0] def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0] def SMBEnum(): return StringToEnum("SMB")[0] @@ -510,4 +513,5 @@ def TemperaturePicardEnum(): return StringToEnum("TemperaturePicard")[0] def ThicknessAbsMisfitEnum(): return StringToEnum("ThicknessAbsMisfit")[0] +def SurfaceAbsMisfitEnum(): return StringToEnum("SurfaceAbsMisfit")[0] def VelEnum(): return StringToEnum("Vel")[0] def VelocityEnum(): return StringToEnum("Velocity")[0] Index: /issm/trunk-jpl/src/m/enum/InversionSurfaceObsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/InversionSurfaceObsEnum.m (revision 18565) +++ /issm/trunk-jpl/src/m/enum/InversionSurfaceObsEnum.m (revision 18565) @@ -0,0 +1,11 @@ +function macro=InversionSurfaceObsEnum() +%INVERSIONSURFACEOBSENUM - Enum of InversionSurfaceObs +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=InversionSurfaceObsEnum() + +macro=StringToEnum('InversionSurfaceObs'); Index: /issm/trunk-jpl/src/m/enum/SurfaceAbsMisfitEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceAbsMisfitEnum.m (revision 18565) +++ /issm/trunk-jpl/src/m/enum/SurfaceAbsMisfitEnum.m (revision 18565) @@ -0,0 +1,11 @@ +function macro=SurfaceAbsMisfitEnum() +%SURFACEABSMISFITENUM - Enum of SurfaceAbsMisfit +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SurfaceAbsMisfitEnum() + +macro=StringToEnum('SurfaceAbsMisfit');