Changeset 17114

Timestamp:

01/14/14 17:42:52 (11 years ago)

Author:

jbondzio

Message:

ADD: Building LSM shell pt2: Adding InputUpdateFromSolution, CreateKMatrix, Artificial Diffusivity etc. CHG: m/py code modified such that ice exists where LSF negative

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (9 diffs)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (1 diff)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/cores/transient_core.cpp (modified) (2 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• m/classes/mask.m (modified) (1 diff)
• m/classes/mask.py (modified) (1 diff)
• m/classes/transient.m (modified) (1 diff)
• m/classes/transient.py (modified) (1 diff)
• m/parameterization/setmask.m (modified) (1 diff)
• m/parameterization/setmask.py (modified) (1 diff)

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

@@ -59 +59 @@
 /*Finite element Analysis*/
 void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
-        _error_("not implemented yet");
+
+        /*parameters: */
+        bool save_results;
+
+        /*activate formulation: */
+        femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
+
+        /*recover parameters: */
+        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+        if(VerboseSolution()) _printf0_("call computational core:\n");
+        solutionsequence_linear(femmodel);
+
+        if(save_results){
+                if(VerboseSolution()) _printf0_("   saving results\n");
+                int outputs = MaskIceLevelsetEnum;
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+        }
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/
@@ -70 +87 @@
 }/*}}}*/
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
         /*Intermediaries */
-
-        int dim = 2; // solve for LSF in horizontal plane only
+        int i,j,dim = 2; // solve for LSF in horizontal plane only
         IssmDouble kappa;
         IssmDouble Jdet, dt, D_scalar;
-        IssmDouble u,v,um,vm,ub,vb,vx,vy;
+        IssmDouble h,hx,hy,hz;
+        IssmDouble vel,vx,vy,bx,by;
         IssmDouble* xyz_list = NULL;
 
@@ -86 +104 @@
         IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble     D[dim][dim];
+        IssmDouble     D[dim][dim], K[dim][dim];
 
         /*Retrieve all inputs and parameters*/
@@ -93 +111 @@
         Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
         Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
-        Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
-        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
         
         /* Start  looping on the number of gaussian points: */
@@ -117 +133 @@
                 GetB(B,element,xyz_list,gauss); 
                 GetBprime(Bprime,element,xyz_list,gauss); 
-                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
-                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
-                ub=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here
-                vb=0.;
-                D[0][0]=D_scalar*(vx+ub);
+                vx_input->GetInputValue(&vx,gauss); // in 3D case, add mesh velocity 
+                vy_input->GetInputValue(&vy,gauss); 
+                bx=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here
+                by=0.;
+                D[0][0]=D_scalar*(vx+bx);
                 D[0][1]=0.;
                 D[1][0]=0.;
-                D[1][1]=D_scalar*(vy+vb);
+                D[1][1]=D_scalar*(vy+by);
                 TripleMultiply(B,dim,numnodes,1,
                                         &D[0][0],dim,dim,0,
@@ -130 +146 @@
                                         &Ke->values[0],1);
 
-                /* Artificial Diffusion */
-                kappa=0.; //FIXME: insert suitable value for kappa
-                GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed?
-                D[0][0]=D_scalar*kappa;
-                D[1][1]=D_scalar*kappa;
-                TripleMultiply(Bprime,dim,numnodes,1,
-                                        &D[0][0],dim,dim,0,
-                                        Bprime,dim,numnodes,0,
-                                        &Ke->values[0],1);
+                /* Stabilization */
+                int stabilization=1;
+                switch(stabilization){
+                        case 0:
+                                // no stabilization, do nothing
+                                break;
+                        case 1:
+                                /* Artificial Diffusion */
+                                element->ElementSizes(&hx,&hy,&hz);
+                                vel=sqrt(vx*vx + vy*vy );
+                                h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); //FIXME: is this correct?
+
+                                kappa=h*vel/2.; //FIXME: insert suitable value for kappa
+                                //GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed?
+                                D[0][0]=D_scalar*kappa;
+                                D[0][1]=0.;
+                                D[1][0]=0.;
+                                D[1][1]=D_scalar*kappa;
+                                TripleMultiply(Bprime,dim,numnodes,1,
+                                                        &D[0][0],dim,dim,0,
+                                                        Bprime,dim,numnodes,0,
+                                                        &Ke->values[0],1);
+                                break;  
+                        case 2:
+                                /* Streamline Upwinding */
+                                element->ElementSizes(&hx,&hy,&hz);
+                                vel=sqrt(vx*vx + vy*vy )+1.e-14;
+                                h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) );
+                                K[0][0]=h/(2.*vel)*vx*vx;  K[0][1]=h/(2.*vel)*vx*vy;
+                                K[1][0]=h/(2.*vel)*vy*vx;  K[1][1]=h/(2.*vel)*vy*vy; 
+                                for(i=0;i<dim;i++) for(j=0;j<dim;j++) K[i][j] = D_scalar*K[i][j];
+
+                                //GetBprime(Bprime,element,xyz_list,gauss); 
+                                TripleMultiply(Bprime,dim,numnodes,1,
+                                                        &K[0][0],dim,dim,0,
+                                                        Bprime,dim,numnodes,0,
+                                                        &Ke->values[0],1);
+                                break;
+                        default:
+                                _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
+                }
         }
 
@@ -161 +209 @@
 
         /*Initialize Element vector*/
-        ElementVector* pe    = element->NewElementVector();
+        ElementVector* pe = element->NewElementVector();
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         
@@ -188 +236 @@
                 delete gauss;
         }
-        else{for(i=0;i<numnodes;i++) pe->values[i]=0.;} 
+        else
+                for(i=0;i<numnodes;i++) 
+                        pe->values[i]=0.; 
         return pe;
 }/*}}}*/
@@ -195 +245 @@
 }/*}}}*/
 void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        int meshtype;
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum);
+                        break;
+                case Mesh3DEnum:
+                        element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum);
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
 }/*}}}*/
 void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -146 +146 @@
                 virtual void   GetNodesLidList(int* sidlist)=0;
 
+                virtual int    Id()=0;
                 virtual int    Sid()=0;
                 virtual bool   IsNodeOnShelfFromFlags(IssmDouble* flags)=0; 

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

@@ -381 +381 @@
         /*clean-up*/
         delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::ElementSizes{{{*/
+void Tria::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){
+
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble xmin,ymin;
+        IssmDouble xmax,ymax;
+
+        /*Get xyz list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
+        ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
+
+        for(int i=1;i<NUMVERTICES;i++){
+                if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
+                if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
+                if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
+                if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
+        }
+
+        *hx=xmax-xmin;
+        *hy=ymax-ymin;
+        *hz=0.;
 }
 /*}}}*/

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

@@ -68 +68 @@
                 void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
                 void        Delta18oParameterization(void);
-                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
                 void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -50 +50 @@
         femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
         femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
+        femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
         femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
         if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
@@ -161 +162 @@
                 if(islevelset){
                         if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
-                        //analysis = new LevelsetAnalysis();
-                        //analysis->Core(femmodel);
-                        //delete analysis;
+                        analysis = new LevelsetAnalysis();
+                        analysis->Core(femmodel);
+                        delete analysis;
                 }
 

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -68 +68 @@
         parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
         parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
         parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
         parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));

issm/trunk-jpl/src/m/classes/mask.m

@@ -25 +25 @@
                         md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices 1]);
                         md = checkfield(md,'fieldname','mask.ice_levelset'        ,'size',[md.mesh.numberofvertices 1]);
-                        isice=(md.mask.ice_levelset>0);
-                        if any(sum(isice(md.mesh.elements),2)==0),
-                                warning('elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0');
-                        end
+                        %isice=(md.mask.ice_levelset>0);
+                        %if any(sum(isice(md.mesh.elements),2)==0),
+                %               warning('elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0');
+                %       end
                 end % }}}
                 function disp(obj) % {{{

issm/trunk-jpl/src/m/classes/mask.py

@@ -33 +33 @@
 
                 md = checkfield(md,'fieldname','mask.ice_levelset'        ,'size',[md.mesh.numberofvertices])
-                isice=numpy.array(md.mask.ice_levelset>0,int)
-                totallyicefree=(numpy.sum(isice[md.mesh.elements-1],axis=1)==0).astype(int)
-                if any(totallyicefree):
-                        raise TypeError("elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0")
+                #isice=numpy.array(md.mask.ice_levelset<0,int)
+                #totallyicefree=(numpy.sum(isice[md.mesh.elements-1],axis=1)==0).astype(int)
+                #if any(totallyicefree):
+                #       raise TypeError("elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0")
 
                 return md

issm/trunk-jpl/src/m/classes/transient.m

@@ -39 +39 @@
                 end % }}}
                 function list = defaultoutputs(self,md) % {{{
-
-                        list = {'SurfaceforcingsMassBalance'};
-
+                        if(self.ismasstransport)
+                                list = {'SurfaceforcingsMassBalance'};
+                        else
+                                list = {};
+                        end
                 end % }}}
                 function md = checkconsistency(obj,md,solution,analyses) % {{{

issm/trunk-jpl/src/m/classes/transient.py

@@ -38 +38 @@
         def defaultoutputs(self,md): # {{{
 
-                return ['SurfaceforcingsMassBalance']
+                if self.ismasstransport:
+                        return ['SurfaceforcingsMassBalance']
+                else:
+                        return []
 
         #}}}

issm/trunk-jpl/src/m/parameterization/setmask.m

@@ -42 +42 @@
 
 %level sets
-md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
+md.mask.ice_levelset=-1.*ones(md.mesh.numberofvertices,1);
 md.mask.groundedice_levelset=vertexongroundedice;
 md.mask.groundedice_levelset(find(vertexongroundedice==0.))=-1.;

issm/trunk-jpl/src/m/parameterization/setmask.py

@@ -42 +42 @@
 
         #level sets
-        md.mask.ice_levelset         = numpy.ones(md.mesh.numberofvertices,bool)
+        md.mask.ice_levelset         = -1.*numpy.ones((md.mesh.numberofvertices,1))
         md.mask.groundedice_levelset = -1.*numpy.ones((md.mesh.numberofvertices,1))
         md.mask.groundedice_levelset[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=1.