Changeset 20117

Timestamp:

02/11/16 10:59:52 (9 years ago)

Author:

jbondzio

Message:

CHG: Enabling extrapolation for level-set-method in 3d; Cleanup of transient_core

Location:

issm/trunk-jpl

Files:

• src/c/Makefile.am (modified) (1 diff)
• src/c/analyses/ExtrapolationAnalysis.cpp (modified) (12 diffs)
• src/c/analyses/ExtrapolationAnalysis.h (modified) (1 diff)
• src/c/classes/Node.cpp (modified) (1 diff)
• src/c/cores/cores.h (modified) (1 diff)
• src/c/cores/levelset_core.cpp (added)
• src/c/cores/transient_core.cpp (modified) (2 diffs)
• test/Archives/Archive802.nc (modified) ( previous)
• test/Archives/Archive803.nc (modified) ( previous)
• test/Archives/Archive805.nc (modified) ( previous)

issm/trunk-jpl/src/c/Makefile.am

@@ -257 +257 @@
                                         ./cores/damage_core.cpp\
                                         ./cores/levelsetfunctionslope_core.cpp\
+                                        ./cores/levelset_core.cpp\
                                         ./modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp\
                                         ./classes/Loads/Riftfront.cpp\

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

@@ -59 +59 @@
         bool save_results;
         int extvar_enum; 
-   femmodel->parameters->FindParam(&extvar_enum, ExtrapolationVariableEnum);
+        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+        femmodel->parameters->FindParam(&extvar_enum, ExtrapolationVariableEnum);
 
         /*activate formulation: */
@@ -67 +68 @@
         solutionsequence_linear(femmodel);
 
-        save_results=true;
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                int outputs[2] = {VxEnum,VyEnum};
-                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+                femmodel->RequestedOutputsx(&femmodel->results,&extvar_enum,1);
         }
 }/*}}}*/
@@ -84 +83 @@
 ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
-        if(!element->IsOnBase()) return NULL;
-        Element* basalelement = element->SpawnBasalElement();
-
         /*Intermediaries */
-        int                domaintype,dim;
-        int        i,row,col,stabilization;
-        bool       extrapolatebydiffusion = true;
+        int     dim, domaintype, extrapolationcase;
+        int     i,row,col,stabilization;
+        bool    extrapolatebydiffusion;
         IssmDouble Jdet,D_scalar,h;
         IssmDouble norm_dlsf;
         IssmDouble hx,hy,hz,kappa;
-        IssmDouble* xyz_list = NULL;
-
-        /*Get problem dimension*/
-        basalelement->FindParam(&domaintype,DomainTypeEnum);
+        IssmDouble*     xyz_list = NULL;
+        Element*                workelement=NULL;
+
+        /*Get problem case*/
+        extrapolationcase=GetExtrapolationCase(element);
+        switch(extrapolationcase){
+                case 0:
+                        if(!element->IsOnBase()) return NULL; 
+                        workelement = element->SpawnBasalElement(); 
+                        break;
+                case 1: case 2: case 3: workelement=element; break;
+        }
+        /* get extrapolation dimension */
+        workelement->FindParam(&domaintype,DomainTypeEnum);
         switch(domaintype){
-                case Domain2DverticalEnum:   dim = 1; break;
-                case Domain2DhorizontalEnum: dim = 2; break;
-                case Domain3DEnum:           dim = 2; break;
-                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+                case Domain2DverticalEnum: dim=1; break;
+                case Domain2DhorizontalEnum: dim=2; break;
+                case Domain3DEnum: dim=3; break;
         }
 
         /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = basalelement->GetNumberOfNodes();
+        int numnodes = workelement->GetNumberOfNodes();
 
         /*Initialize Element vector and other vectors*/
-        ElementMatrix* Ke     = basalelement->NewElementMatrix();
+        ElementMatrix* Ke     = workelement->NewElementMatrix();
         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
@@ -117 +122 @@
 
         /*Retrieve all inputs and parameters*/
-        Input* lsf_slopex_input=basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-        Input* lsf_slopey_input=basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-        basalelement->GetVerticesCoordinates(&xyz_list);
-        h = basalelement->CharacteristicLength();
+        Input* lsf_slopex_input=workelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+        Input* lsf_slopey_input=workelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
+        workelement->GetVerticesCoordinates(&xyz_list);
 
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=basalelement->NewGauss(2);
+        Gauss* gauss=workelement->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
                 gauss->GaussPoint(ig);
 
-                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,basalelement,xyz_list,gauss);
-                GetBprime(Bprime,basalelement,xyz_list,gauss);
+                workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,workelement,xyz_list,gauss,dim);
+                GetBprime(Bprime,workelement,xyz_list,gauss,dim);
                 
                 D_scalar=gauss->weight*Jdet;
 
+                extrapolatebydiffusion=true;
                 if(extrapolatebydiffusion){
                         /* diffuse values outward */
                         for(row=0;row<dim;row++)
                                 for(col=0;col<dim;col++)
-                                        if(row==col)
+                                        if(row==col && row<2) //extrapolate only in xy-plane
                                                 D[row*dim+col] = D_scalar;
                                         else
@@ -151 +156 @@
                         /* Get normal on ice boundary */
                         lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-                        lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+                        if(dim>1)
+                                lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+                        if(dim>2)
+                                dlsf[2]=0.;
                         norm_dlsf=0.;
                         for(i=0;i<dim;i++)      norm_dlsf+=dlsf[i]*dlsf[i]; 
@@ -172 +180 @@
                                                 &Ke->values[0],1);
 
-                        /* Stabilization *//*{{{*/
+                        /* stabilization *//*{{{*/
+                        /* do not use streamline upwinding for extrapolation: it yields oscillating results due to diffusion along normal direction, but none across */
                         stabilization=1;
                         if (stabilization==0){/* no stabilization, do nothing*/}
                         else if(stabilization==1){
                                 /* Artificial Diffusion */
-                                basalelement->ElementSizes(&hx,&hy,&hz);
-                                h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+                                workelement->ElementSizes(&hx,&hy,&hz);
+                                h=sqrt(pow(hx*normal[0],2) + pow(hy*normal[1],2));
                                 kappa=h/2.+1.e-14; 
                                 for(row=0;row<dim;row++)
@@ -191 +200 @@
                                                         &Ke->values[0],1);
                         }
-                        else if(stabilization==2){
-                                /*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
-                                for(row=0;row<dim;row++)
-                                        for(col=0;col<dim;col++)
-                                                D[row*dim+col]=h/(2.*1.)*normal[row]*normal[col];
-
-                                TripleMultiply(Bprime,dim,numnodes,1,
-                                                        D,dim,dim,0,
-                                                        Bprime,dim,numnodes,0,
-                                                        &Ke->values[0],1);
-                        }/*}}}*/
+                        /*}}}*/
                 }
         }/*}}}*/
@@ -213 +212 @@
         xDelete<IssmDouble>(normal);
         delete gauss;
-        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};
         return Ke;
 
@@ -219 +218 @@
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
 
-        if(!element->IsOnBase()) return NULL;
-        Element* basalelement = element->SpawnBasalElement();
-
         /*Intermediaries */
-        int i, domaintype;
-        
+        Element* workelement=NULL;
+
+        /*Get problem dimension*/
+        int extrapolationcase=GetExtrapolationCase(element);
+        switch(extrapolationcase){
+                case 0: 
+                        if(!element->IsOnBase()) return NULL; 
+                        workelement = element->SpawnBasalElement(); 
+                        break;
+                case 1: case 2: case 3: workelement=element; break;
+        }
+
         /*Fetch number of nodes */
-        int numnodes = basalelement->GetNumberOfNodes();
+        int numnodes = workelement->GetNumberOfNodes();
 
         /*Initialize Element vector*/
-        ElementVector* pe = basalelement->NewElementVector();
-        for(i=0;i<numnodes;i++) 
+        ElementVector* pe = workelement->NewElementVector();
+        for(int i=0;i<numnodes;i++) 
                 pe->values[i]=0.; 
 
-        basalelement->FindParam(&domaintype,DomainTypeEnum);
-        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};
         return pe;
 }/*}}}*/
-void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
          * For node i, Bi can be expressed in the actual coordinate system
@@ -256 +261 @@
 
         /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = basis[i];
-                B[numnodes*1+i] = basis[i];
-        }
+        for(int i=0;i<numnodes;i++)
+                for(int d=0;d<dim;d++)
+                        B[numnodes*d+i] = basis[i];
 
         /*Clean-up*/
         xDelete<IssmDouble>(basis);
 }/*}}}*/
-void           ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
          * For node i, Bi' can be expressed in the actual coordinate system
@@ -279 +283 @@
 
         /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B': */
-        for(int i=0;i<numnodes;i++){
-                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
+        for(int i=0;i<numnodes;i++)
+                for(int d=0;d<dim;d++)
+                        Bprime[numnodes*d+i] = dbasis[d*numnodes+i];
 
         /*Clean-up*/
@@ -300 +303 @@
 void           ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
+        int extrapolationvariable, extrapolationcase;
+        extrapolationcase=GetExtrapolationCase(element);
+        element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
+        switch(extrapolationcase){
+                case 0:
+                        element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
+                        break;
+                case 1:
+                        element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
+                        break;
+                case 2:
+                        element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable);
+                        break;
+                case 3:
+                        element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
+                        break;
+        }
+}/*}}}*/
+int                             ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/
+
+        /* Get case of extrapolation, depending on domain quality, and extrapolation variable */
         int domaintype, extrapolationvariable;
+        int extrapolationcase;
         element->FindParam(&domaintype,DomainTypeEnum);
-        element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
         switch(domaintype){
-                case Domain2DhorizontalEnum:
-                        element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
-                        break;
-                case Domain3DEnum:
-                        element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable);
-                        break;
-                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-        }
+                case Domain2DverticalEnum: extrapolationcase=0; break;
+                case Domain2DhorizontalEnum: extrapolationcase=1;break;
+                case Domain3DEnum:  
+                        element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
+                        if(extrapolationvariable==ThicknessEnum) extrapolationcase=2; // scalar fields that are constant along z-axis
+                        else extrapolationcase=3; // scalar fields that vary along z-axis
+                        break;
+        }
+        return extrapolationcase;
 }/*}}}*/
 void           ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/

issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h

@@ -26 +26 @@
         ElementMatrix* CreateKMatrix(Element* element);
         ElementVector* CreatePVector(Element* element);
-        void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-        void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+        void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
+        void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
         void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
         void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
         void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+        int                             GetExtrapolationCase(Element* element);
         void           SetConstraintsOnIce(Element* element);
         void           UpdateConstraints(FemModel* femmodel);

issm/trunk-jpl/src/c/classes/Node.cpp

@@ -95 +95 @@
                                 analysis_enum==HydrologyDCInefficientAnalysisEnum ||
                                 analysis_enum==HydrologyDCEfficientAnalysisEnum ||
-                                analysis_enum==LevelsetAnalysisEnum ||
-                                analysis_enum==ExtrapolationAnalysisEnum
+                                analysis_enum==LevelsetAnalysisEnum
                                 ){
                 if(iomodel->domaintype!=Domain2DhorizontalEnum){

issm/trunk-jpl/src/c/cores/cores.h

@@ -24 +24 @@
 void surfaceslope_core(FemModel* femmodel);
 void levelsetfunctionslope_core(FemModel* femmodel);
+void levelset_core(FemModel* femmodel);
 void bedslope_core(FemModel* femmodel);
 void meshdeformation_core(FemModel* femmodel);

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

@@ -102 +102 @@
                 if(isdamageevolution) damage_core(femmodel);
 
-                if(islevelset){
-                        if(iscalving) Calvingx(femmodel);
-                        if(VerboseSolution()) _printf0_("   computing levelset transport\n");
-                        /* smoothen slope of lsf for computation of normal on ice domain*/
-                        levelsetfunctionslope_core(femmodel);
+                if(islevelset)  levelset_core(femmodel);
 
-                        /* extrapolate velocities onto domain with no ice */
-                        Analysis* extanalysis = new ExtrapolationAnalysis();
-                        const int nvars=3;
-                        int vars[nvars] = {VxEnum, VyEnum, ThicknessEnum};
-                        for(int iv=0;iv<nvars;iv++){
-                                femmodel->parameters->SetParam(vars[iv],ExtrapolationVariableEnum); 
-                                extanalysis->Core(femmodel);
-                        }
-                        delete extanalysis;     
-
-                        /* solve level set equation */
-                        analysis = new LevelsetAnalysis();
-                        analysis->Core(femmodel);
-                        delete analysis;
-
-                        /* update vertices included for next calculation */
-                        GetMaskOfIceVerticesLSMx(femmodel);
-
-                        /* add computation domain mask to outputs */
-                        int outputs[1] = {IceMaskNodeActivationEnum};
-                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
-                }
-
+                /* from here on, prepare geometry for next time step*/
                 if(issmb)smb_core(femmodel);
 
@@ -155 +129 @@
                 }
 
-                /*Calculate new Basal melting on Floating ice*/
+                /*Calculate new basal melting on floating ice*/
                 FloatingiceMeltingRatex(femmodel);