Changeset 25201


Ignore:
Timestamp:
07/03/20 06:23:13 (5 years ago)
Author:
Cheng Gong
Message:

NEW: add an inversion of the friciton coefficient in Schoof's friction law.
In order to reuse the default code of the Budd's law, the coefficient C in Schoof's law is changed to C=sqrt(C_old).
The name of the control_parameter in Schoof's law is called 'FrictionC'.

Location:
issm/trunk-jpl
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

    r24933 r25201  
    11031103        if(control_type!=MaterialsRheologyBbarEnum &&
    11041104                control_type!=FrictionCoefficientEnum   &&
     1105                control_type!=FrictionCEnum   &&
    11051106                control_type!=FrictionAsEnum   &&
    11061107                control_type!=DamageDbarEnum            &&
    … …  
    11261127        switch(control_type){
    11271128                case FrictionCoefficientEnum:
     1129                case FrictionCEnum:
    11281130                        switch(approximation){
    11291131                                case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
    … …  
    16521654        /*Intermediaries*/
    16531655        int      domaintype,dim;
     1656        int frictionlaw;
    16541657        Element* basalelement;
    16551658
    … …  
    16901693        basalelement->GetVerticesCoordinates(&xyz_list);
    16911694        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1692         Input2* dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum);                _assert_(dragcoefficient_input);
     1695
     1696        /* get the friction law: if 11-Schoof, use a special name for the coefficient*/
     1697        element->FindParam(&frictionlaw, FrictionLawEnum);
     1698        Input2* dragcoefficient_input;
     1699        switch(frictionlaw) {
     1700                case 11:
     1701                        dragcoefficient_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoefficient_input);
     1702                        break;
     1703                default:
     1704                        dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
     1705        }
     1706
    16931707        DatasetInput2* weights_input         = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    16941708
    … …  
    19391953        /*Fetch number of vertices for this finite element*/
    19401954        int numvertices = basalelement->GetNumberOfVertices();
     1955        int frictionlaw;
    19411956
    19421957        /*Initialize some vectors*/
    … …  
    19551970        Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    19561971        Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
    1957         Input2* dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1972
     1973        /* get the friction law: 1- Budd, 11-Schoof*/
     1974        element->FindParam(&frictionlaw, FrictionLawEnum);
     1975        Input2* dragcoeff_input;
     1976        switch(frictionlaw) {
     1977                case 1:
     1978                        dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input);
     1979                        break;
     1980                case 11:
     1981                        dragcoeff_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoeff_input);
     1982                        break;
     1983                default:
     1984                        _error_("Friction law "<< frictionlaw <<" not supported in the inversion.");
     1985        }
    19581986
    19591987        /* Start  looping on the number of gaussian points: */

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r24669 r25201  
    5252                        GetAlphaTempComplement(palpha_complement,gauss);
    5353                        break;
     54                case 11:
     55                        GetAlphaSchoofComplement(palpha_complement,gauss);
     56                        break;
    5457          default:
    5558                        _error_("not supported");
    … …  
    162165}
    163166/*}}}*/
     167void Friction::GetAlphaSchoofComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
     168
     169        /* Compute the complement of Schoof's law for inversion
     170         * d alpha2                       
     171         * -------- = |u_b|^(m-1) *(1+ (C/(Cmax N))^(1/m)|u_b|)^(-m-1)
     172         *  dC                           
     173        */
     174        /*diverse: */
     175        IssmDouble  m,Cmax;
     176        IssmDouble  C, coeff;
     177        IssmDouble  alpha_complement;
     178
     179        /*Recover parameters: */
     180        element->GetInputValue(&m,gauss,FrictionMEnum);
     181        element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
     182        element->GetInputValue(&coeff, gauss,FrictionCEnum);
     183
     184        C = coeff*coeff;
     185        /*Get effective pressure*/
     186        IssmDouble Neff = EffectivePressure(gauss);
     187        IssmDouble vmag = VelMag(gauss);
     188
     189        /*Check to prevent dividing by zero if vmag==0*/
     190        if((vmag==0.) || (Neff == 0.)) {
     191                alpha_complement=0.;
     192        }
     193        else {
     194                alpha_complement= pow(vmag, m-1.)*pow((1 + pow(C/(Cmax*Neff),1./m)*vmag), -m-1.);
     195        }
     196        /*Assign output pointers:*/
     197        *palpha_complement=alpha_complement;
     198}/*}}}*/
    164199void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    165200
    … …  
    619654
    620655        /*diverse: */
    621         IssmDouble  C,Cmax,m,alpha2;
     656        IssmDouble  C,coeff,Cmax,m,alpha2;
    622657
    623658        /*Recover parameters: */
    624659        element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
    625         element->GetInputValue(&C,gauss,FrictionCEnum);
     660        element->GetInputValue(&coeff,gauss,FrictionCEnum);
    626661        element->GetInputValue(&m,gauss,FrictionMEnum);
     662
     663        /* scale C for a better inversion */
     664        C = coeff*coeff;
    627665
    628666        /*Get effective pressure and velocity magnitude*/

  • issm/trunk-jpl/src/c/classes/Loads/Friction.h

    r23959 r25201  
    2828                void  GetAlphaTempComplement(IssmDouble* alpha_complement,Gauss* gauss);
    2929                void  GetAlphaViscousComplement(IssmDouble* alpha_complement,Gauss* gauss);
     30                void  GetAlphaSchoofComplement(IssmDouble* alpha_complement,Gauss* gauss);
    3031                void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
    3132                void  GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);

  • issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp

    r24335 r25201  
    3333
    3434        int         domaintype,numcomponents;
     35        int frictionlaw;
    3536        IssmDouble  Jelem=0.;
    3637        IssmDouble  misfit,Jdet;
    … …  
    6162        /*Retrieve all inputs we will be needing: */
    6263        DatasetInput2* weights_input=basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    63         Input2* drag_input   =basalelement->GetInput2(FrictionCoefficientEnum); _assert_(drag_input);
     64
     65        /* get the friction law: if 11-Schoof, which has a different name of C */
     66        element->FindParam(&frictionlaw, FrictionLawEnum);
     67        Input2* drag_input;
     68        switch(frictionlaw) {
     69                case 11:
     70                        drag_input = basalelement->GetInput2(FrictionCEnum); _assert_(drag_input);
     71                        break;
     72                default:
     73                        drag_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(drag_input);
     74        }
    6475
    6576        /* Start  looping on the number of gaussian points: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r24335 r25201  
    122122                        case ThicknessEnum:                           iomodel->FetchData(&independent,&M,&N,"md.geometry.thickness");                              break;
    123123                        case FrictionCoefficientEnum:                 iomodel->FetchData(&independent,&M,&N,"md.friction.coefficient");                            break;
     124                        case FrictionCEnum:                                                                      iomodel->FetchData(&independent,&M,&N,"md.friction.C");                                                            break;
    124125                        case FrictionAsEnum:                          iomodel->FetchData(&independent,&M,&N,"md.friction.As");                                     break;
    125126                        case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.apparent_massbalance");           break;
    … …  
    160161                        case ThicknessEnum:                           iomodel->DeleteData(1,"md.geometry.thickness"); break;
    161162                        case FrictionCoefficientEnum:                 iomodel->DeleteData(1,"md.friction.coefficient"); break;
     163                        case FrictionCEnum:                                        iomodel->DeleteData(1,"md.friction.C"); break;
    162164                        case FrictionAsEnum:                          iomodel->DeleteData(1,"md.friction.As"); break;
    163165                        case BalancethicknessApparentMassbalanceEnum: iomodel->DeleteData(1,"md.balancethickness.apparent_massbalance"); break;

  • issm/trunk-jpl/src/m/inversions/supportedcontrols.m

    r23108 r25201  
    44                'BalancethicknessThickeningRate',...
    55                'FrictionCoefficient',...
     6                'FrictionC',...
    67                'FrictionAs',...
    78                'MaterialsRheologyBbar',...

  • issm/trunk-jpl/test/NightlyRun/test481.m

    r24671 r25201  
    77md.transient.isthermal = 0;
    88md.friction=frictionschoof(md.friction);
    9 md.friction.C    = 20.e4*ones(md.mesh.numberofvertices,1);
     9md.friction.C    = (20.e4)^0.5*ones(md.mesh.numberofvertices,1);
    1010md.friction.Cmax = 0.5*ones(md.mesh.numberofvertices,1);
    1111md.friction.m    = 1./3.*ones(md.mesh.numberofelements,1);
Note: See TracChangeset for help on using the changeset viewer.