Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 18604) @@ -195,25 +195,44 @@ void AdjointBalancethickness2Analysis::GradientJdHdt(Element* element,Vector* gradient,int control_index){/*{{{*/ + /*Intermediaries*/ + IssmDouble lambda,Jdet; + IssmDouble *xyz_list= NULL; + /*Fetch number of vertices for this finite element*/ int numvertices = element->GetNumberOfVertices(); /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); IssmDouble* ge = xNewZeroInit(numvertices); - IssmDouble* lambda = xNew(numvertices); int* vertexpidlist = xNew(numvertices); /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); element->GradientIndexing(&vertexpidlist[0],control_index); - element->GetInputListOnVertices(lambda,AdjointEnum); - for(int i=0;i(ge[i])); - } - gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL); + Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input); + + 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); + + adjoint_input->GetInputValue(&lambda,gauss); + + /*Build gradient vector (actually -dJ/da): */ + for(int i=0;iweight*basis[i]*lambda; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); /*Clean up and return*/ xDelete(ge); - xDelete(lambda); + xDelete(xyz_list); + xDelete(basis); xDelete(vertexpidlist); + delete gauss; }/*}}}*/ void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 18604) @@ -82,8 +82,8 @@ Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); Input* D0_input = element->GetInput(BalancethicknessD0Enum); - if(!D0_input){ + //if(!D0_input){ this->CreateD0(element); D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input); - } + //} /* Start looping on the number of gaussian points: */ @@ -220,5 +220,5 @@ /*Intermediaries */ - IssmDouble Gamma,h,mu0,ds[2],Cmu,B,k,s,b; + IssmDouble Gamma,h,mu0,ds[2],Cmu,B,k,s,b,normds; const int n = 3; const IssmDouble Hstar = 500.; @@ -232,7 +232,9 @@ /*retrieve what we need: */ element->GetVerticesCoordinates(&xyz_list); - Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); - Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); + Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); + Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); /*Calculate damage evolution source term: */ @@ -242,5 +244,22 @@ B_input->GetInputValue(&B,gauss); - surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + if(surfaceslopex_input && surfaceslopey_input){ + surfaceslopex_input->GetInputValue(&ds[0],gauss); + surfaceslopey_input->GetInputValue(&ds[1],gauss); + } + else{ + surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + } + + /*check slopes*/ + normds = sqrt(ds[0]*ds[0]+ds[1]*ds[1]); + if (normds==0.){ + _error_("surface slope is zero"); + } + if(normds<1.e-5){ + ds[0] = ds[0]/normds*1.e+5; + ds[1] = ds[1]/normds*1.e+5; + normds = 1.e-5; + } mu0 = pow(2.,(1-3*n)/(2.*n)) * B; @@ -252,4 +271,8 @@ /*Add input*/ element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType()); + //if(surfaceslopex_input && surfaceslopey_input){ + // element->DeleteInput(SurfaceSlopeXEnum); + // element->DeleteInput(SurfaceSlopeYEnum); + //} /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 18604) @@ -48,5 +48,5 @@ /* Intermediaries */ int domaintype; - IssmDouble Jdet,thickness,l=8.; + IssmDouble Jdet,thickness,l; IssmDouble *xyz_list = NULL; @@ -68,4 +68,5 @@ /*Retrieve all inputs and parameters*/ + element->FindParam(&l,SmoothThicknessMultiplierEnum); _assert_(l>0.); element->GetVerticesCoordinates(&xyz_list); Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp (revision 18604) @@ -15,11 +15,17 @@ bool save_results; - /*activate formulation: */ - femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum); - /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n"); + femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum); + femmodel->SetCurrentConfiguration(SmoothAnalysisEnum); + femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToSmoothEnum); + solutionsequence_linear(femmodel); + femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum); + solutionsequence_linear(femmodel); + if(VerboseSolution()) _printf0_("call computational core:\n"); + femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum); solutionsequence_linear(femmodel); //solutionsequence_nonlinear(femmodel,false); @@ -27,6 +33,8 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - const int numoutputs = 4; - int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum}; + const int numoutputs = 6; + int outputs[numoutputs] = {SurfaceEnum,SurfaceSlopeXEnum,SurfaceSlopeYEnum,VxEnum,VyEnum,VelEnum}; + //const int numoutputs = 4; + //int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum}; femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs); } Index: /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp (revision 18604) @@ -18,5 +18,6 @@ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - if(VerboseSolution()) _printf0_("computing smoothed slopes:\n"); + if(VerboseSolution()) _printf0_("computing smooth driving stress:\n"); + femmodel->parameters->SetParam(8.,SmoothThicknessMultiplierEnum); femmodel->SetCurrentConfiguration(SmoothAnalysisEnum); femmodel->parameters->SetParam(DrivingStressXEnum,InputToSmoothEnum); Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 18604) @@ -13,5 +13,5 @@ int solution_type,n; int num_responses; - IssmDouble j0,j; + IssmDouble j0,j,yts; IssmDouble Ialpha,exponent,alpha; IssmDouble* jlist = NULL; @@ -28,4 +28,5 @@ femmodel->parameters->SetParam(false,SaveResultsEnum); femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); /*Get initial guess*/ @@ -64,5 +65,5 @@ IssmDouble exp0 = 0.; IssmDouble incr = -0.2; - IssmDouble exp1 = -10.; + IssmDouble exp1 = -18.; int num = reCast((exp1-exp0)/incr); @@ -88,5 +89,5 @@ IssmDouble Den = 0.; - for(int i=0;iFetchDataToInput(elements,InversionVyObsEnum,0.); iomodel->FetchDataToInput(elements,InversionThicknessObsEnum,0.); - //iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.); + iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.); iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18603) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18604) @@ -450,4 +450,5 @@ InputToDepthaverageEnum, InputToSmoothEnum, + SmoothThicknessMultiplierEnum, IntParamEnum, IntVecParamEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18604) @@ -447,4 +447,5 @@ case InputToDepthaverageEnum : return "InputToDepthaverage"; case InputToSmoothEnum : return "InputToSmooth"; + case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier"; case IntParamEnum : return "IntParam"; case IntVecParamEnum : return "IntVecParam"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18603) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18604) @@ -456,4 +456,5 @@ else if (strcmp(name,"InputToDepthaverage")==0) return InputToDepthaverageEnum; else if (strcmp(name,"InputToSmooth")==0) return InputToSmoothEnum; + else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; - else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Boundary")==0) return BoundaryEnum; + 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; @@ -628,9 +629,9 @@ else if (strcmp(name,"Time")==0) return TimeEnum; else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; - else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; + if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; + else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum; else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum; - else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum; + if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum; + else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum; else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum; else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18603) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18604) @@ -439,4 +439,5 @@ def InputToDepthaverageEnum(): return StringToEnum("InputToDepthaverage")[0] def InputToSmoothEnum(): return StringToEnum("InputToSmooth")[0] +def SmoothThicknessMultiplierEnum(): return StringToEnum("SmoothThicknessMultiplier")[0] def IntParamEnum(): return StringToEnum("IntParam")[0] def IntVecParamEnum(): return StringToEnum("IntVecParam")[0] Index: /issm/trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m (revision 18604) +++ /issm/trunk-jpl/src/m/enum/SmoothThicknessMultiplierEnum.m (revision 18604) @@ -0,0 +1,11 @@ +function macro=SmoothThicknessMultiplierEnum() +%SMOOTHTHICKNESSMULTIPLIERENUM - Enum of SmoothThicknessMultiplier +% +% 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=SmoothThicknessMultiplierEnum() + +macro=StringToEnum('SmoothThicknessMultiplier'); Index: sm/trunk-jpl/src/m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmoothedSurfaceSlopeXAnalysisEnum.m (revision 18603) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SmoothedSurfaceSlopeXAnalysisEnum() -%SMOOTHEDSURFACESLOPEXANALYSISENUM - Enum of SmoothedSurfaceSlopeXAnalysis -% -% 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=SmoothedSurfaceSlopeXAnalysisEnum() - -macro=StringToEnum('SmoothedSurfaceSlopeXAnalysis'); Index: sm/trunk-jpl/src/m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmoothedSurfaceSlopeYAnalysisEnum.m (revision 18603) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=SmoothedSurfaceSlopeYAnalysisEnum() -%SMOOTHEDSURFACESLOPEYANALYSISENUM - Enum of SmoothedSurfaceSlopeYAnalysis -% -% 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=SmoothedSurfaceSlopeYAnalysisEnum() - -macro=StringToEnum('SmoothedSurfaceSlopeYAnalysis');