Channel Class Reference

#include <Channel.h>

Inheritance diagram for Channel:

Public Member Functions
	Channel ()
	Channel (int numericalflux_id, IssmDouble channelarea, int index, IoModel *iomodel)
	~Channel ()
Object *	copy ()
void	DeepEcho ()
void	Echo ()
int	Id ()
void	Marshall (char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
int	ObjectEnum ()
void	InputUpdateFromConstant (IssmDouble constant, int name)
void	InputUpdateFromConstant (int constant, int name)
void	InputUpdateFromConstant (bool constant, int name)
void	InputUpdateFromIoModel (int index, IoModel *iomodel)
void	InputUpdateFromMatrixDakota (IssmDouble *matrix, int nrows, int ncols, int name, int type)
void	InputUpdateFromVector (IssmDouble *vector, int name, int type)
void	InputUpdateFromVectorDakota (IssmDouble *vector, int name, int type)
void	Configure (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
void	CreateJacobianMatrix (Matrix< IssmDouble > *Jff)
void	CreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
void	CreatePVector (Vector< IssmDouble > *pf)
void	GetNodesLidList (int *lidlist)
void	GetNodesSidList (int *sidlist)
int	GetNumberOfNodes (void)
bool	IsPenalty (void)
void	PenaltyCreateJacobianMatrix (Matrix< IssmDouble > *Jff, IssmDouble kmax)
void	PenaltyCreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *kfs, IssmDouble kmax)
void	PenaltyCreatePVector (Vector< IssmDouble > *pf, IssmDouble kmax)
void	ResetHooks ()
void	SetCurrentConfiguration (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
void	SetwiseNodeConnectivity (int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)
void	SetChannelCrossSectionOld (void)
void	UpdateChannelCrossSection (void)
ElementVector *	CreatePVectorHydrologyGlaDS (void)
ElementMatrix *	CreateKMatrixHydrologyGlaDS (void)
void	WriteChannelCrossSection (IssmPDouble *values)
Public Member Functions inherited from Load
virtual	~Load ()
Public Member Functions inherited from Object
virtual	~Object ()

Data Fields
int	sid
int	id
Hook *	helement
Hook *	hnodes
Hook *	hvertices
Element *	element
Vertex **	vertices
Node **	nodes
Parameters *	parameters

Private Attributes
IssmDouble	S
IssmDouble	Sold
bool	boundary

Detailed Description

Definition at line 18 of file Channel.h.

Constructor & Destructor Documentation

Channel() [1/2]

Channel::Channel

(

)

Definition at line 30 of file Channel.cpp.

                {/*{{{*/
   this->id         = -1;
   this->sid        = -1;
   this->parameters = NULL;
   this->helement   = NULL;
   this->element    = NULL;
   this->hnodes     = NULL;
   this->hvertices  = NULL;
   this->nodes      = NULL;
}

Channel() [2/2]

Channel::Channel	(	int	numericalflux_id,
		IssmDouble	channelarea,
		int	index,
		IoModel *	iomodel
	)

Definition at line 41 of file Channel.cpp.

                                                                                {/*{{{*/
//Channel::Channel(int channel_id,int i,int index,IoModel* iomodel)
 
   this->id=channel_id;
   this->sid=channel_id-1;
   this->parameters = NULL;
   this->element    = NULL;
   this->nodes      = NULL;
 
   /*Set channel cross section to 0*/
   //this->S    = 0.;
   //this->Sold = 0.;
   this->S    = channelarea;
   this->Sold = channelarea;
 
   /*Get edge info*/
   int i1 = iomodel->faces[4*index+0];
   int i2 = iomodel->faces[4*index+1];
   int e1 = iomodel->faces[4*index+2];
   int e2 = iomodel->faces[4*index+3];
 
   if(e2==-1){
      this->boundary = true;
   }
   else{
      this->boundary = false;
   }
 
   /*Set Element hook (4th column may be -1 for boundary edges)*/
   this->helement  = new Hook(&e1,1);
 
   /*Set Vertices hooks (4th column may be -1 for boundary edges)*/
   int channel_vertex_ids[2];
   channel_vertex_ids[0]=i1;
   channel_vertex_ids[1]=i2;
   this->hvertices =new Hook(&channel_vertex_ids[0],2);
 
   /*Set Nodes hooks (!! Assumes P1 CG)*/
   int channel_node_ids[2];
   channel_node_ids[0]=i1;
   channel_node_ids[1]=i2;
   this->hnodes=new Hook(&channel_node_ids[0],2);
 
}/*}}}*/

~Channel()

Channel::~Channel

(

)

Definition at line 85 of file Channel.cpp.

                 {/*{{{*/
   this->parameters=NULL;
   delete helement;
   delete hnodes;
   delete hvertices;
}/*}}}*/

Member Function Documentation

copy()

Object * Channel::copy ( void  )

virtual

Implements Object.

Definition at line 93 of file Channel.cpp.

                      {/*{{{*/
 
   Channel* channel=NULL;
 
   channel=new Channel();
 
   /*copy fields: */
   channel->id=this->id;
   channel->S=this->S;
 
   /*point parameters: */
   channel->parameters=this->parameters;
 
   /*now deal with hooks and objects: */
   channel->hnodes    = (Hook*)this->hnodes->copy();
   channel->hvertices = (Hook*)this->hvertices->copy();
   channel->helement  = (Hook*)this->helement->copy();
 
   /*corresponding fields*/
   channel->nodes    = (Node**)channel->hnodes->deliverp();
   channel->vertices = (Vertex**)channel->hvertices->deliverp();
   channel->element  = (Element*)channel->helement->delivers();
 
   return channel;
}

DeepEcho()

void Channel::DeepEcho ( void  )

virtual

Implements Object.

Definition at line 119 of file Channel.cpp.

                             {/*{{{*/
 
   _printf_("Channel:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   S:  " << S << "\n");
   hnodes->DeepEcho();
   hvertices->DeepEcho();
   helement->DeepEcho();
   _printf_("   parameters\n");
   if(parameters)
    parameters->DeepEcho();
   else
    _printf_("      NULL\n");
}     

Echo()

void Channel::Echo ( void  )

virtual

Implements Object.

Definition at line 134 of file Channel.cpp.

                         {/*{{{*/
   _printf_("Channel:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   S:  " << S << "\n");
   hnodes->Echo();
   hvertices->Echo();
   helement->Echo();
   _printf_("   parameters: " << parameters << "\n");
}

Id()

int Channel::Id ( void  )

virtual

Implements Object.

Definition at line 144 of file Channel.cpp.

                       {/*{{{*/
   return id;
}

Marshall()

void Channel::Marshall ( char **  pmarshalled_data, int *  pmarshalled_data_size, int  marshall_direction  )

virtual

Implements Object.

Definition at line 148 of file Channel.cpp.

                                                                                                   { /*{{{*/
 
   _assert_(this);
 
   /*ok, marshall operations: */
   MARSHALLING_ENUM(ChannelEnum);
   MARSHALLING(id);
   MARSHALLING(S);
 
   if(marshall_direction==MARSHALLING_BACKWARD){
      this->hnodes      = new Hook();
      this->hvertices   = new Hook();
      this->helement    = new Hook();
   }
 
   this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
 
   /*corresponding fields*/
   nodes    =(Node**)this->hnodes->deliverp();
   vertices =(Vertex**)this->hvertices->deliverp();
   element  =(Element*)this->helement->delivers();
 
}

ObjectEnum()

int Channel::ObjectEnum ( void  )

virtual

Implements Object.

Definition at line 174 of file Channel.cpp.

                               {/*{{{*/
   return ChannelEnum;
}/*}}}*/

InputUpdateFromConstant() [1/3]

void Channel::InputUpdateFromConstant ( IssmDouble  constant, int  name  )

inline

Definition at line 55 of file Channel.h.

{/*Do nothing*/};

InputUpdateFromConstant() [2/3]

void Channel::InputUpdateFromConstant ( int  constant, int  name  )

inline

Definition at line 56 of file Channel.h.

{/*Do nothing*/};

InputUpdateFromConstant() [3/3]

void Channel::InputUpdateFromConstant ( bool  constant, int  name  )

inline

Definition at line 57 of file Channel.h.

{_error_("Not implemented yet!");}

InputUpdateFromIoModel()

void Channel::InputUpdateFromIoModel ( int  index, IoModel *  iomodel  )

inline

Definition at line 58 of file Channel.h.

{_error_("not implemented yet");};

InputUpdateFromMatrixDakota()

void Channel::InputUpdateFromMatrixDakota ( IssmDouble *  matrix, int  nrows, int  ncols, int  name, int  type  )

inline

Definition at line 59 of file Channel.h.

{/*Do nothing*/}

InputUpdateFromVector()

void Channel::InputUpdateFromVector ( IssmDouble *  vector, int  name, int  type  )

inline

Definition at line 60 of file Channel.h.

{/*Do nothing*/}

InputUpdateFromVectorDakota()

void Channel::InputUpdateFromVectorDakota ( IssmDouble *  vector, int  name, int  type  )

inline

Definition at line 61 of file Channel.h.

{/*Do nothing*/}

Configure()

void Channel::Configure ( Elements *  elements, Loads *  loads, Nodes *  nodes, Vertices *  vertices, Materials *  materials, Parameters *  parameters  )

virtual

Implements Load.

Definition at line 179 of file Channel.cpp.

                                                                                                                                               {/*{{{*/
 
   /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
    * datasets, using internal ids and offsets hidden in hooks: */
   hnodes->configure((DataSet*)nodesin);
   hvertices->configure((DataSet*)verticesin);
   helement->configure((DataSet*)elementsin);
 
   /*Initialize hooked fields*/
   this->nodes    = (Node**)hnodes->deliverp();
   this->vertices = (Vertex**)hvertices->deliverp();
   this->element  = (Element*)helement->delivers();
 
   /*point parameters to real dataset: */
   this->parameters=parametersin;
}

CreateJacobianMatrix()

void Channel::CreateJacobianMatrix ( Matrix< IssmDouble > *  Jff )

inlinevirtual

Implements Load.

Definition at line 65 of file Channel.h.

{_error_("Not implemented yet");};

CreateKMatrix()

void Channel::CreateKMatrix ( Matrix< IssmDouble > *  Kff, Matrix< IssmDouble > *  Kfs  )

virtual

Implements Load.

Definition at line 196 of file Channel.cpp.

                                                                            {/*{{{*/
 
   /*recover some parameters*/
   ElementMatrix* Ke=NULL;
   int analysis_type;
   this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
   switch(analysis_type){
      case HydrologyGlaDSAnalysisEnum:
         Ke = this->CreateKMatrixHydrologyGlaDS();
         break;
      default:
         _error_("Don't know why we should be here");
   }
 
   /*Add to global matrix*/
   if(Ke){
      Ke->AddToGlobal(Kff,Kfs);
      delete Ke;
   }
 
}

CreatePVector()

void Channel::CreatePVector ( Vector< IssmDouble > *  pf )

virtual

Implements Load.

Definition at line 219 of file Channel.cpp.

                                                  {/*{{{*/
 
   /*recover some parameters*/
   ElementVector* pe=NULL;
   int analysis_type;
   this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
   switch(analysis_type){
      case HydrologyGlaDSAnalysisEnum:
         pe = this->CreatePVectorHydrologyGlaDS();
         break;
      default:
         _error_("Don't know why we should be here");
   }
 
   /*Add to global matrix*/
   if(pe){
      pe->AddToGlobal(pf);
      delete pe;
   }
 
}

GetNodesLidList()

void Channel::GetNodesLidList ( int *  lidlist )

virtual

Implements Load.

Definition at line 242 of file Channel.cpp.

                                          {/*{{{*/
 
   _assert_(lidlist);
   _assert_(nodes);
 
   for(int i=0;i<NUMNODES;i++) lidlist[i]=nodes[i]->Lid();
}

GetNodesSidList()

void Channel::GetNodesSidList ( int *  sidlist )

virtual

Implements Load.

Definition at line 250 of file Channel.cpp.

                                          {/*{{{*/
 
   _assert_(sidlist);
   _assert_(nodes);
 
   for(int i=0;i<NUMNODES;i++) sidlist[i]=nodes[i]->Sid();
}

GetNumberOfNodes()

int Channel::GetNumberOfNodes ( void  )

virtual

Implements Load.

Definition at line 258 of file Channel.cpp.

                                   {/*{{{*/
   return NUMNODES;
}

IsPenalty()

bool Channel::IsPenalty ( void  )

virtual

Implements Load.

Definition at line 262 of file Channel.cpp.

                            {/*{{{*/
   return false;
}

PenaltyCreateJacobianMatrix()

void Channel::PenaltyCreateJacobianMatrix ( Matrix< IssmDouble > *  Jff, IssmDouble  kmax  )

inlinevirtual

Implements Load.

Definition at line 72 of file Channel.h.

{_error_("Not implemented yet");};

PenaltyCreateKMatrix()

void Channel::PenaltyCreateKMatrix ( Matrix< IssmDouble > *  Kff, Matrix< IssmDouble > *  kfs, IssmDouble  kmax  )

virtual

Implements Load.

Definition at line 266 of file Channel.cpp.

                                                                                                   {/*{{{*/
 
   /*No stiffness loads applied, do nothing: */
   return;
 
}

PenaltyCreatePVector()

void Channel::PenaltyCreatePVector ( Vector< IssmDouble > *  pf, IssmDouble  kmax  )

virtual

Implements Load.

Definition at line 273 of file Channel.cpp.

                                                                         {/*{{{*/
 
   /*No penalty loads applied, do nothing: */
   return;
 
}

ResetHooks()

void Channel::ResetHooks ( )

virtual

Implements Load.

Definition at line 280 of file Channel.cpp.

                         {/*{{{*/
 
   this->nodes=NULL;
   this->vertices=NULL;
   this->element=NULL;
   this->parameters=NULL;
 
   /*Get Element type*/
   this->hnodes->reset();
   this->hvertices->reset();
   this->helement->reset();
 
}

SetCurrentConfiguration()

void Channel::SetCurrentConfiguration ( Elements *  elements, Loads *  loads, Nodes *  nodes, Vertices *  vertices, Materials *  materials, Parameters *  parameters  )

virtual

Implements Load.

Definition at line 294 of file Channel.cpp.

                                                                                                                                                             {/*{{{*/
 
}

SetwiseNodeConnectivity()

void Channel::SetwiseNodeConnectivity ( int *  d_nz, int *  o_nz, Node *  node, bool *  flags, int *  flagsindices, int  set1_enum, int  set2_enum  )

virtual

Implements Load.

Definition at line 298 of file Channel.cpp.

                                                                                                                                {/*{{{*/
 
   /*Output */
   int d_nz = 0;
   int o_nz = 0;
 
   /*Loop over all nodes*/
   for(int i=0;i<this->GetNumberOfNodes();i++){
 
      if(!flags[this->nodes[i]->Lid()]){
 
         /*flag current node so that no other element processes it*/
         flags[this->nodes[i]->Lid()]=true;
 
         int counter=0;
         while(flagsindices[counter]>=0) counter++;
         flagsindices[counter]=this->nodes[i]->Lid();
 
         /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
         switch(set2_enum){
            case FsetEnum:
               if(nodes[i]->fsize){
                  if(this->nodes[i]->IsClone())
                   o_nz += 1;
                  else
                   d_nz += 1;
               }
               break;
            case GsetEnum:
               if(nodes[i]->gsize){
                  if(this->nodes[i]->IsClone())
                   o_nz += 1;
                  else
                   d_nz += 1;
               }
               break;
            case SsetEnum:
               if(nodes[i]->ssize){
                  if(this->nodes[i]->IsClone())
                   o_nz += 1;
                  else
                   d_nz += 1;
               }
               break;
            default: _error_("not supported");
         }
      }
   }
 
   /*Assign output pointers: */
   *pd_nz=d_nz;
   *po_nz=o_nz;
}

SetChannelCrossSectionOld()

void Channel::SetChannelCrossSectionOld

(

void

)

Definition at line 596 of file Channel.cpp.

                                                     {/*{{{*/
 
   this->Sold = this->S;
 
} /*}}}*/

UpdateChannelCrossSection()

void Channel::UpdateChannelCrossSection

(

void

)

Definition at line 601 of file Channel.cpp.

                                                     {/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement() || this->boundary){
      this->S = 0.;
      return;
   }
   _assert_(tria->FiniteElement()==P1Enum); 
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
 
   /*Intermediaries */
   IssmDouble  A,B,n,phi,phi_0,ks,Ngrad;
   IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  xyz_list_tria[3][3];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0]     ,this->vertices,NUMVERTICES);
   GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
 
   IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
   IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
   IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
   IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 
   Input2* h_input      = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
   Input2* H_input      = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input      = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* B_input      = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input      = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* ks_input     = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
   Input2* phi_input    = element->GetInput2(HydraulicPotentialEnum);         _assert_(phi_input);
 
   /*Get tangent vector*/
   IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
   IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
   IssmDouble Lt = sqrt(tx*tx+ty*ty);
   tx = tx/Lt;
   ty = ty/Lt;
 
   /*Evaluate fields on center of edge*/
   GaussTria* gauss=new GaussTria();
   gauss->GaussEdgeCenter(index1,index2);
 
   /*Get input values at gauss points*/
   phi_input->GetInputValue(&phi,gauss);
   phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
   h_input->GetInputValue(&h,gauss);
   ks_input->GetInputValue(&ks,gauss);
   B_input->GetInputValue(&B,gauss);
   n_input->GetInputValue(&n,gauss);
   b_input->GetInputValue(&b,gauss);
   b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
   H_input->GetInputValue(&H,gauss);
 
   /*Get values for a few potentials*/
   phi_0   = rho_water*g*b + rho_ice*g*H;
   dphids  = dphi[0]*tx + dphi[1]*ty;
   dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
   Ngrad   = fabs(dphids);
   if(Ngrad<AEPS) Ngrad = AEPS;
 
   /*d(phi - phi_m)/ds*/
   IssmDouble dPw = dphids - dphimds;
 
   /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
   IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
 
   /*Ice rate factor*/
   A=pow(B,-n);
 
   IssmDouble C = C_W*c_t*rho_water;
   IssmDouble Qprime = -kc * pow(Ngrad,BETA_C-2.)*dphids;
   IssmDouble N = phi_0 - phi;
 
   bool converged  = false;
   int  count      = 0;
 
   while(!converged){
 
      IssmDouble Snew = this->S;
 
      /*Compute f factor*/
      IssmDouble fFactor = 0.;
      if(this->S>0. || qc*dPw>0.){
         fFactor = lc * qc;
      }
 
      IssmDouble alpha = 1./(rho_ice*L)*(
               fabs(Qprime*pow(Snew,ALPHA_C-1.)*dphids)
               + C*Qprime*pow(Snew,ALPHA_C-1.)*dPw
               ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
 
      IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
 
      /*Solve ODE*/
      this->S = ODE1(alpha,beta,this->Sold,dt,2);
      if(this->S<0.) this->S = 0.;
      _assert_(!xIsNan<IssmDouble>(this->S)); 
 
      count++;
 
      if(fabs((this->S - Snew)/(Snew+AEPS))<1e-8  || count>=10) converged = true;
   }
 
   /*Clean up and return*/
   delete gauss;
}

CreatePVectorHydrologyGlaDS()

ElementVector * Channel::CreatePVectorHydrologyGlaDS

(

void

)

Definition at line 487 of file Channel.cpp.

                                                       {/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria* tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
   _assert_(tria->FiniteElement()==P1Enum); 
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
 
   /*Intermediaries */
   IssmDouble  Jdet,v2,Afactor,Bfactor,fFactor;
   IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds,dphi[2];
   IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks,Ngrad;
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  xyz_list_tria[3][3];
   const int   numnodes = NUMNODES;
 
   /*Initialize Element vector and other vectors*/
   ElementVector* pe = new ElementVector(this->nodes,NUMNODES,this->parameters);
   IssmDouble     basis[NUMNODES];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
   GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
 
   IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
   IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
 
   Input2* h_input      = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
   Input2* H_input      = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input      = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* B_input      = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input      = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* ks_input     = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
   Input2* phi_input    = element->GetInput2(HydraulicPotentialEnum);         _assert_(phi_input);
 
   /*Get tangent vector*/
   IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
   IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
   IssmDouble Lt = sqrt(tx*tx+ty*ty);
   tx = tx/Lt;
   ty = ty/Lt;
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 
      /*Get input values at gauss points*/
      b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
      phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
      h_input->GetInputValue(&h,gauss);
      ks_input->GetInputValue(&ks,gauss);
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      phi_input->GetInputValue(&phi,gauss);
      b_input->GetInputValue(&b,gauss);
      H_input->GetInputValue(&H,gauss);
 
      /*Get values for a few potentials*/
      phi_0   = rho_water*g*b + rho_ice*g*H;
      dphids  = dphi[0]*tx + dphi[1]*ty;
      dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
      Ngrad   = fabs(dphids);
      if(Ngrad<AEPS) Ngrad = AEPS;
 
      /*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
      IssmDouble Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
      /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
      qc = - Ks * dphids;
 
      /*d(phi - phi_m)/ds*/
      dPw = dphids - dphimds;
 
      /*Compute f factor*/
      fFactor = 0.;
      if(this->S>0. || qc*dPw>0.){
         fFactor = lc * qc;
      }
 
      /*Compute Afactor and Bfactor*/
      Afactor = C_W*c_t*rho_water;
      Bfactor = 1./L * (1./rho_ice - 1./rho_water);
 
      /*Compute closing rate*/
      /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
      A=pow(B,-n);
      v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
 
      for(int i=0;i<numnodes;i++){
         pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
         pe->values[i]+= + Jdet*gauss->weight*Afactor*Bfactor*fFactor*dphimds*basis[i];
      }
   }
 
   /*Clean up and return*/
   delete gauss;
   return pe;
}

CreateKMatrixHydrologyGlaDS()

ElementMatrix * Channel::CreateKMatrixHydrologyGlaDS

(

void

)

Definition at line 354 of file Channel.cpp.

                                                       {/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
   _assert_(tria->FiniteElement()==P1Enum); 
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
 
   /*Intermediaries */
   IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
   IssmDouble  A,B,n,phi_old,phi,phi_0,dPw,ks,Ngrad;
   IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  xyz_list_tria[3][3];
   const int   numnodes = NUMNODES;
 
   /*Initialize Element vector and other vectors*/
   ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
   IssmDouble     basis[NUMNODES];
   IssmDouble     dbasisdx[2*NUMNODES];
   IssmDouble     dbasisds[NUMNODES];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0]     ,this->vertices,NUMVERTICES);
   GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
 
   IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
   IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
   IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
 
   Input2* h_input      = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
   Input2* H_input      = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input      = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* B_input      = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input      = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* ks_input     = element->GetInput2(HydrologySheetConductivityEnum); _assert_(ks_input);
   Input2* phi_input    = element->GetInput2(HydraulicPotentialEnum);         _assert_(phi_input);
 
   /*Get tangent vector*/
   IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
   IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
   IssmDouble Lt = sqrt(tx*tx+ty*ty);
   tx = tx/Lt;
   ty = ty/Lt;
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
      tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
      dbasisds[0] = dbasisdx[0*2+0]*tx + dbasisdx[0*2+1]*ty;
      dbasisds[1] = dbasisdx[1*2+0]*tx + dbasisdx[1*2+1]*ty;
 
      /*Get input values at gauss points*/
      phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
      b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
      phi_input->GetInputValue(&phi,gauss);
      h_input->GetInputValue(&h,gauss);
      ks_input->GetInputValue(&ks,gauss);
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      b_input->GetInputValue(&b,gauss);
      H_input->GetInputValue(&H,gauss);
 
      /*Get values for a few potentials*/
      phi_0   = rho_water*g*b + rho_ice*g*H;
      dphids  = dphi[0]*tx + dphi[1]*ty;
      dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
      Ngrad   = fabs(dphids);
      if(Ngrad<AEPS) Ngrad = AEPS;
 
      /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
      IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
      IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
      /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
      qc = - Ks * dphids;
 
      /*d(phi - phi_m)/ds*/
      dPw = dphids - dphimds;
 
      /*Compute f factor*/
      fFactor = 0.;
      if(this->S>0. || qc*dPw>0.){
         fFactor = lc * qc;
      }
 
      /*Compute Afactor and Bfactor*/
      Afactor = C_W*c_t*rho_water;
      Bfactor = 1./L * (1./rho_ice - 1./rho_water);
      if(dphids>0){
         Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
      }
      else{
         Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
      }
 
      /*Diffusive term*/
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            /*GlaDSCoupledSolver.F90 line 1659*/
            Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
                     +Kc*dbasisds[i]*dbasisds[j]                               /*Diffusion term*/
                     - Afactor * Bfactor* Kc * dPw * basis[i] * dbasisds[j]    /*First part of Pi*/
                     + Afactor * fFactor * Bfactor * basis[i] * dbasisds[j]    /*Second part of Pi*/
                     + Xifactor* basis[i] * dbasisds[j]                        /*Xi term*/
                     );
         }
      }
 
      /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
      A=pow(B,-n);
      v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
         }
      }
   }
 
   /*Clean up and return*/
   delete gauss;
   return Ke;
}

WriteChannelCrossSection()

void Channel::WriteChannelCrossSection

(

IssmPDouble *

values

)

Definition at line 716 of file Channel.cpp.

                                                                   {/*{{{*/
 
   _assert_(values);
   values[this->sid] = reCast<IssmPDouble>(this->S);
}

Field Documentation

S

IssmDouble Channel::S

private

Definition at line 21 of file Channel.h.

Sold

IssmDouble Channel::Sold

private

Definition at line 22 of file Channel.h.

boundary

bool Channel::boundary

private

Definition at line 23 of file Channel.h.

sid

int Channel::sid

Definition at line 26 of file Channel.h.

id

int Channel::id

Definition at line 27 of file Channel.h.

helement

Hook* Channel::helement

Definition at line 30 of file Channel.h.

hnodes

Hook* Channel::hnodes

Definition at line 31 of file Channel.h.

hvertices

Hook* Channel::hvertices

Definition at line 32 of file Channel.h.

element

Element* Channel::element

Definition at line 35 of file Channel.h.

vertices

Vertex** Channel::vertices

Definition at line 36 of file Channel.h.

nodes

Node** Channel::nodes

Definition at line 37 of file Channel.h.

parameters

Parameters* Channel::parameters

Definition at line 38 of file Channel.h.

---

The documentation for this class was generated from the following files:

• src/c/classes/Loads/Channel.h
• src/c/classes/Loads/Channel.cpp