Changeset 15377


Ignore:
Timestamp:
07/01/13 15:34:52 (12 years ago)
Author:
bdef
Message:

CHG:updating Hydro Model to deal with real 3D geometry

Location:
issm/trunk-jpl/src/c/classes
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/ElementResults/TriaP1ElementResult.cpp

    r15128 r15377  
    118118/*FUNCTION TriaP1ElementResult::GetElementVectorFromResults{{{*/
    119119void TriaP1ElementResult::GetElementVectorFromResults(Vector<IssmDouble>* vector,int dof){
    120         _error_("Result " << EnumToStringx(enum_type) << " is a TriaP1ElementResult and should not write vector of size numberofelemenrs");
     120        _error_("Result " << EnumToStringx(enum_type) << " is a TriaP1ElementResult and should not write vector of size numberofelements");
    121121} /*}}}*/

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

    r15375 r15377  
    455455                case MeltingAnalysisEnum:
    456456                        Ke=CreateKMatrixMelting();
     457                        break;
     458                #endif
     459                #ifdef _HAVE_HYDROLOGY_
     460          case HydrologyDCInefficientAnalysisEnum:
     461                        Ke=CreateKMatrixHydrologyDCInefficient();
     462                        break;
     463          case HydrologyDCEfficientAnalysisEnum:
     464                        Ke=CreateKMatrixHydrologyDCEfficient();
    457465                        break;
    458466                #endif
    … …  
    585593                        break;
    586594                #endif
     595                #ifdef _HAVE_HYDROLOGY_
     596          case HydrologyDCInefficientAnalysisEnum:
     597                        pe=CreatePVectorHydrologyDCInefficient();
     598                        break;
     599          case HydrologyDCEfficientAnalysisEnum:
     600                        pe=CreatePVectorHydrologyDCEfficient();
     601                        break;
     602                #endif
    587603                default:
    588604                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    11621178                break;
    11631179        case DiagnosticVertAnalysisEnum:
    1164                 GetSolutionFromInputsDiagnosticVert(solution);
     1180                //GetSolutionFromInputsDiagnosticVert(solution);
     1181                GetSolutionFromInputsOneDof(solution, VzEnum);
    11651182                break;
    11661183        #endif
    11671184        #ifdef _HAVE_THERMAL_
    11681185        case ThermalAnalysisEnum:
    1169                 GetSolutionFromInputsThermal(solution);
     1186                //GetSolutionFromInputsThermal(solution);
     1187                GetSolutionFromInputsOneDof(solution, TemperatureEnum);
    11701188                break;
    11711189        case EnthalpyAnalysisEnum:
    1172                 GetSolutionFromInputsEnthalpy(solution);
     1190                //GetSolutionFromInputsEnthalpy(solution);
     1191                GetSolutionFromInputsOneDof(solution, EnthalpyEnum);
    11731192                break;
    11741193        #endif
     1194        #ifdef _HAVE_HYDROLOGY_
     1195        case HydrologyDCInefficientAnalysisEnum:
     1196                GetSolutionFromInputsOneDof(solution, SedimentHeadEnum);
     1197                break;
     1198        case HydrologyDCEfficientAnalysisEnum:
     1199                GetSolutionFromInputsOneDof(solution, EplHeadEnum);
     1200                break;
     1201  #endif
    11751202        default:
    11761203                _error_("analysis: " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    18801907                break;
    18811908        #endif
     1909        #ifdef _HAVE_HYDROLOGY_
     1910        case HydrologyDCInefficientAnalysisEnum:
     1911                InputUpdateFromSolutionHydrologyDCInefficient(solution);
     1912                break;
     1913        case HydrologyDCEfficientAnalysisEnum:
     1914                InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
     1915                break;
     1916        #endif
    18821917        default:
    18831918                _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    20852120                        return;
    20862121
    2087                 default:
     2122          case NodeSIdEnum:
     2123                        for(int i=0;i<NUMVERTICES;i++){
     2124                                values[i]=vector[nodes[i]->Sid()];
     2125                                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     2126                        }
     2127                        /*Add input to the element: */
     2128                        this->inputs->AddInput(new PentaP1Input(name,values));
     2129                       
     2130                        /*Free ressources:*/
     2131                        xDelete<int>(doflist);
     2132                        return;
     2133                       
     2134          default:
    20882135                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    20892136        }
    … …  
    21522199                                name==QmuMeltingEnum ||
    21532200                                name==GiaWEnum ||
    2154                                 name==GiadWdtEnum
    2155 
     2201                                name==GiadWdtEnum ||
     2202                                name==SedimentHeadEnum ||
     2203                                name==EplHeadEnum ||
     2204                                name==SedimentHeadOldEnum ||
     2205                                name==EplHeadOldEnum ||
     2206                                name==HydrologydcMaskEplactiveEnum ||
     2207                                name==WaterTransferEnum
     2208                               
    21562209                                ) {
    21572210                return true;
    … …  
    92679320#endif
    92689321#ifdef _HAVE_HYDROLOGY_
     9322
     9323/*FUNCTION Penta::CreateKMatrixHydrologyDCInefficient {{{*/
     9324ElementMatrix* Penta::CreateKMatrixHydrologyDCInefficient(void){
     9325
     9326        if (!IsOnBed()) return NULL;
     9327
     9328        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9329        ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCInefficient();
     9330        delete tria->material; delete tria;
     9331
     9332        /*clean up and return*/
     9333        return Ke;
     9334}
     9335/*}}}*/
     9336/*FUNCTION Penta::CreateKMatrixHydrologyDCEfficient {{{*/
     9337ElementMatrix* Penta::CreateKMatrixHydrologyDCEfficient(void){
     9338
     9339        if (!IsOnBed()) return NULL;
     9340
     9341        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9342        ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCEfficient();
     9343        delete tria->material; delete tria;
     9344
     9345        /*clean up and return*/
     9346        return Ke;
     9347}
     9348/*}}}*/
     9349/*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/
     9350ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){
     9351
     9352        if (!IsOnBed()) return NULL;
     9353
     9354        /*Call Tria function*/
     9355        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9356        ElementVector* pe=tria->CreatePVectorHydrologyDCInefficient();
     9357        delete tria->material; delete tria;
     9358
     9359        /*Clean up and return*/
     9360        return pe;
     9361}
     9362/*}}}*/
     9363/*FUNCTION Penta::CreatePVectorHydrologyDCEfficient {{{*/
     9364ElementVector* Penta::CreatePVectorHydrologyDCEfficient(void){
     9365
     9366        if (!IsOnBed()) return NULL;
     9367
     9368        /*Call Tria function*/
     9369        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9370        ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient();
     9371        delete tria->material; delete tria;
     9372
     9373        /*Clean up and return*/
     9374        return pe;
     9375}
     9376/*}}}*/
    92699377/*FUNCTION Penta::GetHydrologyDCInefficientHmax{{{*/
    92709378void  Penta::GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){
    9271         _error_("Hydrological stuff not suported in Penta");
     9379
     9380        if (!IsOnBed()) return;
     9381       
     9382        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9383        tria->GetHydrologyDCInefficientHmax(ph_max,innode);
     9384        delete tria->material; delete tria;
    92729385}
    92739386/*}}}*/
    92749387/*FUNCTION Penta::GetHydrologyTransfer{{{*/
    92759388void  Penta::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
    9276         _error_("Hydrological stuff not suported in Penta");
    9277 }
    9278 /*}}}*/
    9279 
     9389
     9390        if (!IsOnBed()) return;
     9391       
     9392        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9393        tria->GetHydrologyTransfer(transfer);
     9394        delete tria->material; delete tria;
     9395}
     9396/*}}}*/
     9397/*FUNCTION Penta::GetSolutionFromInputsOneDof {{{*/
     9398void Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){
     9399
     9400        const int    numdof=NDOF1*NUMVERTICES;
     9401
     9402        int          i;
     9403        int*         doflist=NULL;
     9404        IssmDouble   values[numdof];
     9405        IssmDouble   enum_value;
     9406        GaussPenta   *gauss=NULL;
     9407
     9408        /*Get dof list: */
     9409        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     9410        Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
     9411
     9412        gauss=new GaussPenta();
     9413        for(i=0;i<NUMVERTICES;i++){
     9414                /*Recover temperature*/
     9415                gauss->GaussVertex(i);
     9416                enum_input->GetInputValue(&enum_value,gauss);
     9417                values[i]=enum_value;
     9418        }
     9419
     9420        /*Add value to global vector*/
     9421        solution->SetValues(numdof,doflist,values,INS_VAL);
     9422
     9423        /*Free ressources:*/
     9424        delete gauss;
     9425        xDelete<int>(doflist);
     9426}
     9427/*{{{*/
    92809428/*FUNCTION Penta::HydrologyEPLGetActive {{{*/
    92819429void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
    9282         _error_("Hydrological stuff not suported in Penta");
     9430
     9431        if (!IsOnBed())return;
     9432       
     9433        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9434        tria->HydrologyEPLGetActive(active_vec);
     9435        delete tria->material; delete tria;
     9436
    92839437}
    92849438/*}}}*/
    92859439/*FUNCTION Penta::HydrologyEPLGetMask{{{*/
    92869440void  Penta::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){
    9287         _error_("Hydrological stuff not suported in Penta");
     9441
     9442        if (!IsOnBed())return;
     9443       
     9444        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     9445        tria->HydrologyEPLGetMask(vec_mask);
     9446        delete tria->material; delete tria;
     9447       
     9448}
     9449/*}}}*/
     9450/*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
     9451void  Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
     9452
     9453        const int   numdof   = NDOF1*NUMVERTICES;
     9454        const int   numdof2d = NDOF1*NUMVERTICES2D;
     9455        int*        doflist  = NULL;
     9456        bool        converged;
     9457        IssmDouble  values[numdof];
     9458        IssmDouble  residual[numdof];
     9459        IssmDouble  penalty_factor;
     9460        IssmDouble  kmax, kappa, h_max;
     9461        Penta      *penta    = NULL;
     9462
     9463        /*If not on bed, return*/
     9464        if (!IsOnBed()) return;
     9465
     9466        /*Get dof list: */
     9467        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     9468
     9469        /*Use the dof list to index into the solution vector and extrude it */
     9470        for(int i=0;i<numdof2d;i++){
     9471                values[i]         =solution[doflist[i]];
     9472                values[i+numdof2d]=values[i];
     9473                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     9474        }
     9475
     9476        /*If converged keep the residual in mind*/
     9477        this->inputs->GetInputValue(&converged,ConvergedEnum);
     9478
     9479        /*Get inputs*/
     9480        if(converged){
     9481                this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
     9482                this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
     9483               
     9484                kappa=kmax*pow(10.,penalty_factor);
     9485               
     9486                Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.       
     9487                for(int i=0;i<NUMVERTICES2D;i++){
     9488                        tria->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
     9489                        if(values[i]>h_max){
     9490                                residual[i]=kappa*(values[i]-h_max);
     9491                        }
     9492                        else{
     9493                                residual[i]=0.0;
     9494                        }
     9495                }
     9496                delete tria->material; delete tria;
     9497        }
     9498
     9499        /*Start looping over all elements above current element and update all inputs*/
     9500        penta=this;
     9501        for(;;){
     9502                /*Add input to the element: */
     9503                penta->inputs->AddInput(new PentaP1Input(SedimentHeadEnum,values));
     9504                penta->inputs->AddInput(new PentaP1Input(SedimentHeadResidualEnum,residual));
     9505
     9506                /*Stop if we have reached the surface*/
     9507                if (penta->IsOnSurface()) break;
     9508
     9509                /* get upper Penta*/
     9510                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     9511        }
     9512
     9513        /*Free ressources:*/
     9514        xDelete<int>(doflist);
    92889515}
    92899516/*}}}*/

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

    r15372 r15377  
    302302
    303303                #ifdef _HAVE_HYDROLOGY_
    304                 void    UpdateConstraints(void);
     304               
     305                ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
     306                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
     307                ElementVector* CreatePVectorHydrologyDCInefficient(void);
     308                ElementVector* CreatePVectorHydrologyDCEfficient(void);
     309               
    305310                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    306311                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
     312                void    GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type);
    307313                void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    308314                void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
     315                void    InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
     316                void    UpdateConstraints(void);
    309317                #endif
    310318                #ifdef _HAVE_THERMAL_

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

    r15372 r15377  
    64066406        int        *doflist  = NULL;
    64076407        bool        converged;
    6408         bool        isefficientlayer;
    64096408        IssmDouble  values[numdof];
    64106409        IssmDouble  residual[numdof];
    … …  
    64156414        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    64166415
    6417         /*Get the flag to know if the efficient layer is present*/
    6418         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    6419 
    64206416        /*Use the dof list to index into the solution vector: */
    64216417        for(int i=0;i<numdof;i++){
    … …  
    64296425        /*Get inputs*/
    64306426        if(converged){
    6431                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    64326427                this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
    64336428                this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
    … …  
    64616456                IssmDouble rho_ice,rho_water;
    64626457                IssmDouble thickness,bed;
    6463 
    64646458                /*Get the flag to the limitation method*/
    64656459                this->parameters->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);

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

    r15230 r15377  
    122122                                analysis_type==BedSlopeAnalysisEnum ||
    123123                                analysis_type==SurfaceSlopeAnalysisEnum ||
    124                                 analysis_type==BalancethicknessAnalysisEnum
     124                                analysis_type==BalancethicknessAnalysisEnum ||
     125                                analysis_type==HydrologyDCInefficientAnalysisEnum ||
     126                                analysis_type==HydrologyDCEfficientAnalysisEnum
    125127                                ){
    126128                if (dim==3){
Note: See TracChangeset for help on using the changeset viewer.