Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 27475) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 27476) @@ -2391,4 +2391,5 @@ case 2: case 11: + case 14: dragcoefficient_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoefficient_input); break; @@ -2484,4 +2485,5 @@ case 2: case 11: + case 14: dragcoeff_input = element->GetInput(FrictionCEnum); _assert_(dragcoeff_input); break; @@ -2593,4 +2595,5 @@ case 2: case 11: + case 14: dragcoeff_input = element->GetInput(FrictionCEnum); _assert_(dragcoeff_input); break; @@ -2714,4 +2717,5 @@ case 2: case 11: + case 14: dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input); break; @@ -2832,4 +2836,5 @@ case 2: case 11: + case 14: dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input); break; Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 27475) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 27476) @@ -125,4 +125,7 @@ GetAlphaSchoofComplement(palpha_complement,gauss); break; + case 14: + GetAlphaRegCoulombComplement(palpha_complement,gauss); + break; default: _error_("not supported"); @@ -295,4 +298,35 @@ *palpha_complement=alpha_complement; }/*}}}*/ +void Friction::GetAlphaRegCoulombComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ + + /* Compute the complement of Schoof's law for inversion + * d alpha2 + * -------- = |u_b|^(1/m-1) * (|u_b + u_0|)^(-1/m) + * dC + */ + + /*diverse: */ + IssmDouble m, u0; + IssmDouble alpha_complement; + + /*Recover parameters: */ + element->GetInputValue(&m,gauss,FrictionMEnum); + element->parameters->FindParam(&u0,FrictionU0Enum); + + /*Get velocity magnitude*/ + IssmDouble ub = VelMag(gauss); + + /*Check to prevent dividing by zero if vmag==0*/ + if(ub==0.) { + alpha_complement=0.; + } + else { + /*Compute friction complement*/ + alpha_complement= (pow(ub,1./m-1.)) / pow(ub/u0 + 1.,1./m); + } + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; +}/*}}}*/ void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ @@ -341,4 +375,7 @@ case 13: GetAlpha2Coulomb2(palpha2,gauss); + break; + case 14: + GetAlpha2RegCoulomb(palpha2,gauss); break; default: @@ -910,5 +947,40 @@ *palpha2=alpha2; }/*}}}*/ - +void Friction::GetAlpha2RegCoulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /*This routine calculates the basal friction coefficient + * + * C |u_b|^(1/m-1) + * alpha2= __________________________ + * (|u_b|/u0 + 1 )^(1/m) + * + * */ + + /*diverse: */ + IssmDouble C,coeff,u0,m,alpha2; + + /*Recover parameters: */ + element->GetInputValue(&coeff,gauss,FrictionCEnum); + element->GetInputValue(&m,gauss,FrictionMEnum); + element->parameters->FindParam(&u0,FrictionU0Enum); + + /* scale C for a better inversion */ + C = coeff*coeff; + + /*Get velocity magnitude*/ + IssmDouble ub = VelMag(gauss); + + /*Check to prevent dividing by zero if vmag==0*/ + if(ub==0.) { + alpha2=0.; + } + else { + /*Compute alpha^2*/ + alpha2= (C*pow(ub,1./m-1.)) / pow(ub/u0 + 1.,1./m); + } + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/ /*Get effective pressure as a function of flag */ @@ -1192,4 +1264,8 @@ } break; + case 14: + iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum); + iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum); + break; default: _error_("friction law "<< frictionlaw <<" not supported"); @@ -1264,4 +1340,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum)); break; + case 14: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.u0",FrictionU0Enum)); + break; default: _error_("Friction law "<GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - /* get the friction law: if 11-Schoof, which has a different name of C */ + /* get the friction law: if 2-Weertman, 11-Schoof or 14-RegularizedCoulomb, which has a different names of C */ element->FindParam(&frictionlaw, FrictionLawEnum); Input* drag_input; @@ -69,4 +69,5 @@ case 2: case 11: + case 14: drag_input = basalelement->GetInput(FrictionCEnum); _assert_(drag_input); break; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27475) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27476) @@ -205,4 +205,5 @@ FrictionLinearizeEnum, FrictionPseudoplasticityExponentEnum, + FrictionU0Enum, FrictionThresholdSpeedEnum, FrictionVoidRatioEnum, Index: /issm/trunk-jpl/test/NightlyRun/test358.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test358.m (revision 27476) +++ /issm/trunk-jpl/test/NightlyRun/test358.m (revision 27476) @@ -0,0 +1,43 @@ +%Test Name: SquareSheetConstrainedCMDragRegCoulombSSA2d +md=triangle(model(),'../Exp/Square.exp',200000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md=setflowequation(md,'SSA','all'); + +%use Regularized Coulomb's law +md.friction = frictionregcoulomb(); +md.friction.m = 3.0*ones(md.mesh.numberofelements,1); +md.friction.u0 = 2000; %m/yr +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.verbose = verbose('all'); +md.debug.valgrind = 0; + +md.cluster=generic('name',oshostname(),'np',1); +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) +};