Changeset 27035 for issm/trunk/src/c/analyses/FreeSurfaceTopAnalysis.cpp

Timestamp:

06/01/22 05:01:48 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 27033

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c/analyses/FreeSurfaceTopAnalysis.cpp (modified) (12 diffs)

issm/trunk/src/c/analyses/FreeSurfaceTopAnalysis.cpp

@@ -80 +80 @@
         iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
         iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
+        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
         if(iomodel->domaintype!=Domain2DhorizontalEnum){
                 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
@@ -111 +112 @@
 }/*}}}*/
 ElementMatrix* FreeSurfaceTopAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
-_error_("Not implemented");
+        _error_("Not implemented");
 }/*}}}*/
 ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
@@ -120 +121 @@
         IssmDouble *xyz_list  = NULL;
         IssmDouble  Jdet,D_scalar,dt,h;
-        IssmDouble  vel,vx,vy;
+        IssmDouble  vel,vx,vy,tau;
 
         /*Get top element*/
@@ -195 +196 @@
                 }
 
-                if(stabilization==2){
+                if(stabilization==1){
+                        /*artifical diffusion*/
+                        if(dim==1){
+                                vx_input->GetInputAverage(&vx);
+                                D[0]=h/2.*fabs(vx);
+                        }
+                        else{
+                                vx_input->GetInputAverage(&vx);
+                                vy_input->GetInputAverage(&vy);
+
+                                D[0*dim+0]=h/2.0*fabs(vx);
+                                D[1*dim+1]=h/2.0*fabs(vy);
+                        }
+                }
+                else if(stabilization==2){
                         /*Streamline upwinding*/
                         if(dim==1){
@@ -209 +224 @@
                         }
                 }
-                else if(stabilization==1){
-                        /*SSA*/
+                else if(stabilization==5){
+                        /*SUPG*/
                         if(dim==1){
                                 vx_input->GetInputAverage(&vx);
-                                D[0]=h/2.*fabs(vx);
+                                tau=h/(2.*fabs(vx)+1e-10);
                         }
                         else{
                                 vx_input->GetInputAverage(&vx);
                                 vy_input->GetInputAverage(&vy);
-
-                                D[0*dim+0]=h/2.0*fabs(vx);
-                                D[1*dim+1]=h/2.0*fabs(vy);
+                                tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
                         }
                 }
@@ -229 +242 @@
                                         for(int j=0;j<numnodes;j++){
                                                 Ke->values[i*numnodes+j] += (
-                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
-                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
-                                                                        );
+                                                                dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
+                                                                dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
+                                                                );
                                         }
                                 }
@@ -237 +250 @@
                         else{
                                 for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
+                        }
+                }
+                else if(stabilization==5){
+                        D_scalar=gauss->weight*Jdet*dt;
+                        if(dim==2){
+                                for(int i=0;i<numnodes;i++){
+                                        for(int j=0;j<numnodes;j++){
+                                                Ke->values[i*numnodes+j]+=tau*D_scalar*
+                                                        (vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
+                                                        (vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
+                                        }
+                                }
+                        }
+                        else{
+                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
                         }
                 }
@@ -251 +279 @@
 }/*}}}*/
 ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
+
         /*Intermediaries*/
-        int         domaintype,dim;
+        int         domaintype,dim,stabilization;
         IssmDouble  Jdet,dt;
-        IssmDouble  ms,surface,vz;
+        IssmDouble  ms,surface,vx,vy,vz,tau;
         Element*    topelement = NULL;
         IssmDouble *xyz_list  = NULL;
@@ -284 +313 @@
         ElementVector* pe    = topelement->NewElementVector();
         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
         /*Retrieve all inputs and parameters*/
         topelement->GetVerticesCoordinates(&xyz_list);
         topelement->FindParam(&dt,TimesteppingTimeStepEnum);
-        Input* ms_input      = topelement->GetInput(SmbMassBalanceEnum);  _assert_(ms_input);
-        Input* surface_input = topelement->GetInput(SurfaceEnum);                     _assert_(surface_input);
-        Input* vz_input      = NULL;
+        Input *ms_input      = topelement->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
+        Input *surface_input = topelement->GetInput(SurfaceEnum);        _assert_(surface_input);
+        Input *vz_input      = NULL;
+        Input *vx_input      = NULL;
+        Input *vy_input      = NULL;
         switch(dim){
-                case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break;
-                case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break;
+                case 1: 
+                        vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
+                        vz_input = topelement->GetInput(VyEnum) ; _assert_(vz_input); 
+                        break;
+                case 2: 
+                        vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
+                        vy_input = topelement->GetInput(VyEnum); _assert_(vy_input);
+                        vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); 
+                        break;
                 default: _error_("not implemented");
         }
+        IssmDouble h = topelement->CharacteristicLength();
 
         /*Initialize mb_correction to 0, do not forget!:*/
@@ -312 +352 @@
         }
 
+        if(stabilization==5){
+                /*SUPG*/
+                topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                if(dim==1){
+                        vx_input->GetInputAverage(&vx);
+                        tau=h/(2.*fabs(vx)+1e-10);
+                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
+                }
+                else{
+                        vx_input->GetInputAverage(&vx);
+                        vy_input->GetInputAverage(&vy);
+                        tau=h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
+                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+                }
+        }
+
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
@@ -321 +378 @@
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-           _error_("not implemented yet");
+        _error_("not implemented yet");
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
@@ -329 +386 @@
 
         element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
+
+        /*Now, we need to do some "processing"*/
+        int numvertices = element->GetNumberOfVertices();
+        int        migration_style;
+
+        IssmDouble* surface = xNew<IssmDouble>(numvertices);
+        IssmDouble* newsurface = xNew<IssmDouble>(numvertices);
+        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+        IssmDouble* base = xNew<IssmDouble>(numvertices);
+        IssmDouble* bed = xNew<IssmDouble>(numvertices);
+        IssmDouble* phi = xNew<IssmDouble>(numvertices);
+        IssmDouble* sealevel = xNew<IssmDouble>(numvertices);
+
+        IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
+        IssmDouble rho_ice      = element->FindParam(MaterialsRhoIceEnum);
+        IssmDouble rho_water    = element->FindParam(MaterialsRhoSeawaterEnum);
+
+        bool isgroundingline;
+        element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
+        element->FindParam(&migration_style,GroundinglineMigrationEnum);
+        if(isgroundingline) element->GetInputListOnVertices(&bed[0],BedEnum);
+
+        element->GetInputListOnVertices(&base[0],BaseEnum);
+        element->GetInputListOnVertices(&surface[0],SurfaceEnum);
+        element->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
+        element->GetInputListOnVertices(&sealevel[0],SealevelEnum);
+
+        for(int i=0;i<numvertices;i++){
+                newsurface[i]=surface[i];
+                thickness[i]=surface[i]-base[i];
+                /*Check solution*/
+                if(xIsNan<IssmDouble>(thickness[i])) _error_("NaN found in solution vector");
+                if(xIsInf<IssmDouble>(thickness[i])) _error_("Inf found in solution vector");
+
+                /* check for thickness<minthickness */
+                if(thickness[i]<minthickness){ 
+                        thickness[i]=minthickness;
+                        if(phi[i]>0.){
+                                if(base[i]<=bed[i]) base[i] = bed[i];
+                                newsurface[i] = base[i]+minthickness;
+                        }else{
+                                // assume floatation condition
+                                newsurface[i] = (1.-rho_ice/rho_water)*minthickness;
+                                base[i] = -rho_ice/rho_water*minthickness;
+                        }
+                }
+
+                /* update thickness */
+                thickness[i]=newsurface[i]-base[i];
+                /* some checks */
+                if(thickness[i]<0.) _error_("thickness<0");
+                if(newsurface[i]<base[i]) _error_("surface<base");
+        }
+
+        /* update inputs */
+        element->AddInput(BaseEnum,base,element->GetElementType());
+        element->AddInput(SurfaceEnum,newsurface,element->GetElementType());
+        element->AddInput(ThicknessEnum,thickness,element->GetElementType());
+
+        /* Free resources */
+        xDelete<IssmDouble>(newsurface);
+        xDelete<IssmDouble>(surface);
+        xDelete<IssmDouble>(thickness);
+        xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(bed);
+        xDelete<IssmDouble>(phi);
+        xDelete<IssmDouble>(sealevel);
 }/*}}}*/
 void           FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/