Index: /issm/trunk/src/c/Container/Inputs.cpp =================================================================== --- /issm/trunk/src/c/Container/Inputs.cpp (revision 5639) +++ /issm/trunk/src/c/Container/Inputs.cpp (revision 5640) @@ -71,4 +71,32 @@ } /*}}}*/ +/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ +void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){ + + vector::iterator object; + Input* input=NULL; + bool found=false; + + /*Go through inputs and check whether any input with the same name is already in: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + input=(Input*)(*object); + if (input->EnumType()==enum_type){ + found=true; + break; + } + } + + if (!found){ + /*we could not find an input with the correct enum type. No defaults values were provided, + * error out: */ + ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type)); + } + + /*Ok, we have an input if we made it here, request the input to return the values: */ + input->GetParameterValue(pvalue,gauss); + +} +/*}}}*/ /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/ void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){ Index: /issm/trunk/src/c/Container/Inputs.h =================================================================== --- /issm/trunk/src/c/Container/Inputs.h (revision 5639) +++ /issm/trunk/src/c/Container/Inputs.h (revision 5640) @@ -16,4 +16,5 @@ class Input; class Node; +class GaussTria; class Inputs: public DataSet{ @@ -43,4 +44,5 @@ void GetParameterValue(double* pvalue,int enum_type); void GetParameterValue(double* pvalue,double* gauss,int enum_type); + void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type); void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue); void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type); Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5639) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5640) @@ -2158,9 +2158,9 @@ switch(analysis_type){ - case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum: + case DiagnosticHorizAnalysisEnum: /*default vx,vy and vz: either observation or 0 */ if(!iomodel->vx){ - if (false && iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM) + if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; else for(i=0;i<6;i++)nodeinputs[i]=0; this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); @@ -2169,5 +2169,5 @@ } if(!iomodel->vy){ - if (false && iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM) + if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; else for(i=0;i<6;i++)nodeinputs[i]=0; this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); @@ -4579,6 +4579,4 @@ double L[numdof]; double l1l2l3[3]; - double alpha2_list[3]; - double basalfriction_list[3]={0.0}; double basalfriction; double epsilon[6]; Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5639) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5640) @@ -688,20 +688,24 @@ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - /*Add Tikhonov regularization term to misfit*/ + /*Add Tikhonov regularization term to misfit: + * + * WARNING: for now, the regularization is only taken into account by the gradient + * the misfit reflects the response only! + * + * */ if (control_type==DragCoefficientEnum){ - if (!shelf) return 0; //only shelf?... TO BE CHECKED ISSMASSERT(drag_input); drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; + //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; } else if (control_type==RheologyBbarEnum){ ISSMASSERT(Bbar_input); Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); - Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; + //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; } else if (control_type==DhDtEnum){ ISSMASSERT(dhdt_input); dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); - Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; + //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; } else{ @@ -1029,6 +1033,4 @@ /*nodal functions: */ double l1l2l3[3]; - double alpha2complement_list[numvertices]; //TO BE DELETED - double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED /* strain rate: */ @@ -1068,9 +1070,4 @@ friction=new Friction("2d",inputs,matpar,analysis_type); - /*COMPUT alpha2_list (TO BE DELETED)*/ - for(i=0;iGetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum); - } - /*Retrieve all inputs we will be needing: */ adjointx_input=inputs->GetInput(AdjointxEnum); @@ -1087,7 +1084,6 @@ /*Build alpha_complement_list: */ - //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED - //else alpha_complement=0; - TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED + if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); + else alpha_complement=0; /*Recover alpha_complement and k: */ @@ -2594,9 +2590,9 @@ switch(analysis_type){ - case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum: + case DiagnosticHorizAnalysisEnum: /*default vx,vy and vz: either observation or 0 */ if(!iomodel->vx){ - if (false && iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM) + if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; else for(i=0;i<3;i++)nodeinputs[i]=0; this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); @@ -2605,5 +2601,5 @@ } if(!iomodel->vy){ - if (false && iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM) + if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; else for(i=0;i<3;i++)nodeinputs[i]=0; this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); @@ -3256,6 +3252,4 @@ Friction *friction = NULL; double alpha2; - double alpha2_list[numvertices]; //TO BE DELETED - double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3295,9 +3289,4 @@ friction=new Friction("2d",inputs,matpar,analysis_type); - /*COMPUT alpha2_list (TO BE DELETED)*/ - for(i=0;iGetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); - } - /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); @@ -3307,6 +3296,5 @@ /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED + friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, @@ -3388,6 +3376,4 @@ Friction *friction = NULL; double alpha2; - double alpha2_list[numvertices]; //TO BE DELETED - double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3426,9 +3412,4 @@ friction=new Friction("2d",inputs,matpar,analysis_type); - /*COMPUT alpha2_list (TO BE DELETED)*/ - for(i=0;iGetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); - } - /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); @@ -3438,6 +3419,5 @@ /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED + friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, @@ -3517,6 +3497,4 @@ Friction *friction = NULL; double alpha2; - double alpha2_list[numvertices]; //TO BE DELETED - double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3555,9 +3533,4 @@ friction=new Friction("2d",inputs,matpar,analysis_type); - /*COMPUT alpha2_list (TO BE DELETED)*/ - for(i=0;iGetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); - } - /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); @@ -3567,6 +3540,5 @@ /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED + friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, @@ -5607,9 +5579,4 @@ double alpha2,vx,vy; double geothermalflux_value; - double alpha2_list[numvertices]; //TO BE DELETED - double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED - double vx_list[numvertices]; //TO BE DELETED - double vy_list[numvertices]; //TO BE DELETED - double basalfriction_list[numvertices]; //TO BE DELETED /* gaussian points: */ @@ -5629,7 +5596,7 @@ /*inputs: */ - Input* vx_input=NULL; //TO BE DELETED - Input* vy_input=NULL; //TO BE DELETED - Input* vz_input=NULL; //TO BE DELETED + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; Input* geothermalflux_input=NULL; @@ -5656,14 +5623,4 @@ geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); - /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/ - for(i=0;iGetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED - } - vx_input->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3); //TO BE DELETED - vy_input->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3); //TO BE DELETED - for(i=0;iGetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum); - TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss); // TO BE DELETED + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); /*Calculate scalar parameter*/ Index: /issm/trunk/src/c/objects/Loads/Friction.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction.cpp (revision 5639) +++ /issm/trunk/src/c/objects/Loads/Friction.cpp (revision 5640) @@ -116,4 +116,65 @@ } /*}}}*/ +/*FUNCTION Friction::GetAlpha2{{{1*/ +void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ + + /*This routine calculates the basal friction coefficient + alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ + + /*diverse: */ + double r,s; + double drag_p, drag_q; + double gravity,rho_ice,rho_water; + double Neff; + double thickness,bed; + double vx,vy,vz,vmag; + double drag_coefficient; + double alpha2; + + + /*Recover parameters: */ + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); + inputs->GetParameterValue(&bed, gauss,BedEnum); + inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); + + /*Get material parameters: */ + gravity=matpar->GetG(); + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + //compute r and q coefficients: */ + r=drag_q/drag_p; + s=1./drag_p; + + //From bed and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*bed); + + /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because + the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour + for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should + replace it by Neff=0 (ie, equival it to an ice shelf)*/ + if (Neff<0)Neff=0; + + if(strcmp(element_type,"2d")==0){ + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)); + } + else if (strcmp(element_type,"3d")==0){ + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + inputs->GetParameterValue(&vz, gauss,vzenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); + } + else ISSMERROR("element_type %s not supported yet",element_type); + + alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +} +/*}}}*/ /*FUNCTION Friction::GetAlphaComplement {{{1*/ void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){ @@ -174,2 +235,60 @@ } /*}}}*/ +/*FUNCTION Friction::GetAlphaComplement {{{1*/ +void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){ + + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. + * FrictionGetAlphaComplement is used in control methods on drag, and it computes: + * alpha_complement= Neff ^r * vel ^s*/ + + /*diverse: */ + int i; + double Neff; + double r,s; + double vx; + double vy; + double vmag; + double drag_p,drag_q; + double drag_coefficient; + double bed,thickness; + double gravity,rho_ice,rho_water; + double alpha_complement; + + /*Recover parameters: */ + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); + inputs->GetParameterValue(&bed, gauss,BedEnum); + inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); + + /*Get material parameters: */ + gravity=matpar->GetG(); + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + + //compute r and q coefficients: */ + r=drag_q/drag_p; + s=1./drag_p; + + //From bed and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*bed); + + /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because + the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour + for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should + replace it by Neff=0 (ie, equival it to an ice shelf)*/ + if (Neff<0)Neff=0; + + //We need the velocity magnitude to evaluate the basal stress: + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)); + + alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Friction.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction.h (revision 5639) +++ /issm/trunk/src/c/objects/Loads/Friction.h (revision 5640) @@ -28,5 +28,7 @@ void Echo(void); void GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum); + void GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); void GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum); + void GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum); };