Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 25200) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 25201) @@ -1103,4 +1103,5 @@ if(control_type!=MaterialsRheologyBbarEnum && control_type!=FrictionCoefficientEnum && + control_type!=FrictionCEnum && control_type!=FrictionAsEnum && control_type!=DamageDbarEnum && @@ -1126,4 +1127,5 @@ switch(control_type){ case FrictionCoefficientEnum: + case FrictionCEnum: switch(approximation){ case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break; @@ -1652,4 +1654,5 @@ /*Intermediaries*/ int domaintype,dim; + int frictionlaw; Element* basalelement; @@ -1690,5 +1693,16 @@ basalelement->GetVerticesCoordinates(&xyz_list); basalelement->GradientIndexing(&vertexpidlist[0],control_index); - Input2* dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoefficient_input); + + /* get the friction law: if 11-Schoof, use a special name for the coefficient*/ + element->FindParam(&frictionlaw, FrictionLawEnum); + Input2* dragcoefficient_input; + switch(frictionlaw) { + case 11: + dragcoefficient_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoefficient_input); + break; + default: + dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoefficient_input); + } + DatasetInput2* weights_input = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); @@ -1939,4 +1953,5 @@ /*Fetch number of vertices for this finite element*/ int numvertices = basalelement->GetNumberOfVertices(); + int frictionlaw; /*Initialize some vectors*/ @@ -1955,5 +1970,18 @@ Input2* adjointx_input = basalelement->GetInput2(AdjointxEnum); _assert_(adjointx_input); Input2* adjointy_input = basalelement->GetInput2(AdjointyEnum); _assert_(adjointy_input); - Input2* dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input); + + /* get the friction law: 1- Budd, 11-Schoof*/ + element->FindParam(&frictionlaw, FrictionLawEnum); + Input2* dragcoeff_input; + switch(frictionlaw) { + case 1: + dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input); + break; + case 11: + dragcoeff_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoeff_input); + break; + default: + _error_("Friction law "<< frictionlaw <<" not supported in the inversion."); + } /* Start looping on the number of gaussian points: */ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 25200) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 25201) @@ -52,4 +52,7 @@ GetAlphaTempComplement(palpha_complement,gauss); break; + case 11: + GetAlphaSchoofComplement(palpha_complement,gauss); + break; default: _error_("not supported"); @@ -162,4 +165,36 @@ } /*}}}*/ +void Friction::GetAlphaSchoofComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ + + /* Compute the complement of Schoof's law for inversion + * d alpha2 + * -------- = |u_b|^(m-1) *(1+ (C/(Cmax N))^(1/m)|u_b|)^(-m-1) + * dC + */ + /*diverse: */ + IssmDouble m,Cmax; + IssmDouble C, coeff; + IssmDouble alpha_complement; + + /*Recover parameters: */ + element->GetInputValue(&m,gauss,FrictionMEnum); + element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum); + element->GetInputValue(&coeff, gauss,FrictionCEnum); + + C = coeff*coeff; + /*Get effective pressure*/ + IssmDouble Neff = EffectivePressure(gauss); + IssmDouble vmag = VelMag(gauss); + + /*Check to prevent dividing by zero if vmag==0*/ + if((vmag==0.) || (Neff == 0.)) { + alpha_complement=0.; + } + else { + alpha_complement= pow(vmag, m-1.)*pow((1 + pow(C/(Cmax*Neff),1./m)*vmag), -m-1.); + } + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; +}/*}}}*/ void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ @@ -619,10 +654,13 @@ /*diverse: */ - IssmDouble C,Cmax,m,alpha2; + IssmDouble C,coeff,Cmax,m,alpha2; /*Recover parameters: */ element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum); - element->GetInputValue(&C,gauss,FrictionCEnum); + element->GetInputValue(&coeff,gauss,FrictionCEnum); element->GetInputValue(&m,gauss,FrictionMEnum); + + /* scale C for a better inversion */ + C = coeff*coeff; /*Get effective pressure and velocity magnitude*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 25200) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 25201) @@ -28,4 +28,5 @@ void GetAlphaTempComplement(IssmDouble* alpha_complement,Gauss* gauss); void GetAlphaViscousComplement(IssmDouble* alpha_complement,Gauss* gauss); + void GetAlphaSchoofComplement(IssmDouble* alpha_complement,Gauss* gauss); void GetAlpha2(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss); Index: /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 25200) +++ /issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (revision 25201) @@ -33,4 +33,5 @@ int domaintype,numcomponents; + int frictionlaw; IssmDouble Jelem=0.; IssmDouble misfit,Jdet; @@ -61,5 +62,15 @@ /*Retrieve all inputs we will be needing: */ DatasetInput2* weights_input=basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input2* drag_input =basalelement->GetInput2(FrictionCoefficientEnum); _assert_(drag_input); + + /* get the friction law: if 11-Schoof, which has a different name of C */ + element->FindParam(&frictionlaw, FrictionLawEnum); + Input2* drag_input; + switch(frictionlaw) { + case 11: + drag_input = basalelement->GetInput2(FrictionCEnum); _assert_(drag_input); + break; + default: + drag_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(drag_input); + } /* Start looping on the number of gaussian points: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 25200) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 25201) @@ -122,4 +122,5 @@ case ThicknessEnum: iomodel->FetchData(&independent,&M,&N,"md.geometry.thickness"); break; case FrictionCoefficientEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.coefficient"); break; + case FrictionCEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.C"); break; case FrictionAsEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.As"); break; case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.apparent_massbalance"); break; @@ -160,4 +161,5 @@ case ThicknessEnum: iomodel->DeleteData(1,"md.geometry.thickness"); break; case FrictionCoefficientEnum: iomodel->DeleteData(1,"md.friction.coefficient"); break; + case FrictionCEnum: iomodel->DeleteData(1,"md.friction.C"); break; case FrictionAsEnum: iomodel->DeleteData(1,"md.friction.As"); break; case BalancethicknessApparentMassbalanceEnum: iomodel->DeleteData(1,"md.balancethickness.apparent_massbalance"); break; Index: /issm/trunk-jpl/src/m/inversions/supportedcontrols.m =================================================================== --- /issm/trunk-jpl/src/m/inversions/supportedcontrols.m (revision 25200) +++ /issm/trunk-jpl/src/m/inversions/supportedcontrols.m (revision 25201) @@ -4,4 +4,5 @@ 'BalancethicknessThickeningRate',... 'FrictionCoefficient',... + 'FrictionC',... 'FrictionAs',... 'MaterialsRheologyBbar',... Index: /issm/trunk-jpl/test/NightlyRun/test356.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test356.m (revision 25201) +++ /issm/trunk-jpl/test/NightlyRun/test356.m (revision 25201) @@ -0,0 +1,41 @@ +%Test Name: SquareSheetConstrainedCMDragSchoofSSA2d +md=triangle(model(),'../Exp/Square.exp',200000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md=setflowequation(md,'SSA','all'); + +%use Schoof's law +Cmax = 0.8; +md.friction = frictionschoof(); +md.friction.m = 1.0/3.0*ones(md.mesh.numberofelements,1); +md.friction.Cmax = Cmax*ones(md.mesh.numberofvertices,1); +md.friction.C = 200*ones(md.mesh.numberofvertices,1); + +%control parameters +md.inversion.iscontrol=1; +md.inversion.control_parameters={'FrictionC'}; +md.inversion.min_parameters=1.*ones(md.mesh.numberofvertices,1); +md.inversion.max_parameters=10000.*ones(md.mesh.numberofvertices,1); +md.inversion.nsteps=2; +md.inversion.cost_functions=[102 501]; +md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2); md.inversion.cost_functions_coefficients(:,2)=2.*10^-7; +md.inversion.gradient_scaling=3.*ones(md.inversion.nsteps,1); +md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1); +md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1); +md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Stressbalance'); + +%Fields and tolerances to track changes +field_names ={'Gradient','Misfits','FrictionC','Pressure','Vel','Vx','Vy'}; +field_tolerances={1e-12,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.StressbalanceSolution.Gradient1),... + (md.results.StressbalanceSolution.J),... + (md.results.StressbalanceSolution.FrictionC),... + (md.results.StressbalanceSolution.Pressure),... + (md.results.StressbalanceSolution.Vel),... + (md.results.StressbalanceSolution.Vx),... + (md.results.StressbalanceSolution.Vy) +}; Index: /issm/trunk-jpl/test/NightlyRun/test481.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test481.m (revision 25200) +++ /issm/trunk-jpl/test/NightlyRun/test481.m (revision 25201) @@ -7,5 +7,5 @@ md.transient.isthermal = 0; md.friction=frictionschoof(md.friction); -md.friction.C = 20.e4*ones(md.mesh.numberofvertices,1); +md.friction.C = (20.e4)^0.5*ones(md.mesh.numberofvertices,1); md.friction.Cmax = 0.5*ones(md.mesh.numberofvertices,1); md.friction.m = 1./3.*ones(md.mesh.numberofelements,1);