Changeset 23968

Timestamp:

05/31/19 10:49:15 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added moulins

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (1 diff)
• classes/Loads/Channel.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

@@ -42 +42 @@
                         loads->AddObject(new Channel(i+1,i,i,iomodel));
                 }
-
-                /*Free data: */
-        }
+        }
+
+        /*Create discrete loads for Moulins*/
+        CreateSingleNodeToElementConnectivity(iomodel);
+        for(int i=0;i<iomodel->numberofvertices;i++){
+                if (iomodel->domaintype!=Domain3DEnum){
+                        /*keep only this partition's nodes:*/
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(i+1,i,iomodel));
+                        }
+                }
+                else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(i+1,i,iomodel));
+                        }       
+                }
+        }
+        iomodel->DeleteData(1,"md.mesh.vertexonbase");
 }/*}}}*/
 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/

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

@@ -585 +585 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
-        IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
+        IssmDouble  A,B,n,phi,phi_0;
         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*/
@@ -610 +602 @@
         IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
         IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
+        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 
         Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
@@ -628 +621 @@
         GaussTria* gauss=new GaussTria();
         gauss->GaussEdgeCenter(index1,index2);
-
-        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*/
@@ -655 +642 @@
 
         /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
-        qc = - Ks * dphids;
+        IssmDouble qc = - Ks * dphids;
 
         /*d(phi - phi_m)/ds*/
-        dPw = dphids - dphimds;
+        IssmDouble dPw = dphids - dphimds;
 
         /*Compute f factor*/
-        fFactor = 0.;
+        IssmDouble 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));
-        }
-
-        _error_("STOP");
-
-        /*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*/
-                                                );
-                }
-        }
+        /*Compute total discharge*/
+        IssmDouble Q = -Kc*dphids;
+
+        /*Compute Pi and Xi*/
+        IssmDouble Pi = -C_W*c_t*rho_water*(Q+fFactor)*dPw;
+        IssmDouble Xi = fabs(Q*dphids) + fabs(lc * qc * dphids);
+
+        /*Compute closing rate*/
+        A=pow(B,-n);
+        IssmDouble vc = 2./pow(n,n)*A*this->S*pow(fabs(phi_0 - phi),n-1.)*(phi_0 - phi);
+
+        /*Compute new S based on Forward Euler (explicit)*/
+        this->S = this->S + dt*( (Xi - Pi)/(rho_ice*L) - vc);
+        if(this->S<0); this->S = 0.;
 
         /*Clean up and return*/