Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 22264) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp (revision 22265) @@ -90,4 +90,10 @@ case SurfaceAbsMisfitEnum: for(i=0;ivalues[i]+=(surfaceobs-surface)*weight*Jdet*gauss->weight*basis[i]; + break; + case OmegaAbsGradientEnum: + /*Nothing in P vector*/ + break; + case EtaDiffEnum: + /*Nothing in P vector*/ break; default: @@ -132,4 +138,6 @@ /*Nothing, \partial J/\partial k = 0*/ break; + case OmegaAbsGradientEnum: GradientJOmegaGradient(element,gradient,control_index); break; + case EtaDiffEnum: GradientJEtaDiff(element,gradient,control_index); break; default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); } @@ -193,5 +201,5 @@ /*Intermediaries*/ int n=3; - IssmDouble dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet; + IssmDouble dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet,velobs; IssmDouble *xyz_list= NULL; @@ -212,5 +220,5 @@ Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); - IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum); + Input* velobs_input = element->GetInput(InversionVelObsEnum); _assert_(velobs_input); Gauss* gauss=element->NewGauss(2); @@ -226,10 +234,12 @@ surfaceslopex_input->GetInputValue(&slopex,gauss); surfaceslopey_input->GetInputValue(&slopey,gauss); + velobs_input->GetInputValue(&velobs,gauss); slope = sqrt(slopex*slopex + slopey*slopey); + //if(slope<1.e-5) slope = 1.e-5; /*Build gradient vector (actually -dJ/da): */ for(int i=0;iweight*basis[i]*pow(rhog,n)*exp(omega)*pow(slope,n-1)*(ds[0]*dlambda[0] + ds[1]*dlambda[1]); + ge[i]+= - Jdet*gauss->weight*basis[i]*velobs/slope*(ds[0]*dlambda[0] + ds[1]*dlambda[1]); _assert_(!xIsNan(ge[i])); } @@ -243,4 +253,104 @@ xDelete(vertexpidlist); delete gauss; +}/*}}}*/ +void AdjointBalancethickness2Analysis::GradientJOmegaGradient(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble dk[3]; + IssmDouble *xyz_list= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* dbasis = xNew(2*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* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_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->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss); + weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum); + + /*Build alpha_complement_list: */ + omega_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss); + + /*Build gradient vector (actually -dJ/ddrag): */ + for(int i=0;iweight*2*(dk[0]*dk[0] + dk[1]*dk[1])*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(dbasis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + +}/*}}}*/ +void AdjointBalancethickness2Analysis::GradientJEtaDiff(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble omega,omega0; + 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); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); + Input* omega0_input = element->GetInput(BalancethicknessOmega0Enum); _assert_(omega0_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,EtaDiffEnum); + + /*Build alpha_complement_list: */ + omega_input->GetInputValue(&omega,gauss); + omega0_input->GetInputValue(&omega0,gauss); + + /*Build gradient vector (actually -dJ/ddrag): */ + for(int i=0;iweight*(omega - omega0)*basis[i]; + _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 AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 22264) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h (revision 22265) @@ -30,4 +30,6 @@ void GradientJdHdt(Element* element,Vector* gradient,int control_index); void GradientJOmega(Element* element,Vector* gradient,int control_index); + void GradientJOmegaGradient(Element* element,Vector* gradient,int control_index); + void GradientJEtaDiff(Element* element,Vector* gradient,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 22264) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 22265) @@ -176,4 +176,5 @@ if(control_type!=VxEnum && control_type!=VyEnum && + control_type!=BalancethicknessSpcthicknessEnum && control_type!=BalancethicknessThickeningRateEnum){ _error_("Control "<(responses); +}/*}}}*/ +void AdjointBalancethicknessAnalysis::GradientJDirichlet(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* ge = xNewZeroInit(numvertices); + IssmDouble* lambda = xNew(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GradientIndexing(&vertexpidlist[0],control_index); + element->GetInputListOnVertices(lambda,AdjointEnum); + + BalancethicknessAnalysis* analysis = new BalancethicknessAnalysis(); + ElementMatrix* Ke = analysis->CreateKMatrix(element); + delete analysis; + + /*Transpose and return Ke*/ + Ke->Transpose(); + _assert_(Ke->nrows == numvertices && Ke->ncols == numvertices); + + for(int i=0;ivalues[i*Ke->ncols + j] * lambda[j]; + } + //ge[i]=-lambda[i]; + _assert_(!xIsNan(ge[i])); + } + gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL); + + /*Clean up and return*/ + xDelete(ge); + xDelete(lambda); + xDelete(vertexpidlist); + delete Ke; }/*}}}*/ void AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector* gradient,int control_index){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 22264) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 22265) @@ -28,4 +28,5 @@ void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void GradientJDirichlet(Element* element,Vector* gradient,int control_index); void GradientJDhDt(Element* element,Vector* gradient,int control_index); void GradientJVx(Element* element,Vector* gradient,int control_index); Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 22264) +++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 22265) @@ -35,10 +35,9 @@ iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0); iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum); + iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum); iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum); iomodel->FetchDataToInput(elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum); - iomodel->FetchDataToInput(elements,"md.balancethickness.omega",BalancethicknessOmegaEnum); - iomodel->FetchDataToInput(elements,"md.balancethickness.slopex",SurfaceSlopeXEnum); - iomodel->FetchDataToInput(elements,"md.balancethickness.slopey",SurfaceSlopeYEnum); /*Update elements: */ @@ -71,6 +70,5 @@ /*Intermediaries */ - IssmDouble yts = 365*24*3600.; - IssmDouble Jdet,vx,vy,vel; + IssmDouble Jdet,ds[2],slope,velobs,omega; IssmDouble* xyz_list = NULL; @@ -84,18 +82,8 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - Input* vx_input = element->GetInput(SurfaceSlopeXEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(SurfaceSlopeYEnum); _assert_(vy_input); - - /*make sure are diffusivisty is large enough*/ - vel = sqrt(vx*vx+vy*vy); - if(sqrt(vx*vx+vy*vy)==0.){ - vx = 0.1/yts; - vy = 0.1/yts; - } - else if(vel<0.1/yts){ - vx = vx/vel*0.1; - vy = vy/vel*0.1; - - } + Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input); + Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + Input* velobs_input = element->GetInput(InversionVelObsEnum); _assert_(velobs_input); /* Start looping on the number of gaussian points: */ @@ -105,10 +93,15 @@ element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); element->JacobianDeterminant(&Jdet,xyz_list,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + surfaceslopex_input->GetInputValue(&ds[0],gauss); + surfaceslopey_input->GetInputValue(&ds[1],gauss); + velobs_input->GetInputValue(&velobs,gauss); + omega_input->GetInputValue(&omega,gauss); + + slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]); + //if(slope<1.e-5) slope = 1.e-5; for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*(vx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + Ke->values[i*numnodes+j] += velobs/slope*omega*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); } } @@ -124,5 +117,5 @@ /*Intermediaries */ - IssmDouble dhdt[2],mb[2],ms[2],Jdet; + IssmDouble dhdt,mb,ms,Jdet; IssmDouble* xyz_list = NULL; @@ -148,12 +141,10 @@ element->NodalFunctions(basis,gauss); - ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss); - - ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss); - mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss); - dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss); + ms_input->GetInputValue(&ms,gauss); + mb_input->GetInputValue(&mb,gauss); + dhdt_input->GetInputValue(&dhdt,gauss); for(int i=0;ivalues[i]+=Jdet*gauss->weight*( - (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i] + (ms-mb-dhdt)*basis[i] ); } @@ -166,5 +157,5 @@ }/*}}}*/ void Balancethickness2Analysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - element->GetSolutionFromInputsOneDof(solution,ThicknessEnum); + element->GetSolutionFromInputsOneDof(solution,SurfaceEnum); }/*}}}*/ void Balancethickness2Analysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ @@ -173,6 +164,55 @@ void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum); - + /*Intermediaries*/ + IssmDouble ds[2],s,b,D; + IssmDouble* xyz_list = NULL; + + //element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum); + + /*Fetch number of vertices and allocate velocity vectors*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* vel_list = xNew(numvertices); + IssmDouble* vx_list = xNew(numvertices); + IssmDouble* vy_list = xNew(numvertices); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* s_input = element->GetInput(SurfaceEnum); _assert_(s_input); + Input* b_input = element->GetInput(BaseEnum); _assert_(b_input); + + /*Calculate velocities*/ + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + + if(D_input){ + D_input->GetInputValue(&D,gauss); + } + else{ + D = 0.; + } + b_input->GetInputValue(&b,gauss); + s_input->GetInputValue(&s,gauss); + s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss); + + vx_list[iv] = -1./(s-b)*D*ds[0]; + vy_list[iv] = -1./(s-b)*D*ds[1]; + vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2)); + } + + /*Add vx and vy as inputs to the tria element: */ + element->AddInput(VxEnum,vx_list,P1Enum); + element->AddInput(VyEnum,vy_list,P1Enum); + element->AddInput(VelEnum,vel_list,P1Enum); + + /*Free ressources:*/ + delete gauss; + xDelete(vy_list); + xDelete(vx_list); + xDelete(vel_list); + xDelete(xyz_list); }/*}}}*/ void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/