Changeset 14591


Ignore:
Timestamp:
04/15/13 16:54:51 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added sediment hydrology (temporary simplified version)

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

Legend:

Unmodified
Added
Removed

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

    r14589 r14591  
    34873487                        GetNodalFunctionsP1(&L[0], gauss);
    34883488                        D_scalar_trans=gauss->weight*Jdet;
    3489                         D_scalar_trans=D_scalar_trans;
    34903489
    34913490                        TripleMultiply(&L[0],numdof,1,0,
    … …  
    37173716                        GetNodalFunctionsP1(&L[0], gauss);
    37183717                        D_scalar_trans=gauss->weight*Jdet;
    3719                         D_scalar_trans=D_scalar_trans;
    37203718
    37213719                        TripleMultiply(&L[0],numdof,1,0,

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

    r14589 r14591  
    533533        IssmDouble  D,Jdet;
    534534        IssmDouble  xyz_list[NUMVERTICES][3];
    535         IssmDouble  L[1][3];
     535        IssmDouble  L[1][numdof];
    536536        GaussTria  *gauss = NULL;
    537537
    … …  
    56775677ElementMatrix* Tria::CreateKMatrixHydrology(void){
    56785678
     5679        int hydrology_model;
     5680        parameters->FindParam(&hydrology_model,HydrologyEnum);
     5681
     5682        switch(hydrology_model){
     5683                case HydrologyshreveEnum:
     5684                        return CreateKMatrixHydrologyShreve();
     5685                case HydrologydcEnum:
     5686                        return CreateKMatrixHydrologyDC();
     5687                default:
     5688                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
     5689        }
     5690
     5691}
     5692/*}}}*/
     5693/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/
     5694ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
     5695
    56795696        /*Constants*/
    56805697        const int  numdof=NDOF1*NUMVERTICES;
    … …  
    57785795}
    57795796/*}}}*/
    5780 /*FUNCTION Tria::CreatePVectorHydrology {{{*/
     5797/*FUNCTION Tria::CreateKMatrixHydrologyDC{{{*/
     5798ElementMatrix* Tria::CreateKMatrixHydrologyDC(void){
     5799
     5800        /*constants: */
     5801        const int    numdof=NDOF1*NUMVERTICES;
     5802
     5803        /* Intermediaries */
     5804        IssmDouble  D_scalar,Jdet;
     5805        IssmDouble  sediment_porosity,dt;
     5806        IssmDouble  xyz_list[NUMVERTICES][3];
     5807        IssmDouble  B[2][numdof];
     5808        IssmDouble L[NUMVERTICES];
     5809        IssmDouble  D[2][2];
     5810        GaussTria  *gauss = NULL;
     5811
     5812        /*Initialize Element matrix*/
     5813        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5814
     5815        /*Retrieve all inputs and parameters*/
     5816        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
     5817        sediment_porosity = matpar->GetSedimentPorosity();
     5818        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     5819
     5820        /* Start looping on the number of gaussian points: */
     5821        gauss=new GaussTria(2);
     5822        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5823
     5824                gauss->GaussPoint(ig);
     5825
     5826                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     5827
     5828                /*Diffusivity*/
     5829                D_scalar=sediment_porosity*gauss->weight*Jdet;
     5830                if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
     5831                D[0][0]=D_scalar; D[0][1]=0.;
     5832                D[1][0]=0.;       D[1][1]=D_scalar;
     5833                GetBHydro(&B[0][0],&xyz_list[0][0],gauss);
     5834                TripleMultiply(&B[0][0],2,numdof,1,
     5835                                        &D[0][0],2,2,0,
     5836                                        &B[0][0],2,numdof,0,
     5837                                        &Ke->values[0],1);
     5838
     5839                /*Transient*/
     5840                if(reCast<bool,IssmDouble>(dt)){
     5841                        GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     5842                        D_scalar=gauss->weight*Jdet;
     5843
     5844                        TripleMultiply(&L[0],numdof,1,0,
     5845                                                &D_scalar,1,1,0,
     5846                                                &L[0],1,numdof,0,
     5847                                                &Ke->values[0],1);
     5848                }
     5849        }
     5850
     5851        /*Clean up and return*/
     5852        delete gauss;
     5853        return Ke;
     5854}
     5855/*}}}*/
     5856/*FUNCTION Tria::CreatePVectorHydrology{{{*/
    57815857ElementVector* Tria::CreatePVectorHydrology(void){
     5858
     5859        int hydrology_model;
     5860        parameters->FindParam(&hydrology_model,HydrologyEnum);
     5861
     5862        switch(hydrology_model){
     5863                case HydrologyshreveEnum:
     5864                        return CreatePVectorHydrologyShreve();
     5865                case HydrologydcEnum:
     5866                        return CreatePVectorHydrologyDC();
     5867                default:
     5868                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
     5869        }
     5870
     5871}
     5872/*}}}*/
     5873/*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/
     5874ElementVector* Tria::CreatePVectorHydrologyShreve(void){
    57825875
    57835876        /*Constants*/
    … …  
    58205913                if(reCast<int,IssmDouble>(dt))for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
    58215914                else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
     5915        }
     5916
     5917        /*Clean up and return*/
     5918        delete gauss;
     5919        return pe;
     5920}
     5921/*}}}*/
     5922/*FUNCTION Tria::CreatePVectorHydrologyDC {{{*/
     5923ElementVector* Tria::CreatePVectorHydrologyDC(void){
     5924
     5925        /*Constants*/
     5926        const int    numdof=NDOF1*NUMVERTICES;
     5927
     5928        /*Intermediaries */
     5929        IssmDouble Jdet;
     5930        IssmDouble xyz_list[NUMVERTICES][3];
     5931        IssmDouble dt,gravity,scalar,water_head;
     5932        IssmDouble basis[numdof];
     5933        GaussTria* gauss=NULL;
     5934
     5935        /*Initialize Element vector*/
     5936        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     5937
     5938        /*Retrieve all inputs and parameters*/
     5939        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     5940        gravity = matpar->GetG();
     5941        Input* old_wh_input=NULL;
     5942        if(reCast<bool,IssmDouble>(dt)){
     5943                old_wh_input=inputs->GetInput(SedimentHeadEnum); _assert_(old_wh_input);
     5944        }
     5945
     5946        /* Start  looping on the number of gaussian points: */
     5947        gauss=new GaussTria(2);
     5948        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5949
     5950                gauss->GaussPoint(ig);
     5951
     5952                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     5953                GetNodalFunctions(basis, gauss);
     5954
     5955                /*Gravity term*/
     5956                scalar = Jdet*gauss->weight*gravity;
     5957                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
     5958                for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     5959
     5960                /*Transient term*/
     5961                if(reCast<bool,IssmDouble>(dt)){
     5962                        old_wh_input->GetInputValue(&water_head,gauss);
     5963                        scalar = Jdet*gauss->weight*water_head;
     5964                        for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     5965                }
    58225966        }
    58235967

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

    r14589 r14591  
    239239                #ifdef _HAVE_HYDROLOGY_
    240240                ElementMatrix* CreateKMatrixHydrology(void);
     241                ElementMatrix* CreateKMatrixHydrologyShreve(void);
     242                ElementMatrix* CreateKMatrixHydrologyDC(void);
    241243                ElementVector* CreatePVectorHydrology(void);
     244                ElementVector* CreatePVectorHydrologyShreve(void);
     245                ElementVector* CreatePVectorHydrologyDC(void);
    242246                void      CreateHydrologyWaterVelocityInput(void);
    243247                void      GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution);

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

    r14384 r14591  
    5656
    5757/*Reference Element numerics*/
     58/*FUNCTION TriaRef::GetBHydro {{{*/
     59void TriaRef::GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
     60        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     61         * For node i, Bi can be expressed in the actual coordinate system
     62         * by:
     63         *       Bi=[ dh/dx ]
     64         *          [ dh/dy ]
     65         * where h is the interpolation function for node i.
     66         *
     67         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
     68         */
     69
     70        int i;
     71        IssmDouble dbasis[NDOF2][NUMNODES];
     72
     73        /*Get dh1dh2dh3 in actual coordinate system: */
     74        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
     75
     76        /*Build B: */
     77        for (i=0;i<NUMNODES;i++){
     78                B[NDOF1*NUMNODES*0+NDOF1*i]=dbasis[0][i];
     79                B[NDOF1*NUMNODES*1+NDOF1*i]=dbasis[1][i];
     80        }
     81}
     82/*}}}*/
    5883/*FUNCTION TriaRef::GetBMacAyeal {{{*/
    5984void TriaRef::GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){

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

    r14384 r14591  
    2929                void GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss);
    3030                void GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss);
     31                void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
    3132                void GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof);
    3233                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss);

  • issm/trunk-jpl/src/c/solutions/hydrology_core.cpp

    r14577 r14591  
    6060
    6161                if (hydrology_model==HydrologyshreveEnum){
    62                         /*Compute hydrology solution: */
    6362                        if(VerboseSolution()) _pprintLine_("   computing water column");
    6463                        femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum);
    … …  
    8079                }
    8180                else if (hydrology_model==HydrologydcEnum){
    82                         //solver_linear(femmodel,modify_loads);
    83                         _error_("not implemented yet");
     81                        if(VerboseSolution()) _pprintLine_("   computing water head");
     82                        femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum);
     83                        solver_linear(femmodel);
    8484                        if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
    8585                                if(VerboseSolution()) _pprintLine_("   saving results ");
Note: See TracChangeset for help on using the changeset viewer.