Changeset 14967


Ignore:
Timestamp:
05/08/13 13:48:23 (12 years ago)
Author:
bdef
Message:

BUG: Modification in the residual computation

Location:
issm/trunk-jpl/src
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r14960 r14967  
    14591459#ifdef  _HAVE_DAKOTA_
    14601460void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
    1461 
     1461       
    14621462        int        i,j;
    14631463        int        my_rank;
    … …  
    15621562void FemModel::Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
    15631563
    1564         int      i;
    1565 
     1564  int      i;
     1565       
    15661566        /*intermediary: */
    15671567        Element *element     = NULL;
    1568 
     1568       
    15691569        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    15701570        for (i=0;i<elements->Size();i++){
    … …  
    15751575/*}}}*/
    15761576#endif
     1577
     1578void FemModel::BasisIntegralsx(void){ /*{{{*/
     1579
     1580        Vector<IssmDouble>* basisg=NULL;
     1581       
     1582        /*Vector allocation*/
     1583        basisg=new Vector<IssmDouble>(nodes->NumberOfNodes());
     1584       
     1585        for (int i=0;i<elements->Size();i++){
     1586                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1587                element->BasisIntegral(basisg);
     1588        }
     1589        /*Assemble*/
     1590        basisg->Assemble();
     1591       
     1592        /*Update Inputs*/
     1593        InputUpdateFromVectorx(elements,nodes,vertices,loads,materials,parameters,basisg,BasisIntegralEnum,NodesEnum);
     1594       
     1595}
     1596/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r14807 r14967  
    100100                void UpdateConstraintsx(void);
    101101                int  UpdateVertexPositionsx(void);
     102                void BasisIntegralsx(void);
    102103};
    103104

  • issm/trunk-jpl/src/c/classes/objects/Elements/Element.h

    r14953 r14967  
    128128
    129129                #ifdef _HAVE_HYDROLOGY_
    130                 virtual void  GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
     130                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
     131                virtual void BasisIntegral(Vector<IssmDouble>* basisg)=0;
    131132                #endif
    132133};

  • issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

    r14953 r14967  
    93479347}
    93489348/*}}}*/
     9349/*FUNCTION Tria::BasisIntegral {{{*/
     9350void Penta::BasisIntegral(Vector<IssmDouble>* basisg){
     9351        _error_("Hydrological stuff not suported in Penta");
     9352}
     9353/*}}}*/
     9354
    93499355#endif

  • issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h

    r14960 r14967  
    305305                void    CreateHydrologyWaterVelocityInput(void);
    306306                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
     307                void    BasisIntegral(Vector<IssmDouble>* gbasis);
    307308                #endif
    308309                #ifdef _HAVE_THERMAL_

  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

    r14953 r14967  
    18811881void  Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type){
    18821882
     1883       
     1884        const int   numdof         = NDOF1 *NUMVERTICES;
     1885        int        *doflist        = NULL;
     1886        IssmDouble  values[numdof];
     1887
    18831888        /*Check that name is an element input*/
    18841889        if (!IsInput(name)) return;
    18851890
    18861891        switch(type){
    1887 
    1888                 case VertexEnum: {
    1889 
    1890                         /*New TriaP1Input*/
    1891                         IssmDouble values[3];
    1892 
    1893                         /*Get values on the 3 vertices*/
    1894                         for (int i=0;i<3;i++){
    1895                                 values[i]=vector[this->nodes[i]->GetVertexPid()];
    1896                         }
    1897 
    1898                         /*update input*/
    1899                         if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    1900                                 material->inputs->AddInput(new TriaP1Input(name,values));
    1901                         }
    1902                         else{
    1903                                 this->inputs->AddInput(new TriaP1Input(name,values));
    1904                         }
    1905                         return;
     1892        case VertexEnum:
     1893                /*Get values on the 3 vertices*/
     1894                for (int i=0;i<3;i++){
     1895                        values[i]=vector[this->nodes[i]->GetVertexPid()];
    19061896                }
    1907                 default:
     1897                /*update input*/
     1898                if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
     1899                        material->inputs->AddInput(new TriaP1Input(name,values));
     1900                }
     1901                else{
     1902                        this->inputs->AddInput(new TriaP1Input(name,values));
     1903                }
     1904                return;
     1905
     1906        case NodesEnum:
     1907                /*Get dof list: */
     1908                GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1909               
     1910                /*Use the dof list to index into the vector: */
     1911                for(int i=0;i<numdof;i++){
     1912                        values[i]=vector[doflist[i]];
     1913                        if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1914                }
     1915                /*Add input to the element: */
     1916                this->inputs->AddInput(new TriaP1Input(BasisIntegralEnum,values));
     1917               
     1918                /*Free ressources:*/
     1919                xDelete<int>(doflist);
     1920        default:
     1921       
    19081922                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    19091923        }
    … …  
    20532067                                name==OldGradientEnum ||
    20542068        name==ConvergedEnum ||
     2069                                name==BasisIntegralEnum ||
    20552070                                name==QmuVxEnum ||
    20562071                                name==QmuVyEnum ||
    … …  
    61756190                /*Loading term*/
    61766191                water_input->GetInputValue(&water_load,gauss);
     6192                if(this->id==1) printf("Qin %e \n ", water_load);
    61776193                scalar = Jdet*gauss->weight*water_load;
    61786194                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    … …  
    62326248                /*Loading term*/
    62336249                residual_input->GetInputValue(&residual_load,gauss);
     6250                if(this->id==1) printf("Qres %e \n ", residual_load);
    62346251                scalar = Jdet*gauss->weight*residual_load;
    62356252                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    … …  
    63166333        const int   numdof         = NDOF1 *NUMVERTICES;
    63176334        int        *doflist        = NULL;
    6318         bool  converged;
     6335        bool        converged;
    63196336        IssmDouble  values[numdof];
    63206337        IssmDouble  residual[numdof];
     6338        IssmDouble  intbasis[numdof];   
    63216339        IssmDouble  sediment_storing;
    6322         IssmDouble  penalty_factor;
     6340        IssmDouble  penalty_factor, dt;
    63236341        IssmDouble  kmax, kappa, h_max;
    63246342
    … …  
    63356353        /*If converged keep the residual in mind*/
    63366354        this->inputs->GetInputValue(&converged,ConvergedEnum);
     6355        GetInputListOnVertices(&intbasis[0],BasisIntegralEnum);
     6356
     6357        if(this->id==1) printf("val %e \n ", values[1]);
     6358       
    63376359        if(converged){
     6360                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    63386361                this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
    63396362                this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
    … …  
    63446367                        this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    63456368                        if(values[i]>h_max)
    6346                                 residual[i]=kappa*(values[i]-h_max)*sediment_storing;
     6369                                residual[i]=kappa*(values[i]-h_max)*sediment_storing/(dt*intbasis[i]);
    63476370                        else
    63486371                                residual[i]=0.0;
    63496372                }
    6350         }
     6373                if(this->id==1){
     6374                        printf("res  %e val %e h-max %e Stor %e \n ", residual[1], values[1], h_max, sediment_storing);
     6375                        printf("kappa %e kmax %e pen %e dt %e \n", kappa, kmax,penalty_factor,dt);
     6376                }
     6377        }       
    63516378
    63526379        /*Add input to the element: */
    … …  
    64216448}
    64226449/*}}}*/
     6450/*FUNCTION Tria::BasisIntegral {{{*/
     6451void Tria::BasisIntegral(Vector<IssmDouble>* basisg){
     6452
     6453        /*Constants*/
     6454        const int    numdof=NDOF1*NUMVERTICES;
     6455
     6456        /*Intermediaries */
     6457        IssmDouble Jdet;
     6458        IssmDouble xyz_list[NUMVERTICES][3];
     6459        IssmDouble basis[numdof];
     6460        IssmDouble basisint[numdof]={0.};
     6461        int       *doflist=NULL;
     6462        GaussTria* gauss=NULL;
     6463
     6464        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6465
     6466        /* Start looping on the number of gaussian points: */
     6467        gauss=new GaussTria(2);
     6468        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6469
     6470                gauss->GaussPoint(ig);
     6471
     6472                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6473                GetNodalFunctions(&basis[0], gauss);
     6474
     6475                for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
     6476        }
     6477
     6478        basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
     6479
     6480        /*Clean up and return*/
     6481        delete gauss;
     6482        xDelete<int>(doflist);
     6483}
     6484/*}}}*/
     6485
    64236486#endif
    64246487

  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h

    r14960 r14967  
    253253                void      InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
    254254                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
     255                void    BasisIntegral(Vector<IssmDouble>* gbasis);
    255256                #endif
    256257                #ifdef _HAVE_BALANCED_

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

    r14960 r14967  
    108108        HydrologyEfficientEnum,
    109109        HydrologySedimentKmaxEnum,
     110        BasisIntegralEnum,
    110111        IndependentObjectEnum,
    111112        InversionControlParametersEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r14960 r14967  
    116116                case HydrologyEfficientEnum : return "HydrologyEfficient";
    117117                case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax";
     118                case BasisIntegralEnum : return "BasisIntegral";
    118119                case IndependentObjectEnum : return "IndependentObject";
    119120                case InversionControlParametersEnum : return "InversionControlParameters";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r14960 r14967  
    116116              else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
    117117              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
     118              else if (strcmp(name,"BasisIntegral")==0) return BasisIntegralEnum;
    118119              else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
    119120              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    … …  
    136137              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    137138              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    138               else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
     142              if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
     143              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    143144              else if (strcmp(name,"MaskElementonfloatingice")==0) return MaskElementonfloatingiceEnum;
    144145              else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
    … …  
    259260              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    260261              else if (strcmp(name,"TransientIsprognostic")==0) return TransientIsprognosticEnum;
    261               else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
     265              if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum;
     266              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    266267              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
    267268              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    … …  
    382383              else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
    383384              else if (strcmp(name,"Vertex")==0) return VertexEnum;
    384               else if (strcmp(name,"Air")==0) return AirEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"Ice")==0) return IceEnum;
     388              if (strcmp(name,"Air")==0) return AirEnum;
     389              else if (strcmp(name,"Ice")==0) return IceEnum;
    389390              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    390391              else if (strcmp(name,"Water")==0) return WaterEnum;
    … …  
    505506              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    506507              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    507               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
     511              if (strcmp(name,"MinVx")==0) return MinVxEnum;
     512              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    512513              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    513514              else if (strcmp(name,"MinVy")==0) return MinVyEnum;

  • issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp

    r14960 r14967  
    7171                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
    7272                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    73                                 /*InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);*/
     73                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);
    7474                                break;
    7575                        }

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r14963 r14967  
    13951395        return StringToEnum('HydrologySedimentKmax')[0]
    13961396
     1397def BasisIntegralEnum():
     1398        """
     1399        BASISINTEGRALENUM - Enum of BasisIntegral
     1400
     1401        WARNING: DO NOT MODIFY THIS FILE
     1402                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1403                                Please read src/c/shared/Enum/README for more information
     1404
     1405           Usage:
     1406              macro=BasisIntegralEnum()
     1407        """
     1408
     1409        return StringToEnum('BasisIntegral')[0]
     1410
    13971411def IndependentObjectEnum():
    13981412        """
    … …  
    76657679        """
    76667680
    7667         return 546
    7668 
     7681        return 547
     7682

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r14963 r14967  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=546;
     11macro=547;
Note: See TracChangeset for help on using the changeset viewer.