Changeset 26457


Ignore:
Timestamp:
09/23/21 08:46:01 (4 years ago)
Author:
Cheng Gong
Message:

CHG: add VxSurface, VxBase, VySurface, VyBase to all 2D models: SSA, L1L2 and MLHO. Modified friction and misfit functions to use these variables instead of Vx and Vy.

Location:
issm/trunk-jpl/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r26442 r26457  
    792792                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.);
    793793                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.);
     794        }
     795   if(iomodel->domaintype==Domain2DhorizontalEnum){
    794796                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.);
    795797                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.);
     798                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxSurfaceEnum,0.);
     799                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VySurfaceEnum,0.);
    796800        }
    797801        iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum);
    … …  
    21202124
    21212125        /*Add vx and vy as inputs to the tria element: */
     2126        /*Also add surface vx and vy for the misfit, and base vx and vy for friction*/
    21222127        element->AddBasalInput(VxEnum,vx,element->GetElementType());
    2123         if(dim==2)element->AddBasalInput(VyEnum,vy,element->GetElementType());
     2128        element->AddBasalInput(VxSurfaceEnum,vx,element->GetElementType());
     2129        element->AddBasalInput(VxBaseEnum,vx,element->GetElementType());
     2130        if(dim==2) {
     2131                element->AddBasalInput(VyEnum,vy,element->GetElementType());
     2132                element->AddBasalInput(VySurfaceEnum,vy,element->GetElementType());
     2133                element->AddBasalInput(VyBaseEnum,vy,element->GetElementType());
     2134        }
    21242135        element->AddBasalInput(VelEnum,vel,element->GetElementType());
    2125 
     2136       
    21262137        /*Free ressources:*/
    21272138        xDelete<IssmDouble>(vel);
    … …  
    27242735        element->AddBasalInput(VelEnum,vel,element->GetElementType());
    27252736
     2737        /*Also add surface vx and vy for the misfit, and base vx and vy for friction*/
     2738        element->AddBasalInput(VxSurfaceEnum,vx,element->GetElementType());
     2739        element->AddBasalInput(VySurfaceEnum,vy,element->GetElementType());
     2740        element->AddBasalInput(VxBaseEnum,vx,element->GetElementType());
     2741        element->AddBasalInput(VyBaseEnum,vy,element->GetElementType());
     2742
    27262743        /*Free ressources:*/
    27272744        xDelete<IssmDouble>(vel);
    … …  
    27892806
    27902807        /*build friction object, used later on: */
    2791         Friction* friction=new Friction(element,dim,true);
     2808        Friction* friction=new Friction(element, dim);
    27922809
    27932810        /*Recover portion of element that is grounded*/
    … …  
    32803297                vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2);
    32813298                vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
    3282                 // HOTFIX: set surface velocity to VxEnum and VyEnum
    3283         //      vx[i]=vbx[i]+vshx[i];
    3284         //      vy[i]=vby[i]+vshy[i];
    32853299        }
    32863300               
    … …  
    34743488
    34753489        /*build friction object, used later on: */
    3476         Friction* friction=new Friction(element,dim==3?2:1);
     3490        /*dim=4 is special for HO, which is actually 2.5D*/
     3491        Friction* friction=new Friction(element,dim==3?4:1);
    34773492
    34783493        /*Recover portion of element that is grounded*/

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

    r26206 r26457  
    2020        this->dim=0;
    2121        this->law=0;
    22         this->isMLHO=false;
    23 
    2422}
    2523/*}}}*/
    26 Friction::Friction(Element* element_in,int dim_in,bool isMLHO_in){/*{{{*/
     24Friction::Friction(Element* element_in,int dim_in){/*{{{*/
    2725
    2826        this->element=element_in;
    2927        this->dim=dim_in;
    30         this->isMLHO=isMLHO_in;
    3128        element_in->FindParam(&this->law,FrictionLawEnum);
    3229}
    … …  
    859856        IssmDouble vx,vy,vz,vmag;
    860857       
    861         if(this->isMLHO) _assert_(dim==2);
    862 
    863858        switch(dim){
    864859                case 1:
    … …  
    867862                        break;
    868863                case 2:
    869                         if(this->isMLHO){
    870                                 element->GetInputValue(&vx,gauss,VxBaseEnum);
    871                                 element->GetInputValue(&vy,gauss,VyBaseEnum);
    872                         }
    873                         else{
    874                                 element->GetInputValue(&vx,gauss,VxEnum);
    875                                 element->GetInputValue(&vy,gauss,VyEnum);
    876                         }               
     864                        element->GetInputValue(&vx,gauss,VxBaseEnum);
     865                        element->GetInputValue(&vy,gauss,VyBaseEnum);
    877866                        vmag=sqrt(vx*vx+vy*vy);
    878867                        break;
    … …  
    883872                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
    884873                        break;
     874                case 4:
     875                        /* This is for HO 3D case, since it requires the horizontal velocities */
     876                        element->GetInputValue(&vx,gauss,VxEnum);
     877                        element->GetInputValue(&vy,gauss,VyEnum);
     878                        vmag=sqrt(vx*vx+vy*vy);
     879                        break;
    885880                default:
    886881                        _error_("not supported");

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

    r26206 r26457  
    1717                int      dim;
    1818                int      law;
    19                 bool            isMLHO=false;
    2019
    2120                /*methods: */
    2221                Friction();
    23                 Friction(Element* element_in,int dim_in,bool isMLHO=0);
     22                Friction(Element* element_in,int dim_in);
    2423                ~Friction();
    2524

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

    r26090 r26457  
    6161        /*Retrieve all inputs we will be needing: */
    6262        DatasetInput* weights_input=topelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    63         Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    64         Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     63        Input* vx_input     = NULL;
     64        Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
    6565        Input* vy_input     = NULL;
    6666        Input* vyobs_input  = NULL;
    67         if(numcomponents==2){
    68                 vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
    69                 vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     67
     68        /*Read SurfaceEnum from 2D models:SSA, L1L2, MLHO*/
     69        if (domaintype == Domain2DhorizontalEnum) {
     70                vx_input = topelement->GetInput(VxSurfaceEnum);                                 _assert_(vx_input);
     71                if(numcomponents==2){
     72                        vy_input    =topelement->GetInput(VySurfaceEnum);        _assert_(vy_input);
     73                        vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     74                }
     75        }
     76        else {
     77                vx_input = topelement->GetInput(VxEnum);                                                        _assert_(vx_input);
     78                if(numcomponents==2){
     79                        vy_input    =topelement->GetInput(VyEnum);                                      _assert_(vy_input);
     80                        vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     81                }
    7082        }
    7183

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

    r26090 r26457  
    6969        DatasetInput* weights_input=topelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    7070        Input* S_input      = topelement->GetInput(SurfaceAreaEnum);     _assert_(S_input);
    71         Input* vx_input     = topelement->GetInput(VxEnum);              _assert_(vx_input);
    72         Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum);  _assert_(vxobs_input);
    73         Input* vy_input     = NULL;
    74         Input* vyobs_input  = NULL;
    75         if(numcomponents==2){
    76                 vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
    77                 vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     71   Input* vx_input     = NULL;
     72   Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
     73   Input* vy_input     = NULL;
     74   Input* vyobs_input  = NULL;
     75
     76   /*Read SurfaceEnum from 2D models:SSA, L1L2, MLHO*/
     77   if (domaintype == Domain2DhorizontalEnum) {
     78      vx_input = topelement->GetInput(VxSurfaceEnum);             _assert_(vx_input);
     79      if(numcomponents==2){
     80         vy_input    =topelement->GetInput(VySurfaceEnum);        _assert_(vy_input);
     81         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     82      }
    7883        }
     84   else {
     85      vx_input = topelement->GetInput(VxEnum);                    _assert_(vx_input);
     86      if(numcomponents==2){
     87         vy_input    =topelement->GetInput(VyEnum);               _assert_(vy_input);
     88         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     89      }
     90   }
    7991
    8092        /* Start  looping on the number of gaussian points: */

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

    r26090 r26457  
    6363        /*Retrieve all inputs we will be needing: */
    6464        DatasetInput* weights_input=topelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    65         Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    66         Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    67         Input* vy_input     = NULL;
    68         Input* vyobs_input  = NULL;
    69         if(numcomponents==2){
    70                 vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
    71                 vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     65   Input* vx_input     = NULL;
     66   Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum);   _assert_(vxobs_input);
     67   Input* vy_input     = NULL;
     68   Input* vyobs_input  = NULL;
     69
     70   /*Read SurfaceEnum from 2D models:SSA, L1L2, MLHO*/
     71   if (domaintype == Domain2DhorizontalEnum) {
     72      vx_input = topelement->GetInput(VxSurfaceEnum);             _assert_(vx_input);
     73      if(numcomponents==2){
     74         vy_input    =topelement->GetInput(VySurfaceEnum);        _assert_(vy_input);
     75         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     76      }
    7277        }
     78   else {
     79      vx_input = topelement->GetInput(VxEnum);                    _assert_(vx_input);
     80      if(numcomponents==2){
     81         vy_input    =topelement->GetInput(VyEnum);               _assert_(vy_input);
     82         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     83      }
     84   }
    7385
    7486        /* Start  looping on the number of gaussian points: */

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

    r26090 r26457  
    6363        /*Retrieve all inputs we will be needing: */
    6464        DatasetInput* weights_input=topelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    65         Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    66         Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    67         Input* vy_input     = NULL;
    68         Input* vyobs_input  = NULL;
    69         if(numcomponents==2){
    70                 vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
    71                 vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     65   Input* vx_input     = NULL;
     66   Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum);   _assert_(vxobs_input);
     67   Input* vy_input     = NULL;
     68   Input* vyobs_input  = NULL;
     69
     70   /*Read SurfaceEnum from 2D models:SSA, L1L2, MLHO*/
     71   if (domaintype == Domain2DhorizontalEnum) {
     72      vx_input = topelement->GetInput(VxSurfaceEnum);             _assert_(vx_input);
     73      if(numcomponents==2){
     74         vy_input    =topelement->GetInput(VySurfaceEnum);        _assert_(vy_input);
     75         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     76      }
    7277        }
     78   else {
     79      vx_input = topelement->GetInput(VxEnum);                    _assert_(vx_input);
     80      if(numcomponents==2){
     81         vy_input    =topelement->GetInput(VyEnum);               _assert_(vy_input);
     82         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     83      }
     84   }
    7385
    7486        /* Start  looping on the number of gaussian points: */

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

    r26090 r26457  
    6363        /*Retrieve all inputs we will be needing: */
    6464        DatasetInput* weights_input=topelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    65         Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    66         Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    67         Input* vy_input     = NULL;
    68         Input* vyobs_input  = NULL;
    69         if(numcomponents==2){
    70                 vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
    71                 vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     65   Input* vx_input     = NULL;
     66   Input* vxobs_input  = topelement->GetInput(InversionVxObsEnum);   _assert_(vxobs_input);
     67   Input* vy_input     = NULL;
     68   Input* vyobs_input  = NULL;
     69
     70   /*Read SurfaceEnum from 2D models:SSA, L1L2, MLHO*/
     71   if (domaintype == Domain2DhorizontalEnum) {
     72      vx_input = topelement->GetInput(VxSurfaceEnum);             _assert_(vx_input);
     73      if(numcomponents==2){
     74         vy_input    =topelement->GetInput(VySurfaceEnum);        _assert_(vy_input);
     75         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     76      }
    7277        }
     78   else {
     79      vx_input = topelement->GetInput(VxEnum);                    _assert_(vx_input);
     80      if(numcomponents==2){
     81         vy_input    =topelement->GetInput(VyEnum);               _assert_(vy_input);
     82         vyobs_input =topelement->GetInput(InversionVyObsEnum);   _assert_(vyobs_input);
     83      }
     84   }
    7385
    7486        /* Start  looping on the number of gaussian points: */

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26445 r26457  
    536536        AdjointpEnum,
    537537        AdjointxEnum,
     538        AdjointxBaseEnum,
     539        AdjointxShearEnum,
    538540        AdjointyEnum,
     541        AdjointyBaseEnum,
     542        AdjointyShearEnum,
    539543        AdjointzEnum,
    540544        AirEnum,
Note: See TracChangeset for help on using the changeset viewer.