Changeset 23962


Ignore:
Timestamp:
05/30/19 13:56:00 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on channels

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

Legend:

Unmodified
Added
Removed

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

    r22983 r23962  
    328328        /*Clean up*/
    329329        xDelete<IssmDouble>(triabasis);
     330}
     331/*}}}*/
     332void TriaRef::GetSegmentNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list_tria,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
     333        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     334
     335        _assert_(index1>=0 && index1<3);
     336        _assert_(index2>=0 && index2<3);
     337
     338        /*Fetch number of nodes for this finite element*/
     339        int numnodes = this->NumberofNodes(finiteelement);
     340
     341        /*Get nodal functions*/
     342        IssmDouble* dtriabasis=xNew<IssmDouble>(2*numnodes);
     343        GetNodalFunctionsDerivatives(dtriabasis,xyz_list_tria,gauss,finiteelement);
     344
     345        switch(finiteelement){
     346                case P1Enum: case P1DGEnum:
     347                        dbasis[2*0+0] = dtriabasis[numnodes*0+index1];
     348                        dbasis[2*0+1] = dtriabasis[numnodes*1+index1];
     349                        dbasis[2*1+0] = dtriabasis[numnodes*0+index2];
     350                        dbasis[2*1+1] = dtriabasis[numnodes*1+index2];
     351                        break;
     352                default:
     353                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     354        }
     355
     356        /*Clean up*/
     357        xDelete<IssmDouble>(dtriabasis);
    330358}
    331359/*}}}*/

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

    r22983 r23962  
    2929                void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
    3030                void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement);
     31                void GetSegmentNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list_tria,Gauss* gauss, int index1,int index2,int finiteelement);
    3132                void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */};
    3233                void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);

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

    r23960 r23962  
    1414#include "../classes.h"
    1515/*}}}*/
    16 #define NUMNODES 2
     16#define NUMNODES    2
     17#define NUMVERTICES 2
    1718
    1819/*Channel constructors and destructor*/
    … …  
    3334        this->nodes      = NULL;
    3435
     36        /*Set channel cross section to 0*/
     37        this->S = 0.;
     38
    3539        /*Get edge info*/
    3640        int i1 = iomodel->faces[4*index+0];
    … …  
    7175        /*copy fields: */
    7276        channel->id=this->id;
     77        channel->S=this->S;
    7378
    7479        /*point parameters: */
    … …  
    9297        _printf_("Channel:\n");
    9398        _printf_("   id: " << id << "\n");
     99        _printf_("   S:  " << S << "\n");
    94100        hnodes->DeepEcho();
    95101        hvertices->DeepEcho();
    … …  
    105111        _printf_("Channel:\n");
    106112        _printf_("   id: " << id << "\n");
     113        _printf_("   S:  " << S << "\n");
    107114        hnodes->Echo();
    108115        hvertices->Echo();
    … …  
    122129        MARSHALLING_ENUM(ChannelEnum);
    123130        MARSHALLING(id);
     131        MARSHALLING(S);
    124132
    125133        if(marshall_direction==MARSHALLING_BACKWARD){
    … …  
    141149/*}}}*/
    142150int     Channel::ObjectEnum(void){/*{{{*/
    143 
    144151        return ChannelEnum;
    145 
    146 }
    147 /*}}}*/
     152}/*}}}*/
    148153
    149154/*Load virtual functions definitions:*/
    … …  
    323328
    324329/*Channel specific functions*/
    325 ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
    326 
    327         _error_("not implemented :( ");
     330ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
    328331
    329332        /*Initialize Element matrix and return if necessary*/
    330333        Tria*  tria=(Tria*)element;
    331334        if(!tria->IsIceInElement()) return NULL;
    332         ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters);
     335        _assert_(tria->FiniteElement()==P1Enum);
     336        int index1=tria->GetNodeIndex(nodes[0]);
     337        int index2=tria->GetNodeIndex(nodes[1]);
     338
     339        /*Intermediaries */
     340        IssmDouble  Jdet,v1;
     341        IssmDouble  A,B,n,phi_old,phi,phi_0;
     342        IssmDouble  H,h,b;
     343        IssmDouble  xyz_list[NUMVERTICES][3];
     344        IssmDouble  xyz_list_tria[3][3];
     345        const int   numnodes = NUMNODES;
     346
     347        /*Initialize Element vector and other vectors*/
     348        ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
     349        IssmDouble     basis[NUMNODES];
     350        IssmDouble     dbasisdx[2*NUMNODES];
     351        IssmDouble     dbasisds[NUMNODES];
     352
     353        /*Retrieve all inputs and parameters*/
     354        GetVerticesCoordinates(&xyz_list[0][0]     ,this->vertices,NUMVERTICES);
     355        GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
     356        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
     357        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
     358        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     359        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
     360        IssmDouble g         = element->FindParam(ConstantsGEnum);
     361        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
     362        Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
     363        Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
     364        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     365        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     366        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
     367        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
     368
     369        /*Get tangent vector*/
     370        IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
     371        IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
     372        tx = tx/sqrt(tx*tx+ty*ty);
     373        ty = ty/sqrt(tx*tx+ty*ty);
     374
     375        /* Start  looping on the number of gaussian points: */
     376        Gauss* gauss=new GaussTria(index1,index2,2);
     377        for(int ig=gauss->begin();ig<gauss->end();ig++){
     378                gauss->GaussPoint(ig);
     379
     380                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     381                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
     382                tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
     383
     384                /*Get input values at gauss points*/
     385                h_input->GetInputValue(&h,gauss);
     386                B_input->GetInputValue(&B,gauss);
     387                n_input->GetInputValue(&n,gauss);
     388                phi_input->GetInputValue(&phi,gauss);
     389                b_input->GetInputValue(&b,gauss);
     390                H_input->GetInputValue(&H,gauss);
     391
     392                /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
     393                phi_0 = rho_water*g*b + rho_ice*g*H;
     394                A=pow(B,-n);
     395                v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));
     396                for(int i=0;i<numnodes;i++){
     397                        for(int j=0;j<numnodes;j++){
     398                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
     399                        }
     400                }
     401
     402                /*Transient term if dt>0*/
     403                if(dt>0.){
     404                }
     405        }
    333406
    334407        /*Clean up and return*/
     408        delete gauss;
     409        return Ke;
     410}
     411/*}}}*/
     412ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/
     413
     414        /*Initialize Element matrix and return if necessary*/
     415        Tria* tria=(Tria*)element;
     416        if(!tria->IsIceInElement()) return NULL;
     417        _assert_(tria->FiniteElement()==P1Enum);
     418        int index1=tria->GetNodeIndex(nodes[0]);
     419        int index2=tria->GetNodeIndex(nodes[1]);
     420
     421        /*Intermediaries */
     422        IssmDouble  Jdet,v2;
     423        IssmDouble  A,B,n,phi_old,phi,phi_0;
     424        IssmDouble  H,h,b;
     425        IssmDouble  xyz_list[NUMVERTICES][3];
     426        const int   numnodes = NUMNODES;
     427
     428        /*Initialize Element vector and other vectors*/
     429        ElementVector* pe = new ElementVector(this->nodes,NUMNODES,this->parameters);
     430        IssmDouble     basis[NUMNODES];
     431
     432        /*Retrieve all inputs and parameters*/
     433        GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     434        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
     435        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
     436        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     437        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
     438        IssmDouble g         = element->FindParam(ConstantsGEnum);
     439        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
     440        Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
     441        Input* b_input      = element->GetInput(BedEnum); _assert_(b_input);
     442        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
     443        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     444        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
     445        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
     446
     447        /* Start  looping on the number of gaussian points: */
     448        Gauss* gauss=new GaussTria(index1,index2,2);
     449        for(int ig=gauss->begin();ig<gauss->end();ig++){
     450                gauss->GaussPoint(ig);
     451
     452                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     453                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
     454
     455                /*Get input values at gauss points*/
     456                h_input->GetInputValue(&h,gauss);
     457                B_input->GetInputValue(&B,gauss);
     458                n_input->GetInputValue(&n,gauss);
     459                phi_input->GetInputValue(&phi,gauss);
     460                b_input->GetInputValue(&b,gauss);
     461                H_input->GetInputValue(&H,gauss);
     462
     463                /*Compute closing rate*/
     464                /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
     465                phi_0 = rho_water*g*b + rho_ice*g*H;
     466                A=pow(B,-n);
     467                v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
     468
     469                for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
     470
     471                /*Transient term if dt>0*/
     472                if(dt>0.){
     473                        //phiold_input->GetInputValue(&phi_old,gauss);
     474                        //for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
     475                }
     476        }
     477
     478        /*Clean up and return*/
     479        delete gauss;
    335480        return pe;
    336481}
    337482/*}}}*/
    338 ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
    339 
    340         _error_("not implemented :( ");
    341 
    342         /*Initialize Element matrix and return if necessary*/
    343         Tria*  tria=(Tria*)element;
    344         if(!tria->IsIceInElement()) return NULL;
    345         ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters);
    346 
    347         /*Clean up and return*/
    348         return Ke;
    349 }
    350 /*}}}*/

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

    r23960 r23962  
    1717
    1818class Channel: public Load {
     19
     20        private:
     21                IssmDouble S;
    1922
    2023        public:
Note: See TracChangeset for help on using the changeset viewer.