Changeset 17503

Timestamp:

03/20/14 11:02:44 (11 years ago)

Author:

jbondzio

Message:

CHG: Levelsetfunction initialization routine to solve penalization scheme in a transient fashion with artificial time steps

Location:

issm/trunk-jpl/src/c

Files:

• analyses/LsfReinitializationAnalysis.cpp (modified) (20 diffs)
• analyses/LsfReinitializationAnalysis.h (modified) (1 diff)
• modules/ModelProcessorx/ModelProcessorx.cpp (modified) (1 diff)

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

@@ -4 +4 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../cores/cores.h"
 
 #include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
@@ -29 +31 @@
                 }
         }
+
+        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
         int finiteelement=P1Enum;
+        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
         ::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
+        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -46 +52 @@
         /*parameters: */
         bool save_results;
+        int maxiter = 3;
+        int step;
+        IssmDouble reltol = 0.05;
+
+        Vector<IssmDouble>* lsfg     = NULL;
+        Vector<IssmDouble>* lsfg_old = NULL;
+
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
@@ -53 +66 @@
         /* set spcs for reinitialization */
         if(VerboseSolution()) _printf0_("Update spcs for reinitialization:\n");
-        UpdateReinitSPCs(femmodel);
-
-        if(VerboseSolution()) _printf0_("call computational core for reinitialization:\n");
-//      solutionsequence_lsfreinit_linear(femmodel);
+        SetReinitSPCs(femmodel);
+
+        step = 1;
+        for(;;){
+
+                _printf_("smoothing lsf slope\n");
+                /* smoothen slope of lsf for computation of normal on ice domain*/
+                levelsetfunctionslope_core(femmodel);
+
+                //solve current artificial time step
+                if(VerboseSolution()) _printf0_("call computational core for reinitialization in step " << step << ":\n");
+                solutionsequence_linear(femmodel);
+                GetSolutionFromInputsx(&lsfg,femmodel);
+
+                if(step>1){
+                        if(VerboseSolution()) _printf0_("   checking reinitialization convergence\n");
+                        if(ReinitConvergence(lsfg,lsfg_old,reltol)) break;
+                }
+                if(step>maxiter){
+                        if(VerboseSolution()) _printf0_("   maximum number reinitialization iterations " << maxiter << " reached\n");
+                        break;
+                }
+
+                /*update results and increase counter*/
+                delete lsfg_old;lsfg_old=lsfg;
+                step++;
+        }
 
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                int outputs = MaskIceLevelsetEnum;
-                femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+                int outputs[1] = {MaskIceLevelsetEnum};
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
         }
 
@@ -73 +109 @@
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
-
+        
         /*Intermediaries */
         const int dim = 2;
         int        i,row,col,stabilization;
-        IssmDouble Jdet,D_scalar,h;
-        IssmDouble dlsf[3],normal[3];
-        IssmDouble norm_dlsf;
-        IssmDouble hx,hy,hz,kappa;
+        IssmDouble Jdet,D_scalar;
+        IssmDouble dtau = 1.;
+        IssmDouble mu = 1.;
         IssmDouble* xyz_list = NULL;
 
@@ -88 +123 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix();
-        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble     D[dim][dim];
-
-        /*Retrieve all inputs and parameters*/
-        Input* lsfpicard_input=element->GetInput(LevelsetfunctionPicardEnum); _assert_(lsfpicard_input);
+        IssmDouble*    D                = xNew<IssmDouble>(dim*dim);
+        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
+        IssmDouble*    normal= xNew<IssmDouble>(dim);
+
         element->GetVerticesCoordinates(&xyz_list);
-        h = element->CharacteristicLength();
 
         /* Start  looping on the number of gaussian points: */
@@ -103 +137 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,element,xyz_list,gauss);
+                D_scalar=gauss->weight*Jdet;
+
+                if(dtau!=0.){
+                        element->NodalFunctions(basis,gauss);
+                        TripleMultiply(basis,numnodes,1,0,
+                                                &D_scalar,1,1,0,
+                                                basis,1,numnodes,0,
+                                                &Ke->values[0],1);
+                        D_scalar*=dtau;
+                }
+
                 GetBprime(Bprime,element,xyz_list,gauss);
-
-                /* Get normal from last iteration on lsf */
-                lsfpicard_input->GetInputDerivativeValue(&dlsf[0],xyz_list,gauss);
-
-                norm_dlsf=0.;
-                for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i]; 
-                norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
-                for(i=0;i<dim;i++)
-                        normal[i]=dlsf[i]/norm_dlsf;
-                
-                D_scalar=gauss->weight*Jdet;
 
                 for(row=0;row<dim;row++)
                         for(col=0;col<dim;col++)
                                 if(row==col)
-                                        D[row][col]=D_scalar*normal[row];
+                                        D[row*dim+col]=mu*D_scalar;
                                 else
-                                        D[row][col]=0.;
-                TripleMultiply(B,dim,numnodes,1,
-                                        &D[0][0],dim,dim,0,
+                                        D[row*dim+col]=0.;
+                TripleMultiply(Bprime,dim,numnodes,1,
+                                        D,dim,dim,0,
                                         Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
 
-                /* Stabilization *//*{{{*/
-                stabilization=1;
+                /* Stabilization */
+                stabilization=0;
                 if (stabilization==0){/* no stabilization, do nothing*/}
-                else if(stabilization==1){
-                        /* Artificial Diffusion */
-                        element->ElementSizes(&hx,&hy,&hz);
-                        h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
-                        kappa=h/2.; 
-                        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);
-                }
-                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][col]=h/(2.*1.)*normal[row]*normal[col];
-
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                &D[0][0],dim,dim,0,
-                                                Bprime,dim,numnodes,0,
-                                                &Ke->values[0],1);
-                }/*}}}*/
+                
         }/*}}}*/
 
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(D);
         delete gauss;
         return Ke;
@@ -168 +178 @@
         
         /*Intermediaries */
-        IssmDouble Jdet;
-        IssmDouble* xyz_list = NULL;
-        
+        int i,k;
+        int dim = 2;
+        IssmDouble dtau = 1.;
+        IssmDouble mu = 1.;
+        IssmDouble Jdet, D_scalar;
+        IssmDouble lsf;
+        IssmDouble norm_dlsf;
+        IssmDouble dbasis_normal;
+
         /*Fetch number of nodes */
         int numnodes = element->GetNumberOfNodes();
 
+        IssmDouble* xyz_list = NULL;
         IssmDouble* basis = xNew<IssmDouble>(numnodes);
+        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
+        IssmDouble* dlsf = xNew<IssmDouble>(dim);
+        IssmDouble* normal = xNew<IssmDouble>(dim);
         element->GetVerticesCoordinates(&xyz_list);
 
         /*Initialize Element vector*/
         ElementVector* pe = element->NewElementVector();
+
+        /*Retrieve all inputs and parameters*/
+        Input* lsf_input = element->GetInput(MaskIceLevelsetEnum); _assert_(lsf_input);
+        Input* lsf_slopex_input=element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+        Input* lsf_slopey_input=element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
 
         Gauss* gauss=element->NewGauss(2);
@@ -186 +211 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
-
-                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*basis[i]; 
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+                D_scalar=Jdet*gauss->weight;
+
+                if(dtau!=0.){
+                        /* old function value */
+                        lsf_input->GetInputValue(&lsf,gauss);
+                        for(i=0;i<numnodes;i++) pe->values[i]+=D_scalar*lsf*basis[i];
+                        D_scalar*=dtau;
+                }
+
+                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+                lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+
+                /*get normal*/
+                norm_dlsf=0.;
+                for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i]; 
+                norm_dlsf=sqrt(norm_dlsf);
+                if(norm_dlsf>0.)
+                        for(i=0;i<dim;i++)      normal[i]=dlsf[i]/norm_dlsf;
+                else
+                        for(i=0;i<dim;i++)      normal[i]=0.;
+
+                /* multiply normal and dbasis */
+                for(i=0;i<numnodes;i++){
+                        dbasis_normal=0.;
+                        for(k=0;k<dim;k++) dbasis_normal+=dbasis[k*numnodes+i]*normal[k];
+                        pe->values[i]+=D_scalar*mu*dbasis_normal; 
+                }
         }
 
         xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(dlsf);
+        xDelete<IssmDouble>(normal);
         return pe;
         }/*}}}*/
 void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        IssmDouble   lsf;
+        int          meshtype,dim,approximation,dofpernode;
+        int*         doflist = NULL;
+
+        /*Get some parameters*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum: dim = 2; dofpernode = 1; break;
+                case Mesh2DverticalEnum:   dim = 2; dofpernode = 1; break;
+                case Mesh3DEnum:           dim = 3; dofpernode = 1; break;
+                case Mesh3DtetrasEnum:     dim = 3; dofpernode = 1; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int numdof   = numnodes*dofpernode;
+
+        /*Fetch dof list and allocate solution vector*/
+        element->GetDofList(&doflist,approximation,GsetEnum);
+        IssmDouble* values = xNew<IssmDouble>(numdof);
+
+        /*Get inputs*/
+        Input* lsf_input=element->GetInput(MaskIceLevelsetEnum); _assert_(lsf_input);
+
+        Gauss* gauss=element->NewGauss();
+        for(int i=0;i<numnodes;i++){
+                gauss->GaussNode(element->FiniteElement(),i);
+
+                lsf_input->GetInputValue(&lsf,gauss);
+                values[i*dofpernode+0]=lsf;
+        }
+
+        solution->SetValues(numdof,doflist,values,INS_VAL);
+
+        /*Free ressources:*/
+        delete gauss;
+        xDelete<IssmDouble>(values);
+        xDelete<int>(doflist);
+
 }/*}}}*/
 void LsfReinitializationAnalysis::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 LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-        _error_("not implemented yet");
+        /* Do nothing for now */
 }/*}}}*/
 void LsfReinitializationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
@@ -260 +366 @@
 
 /* Other */
-void LsfReinitializationAnalysis::UpdateReinitSPCs(FemModel* femmodel){/*{{{*/
+void LsfReinitializationAnalysis::SetReinitSPCs(FemModel* femmodel){/*{{{*/
 
         int i,k, numnodes;
@@ -270 +376 @@
                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 for(k=0;k<element->GetNumberOfNodes();k++){
-                                node=element->GetNode(k);
+                        node=element->GetNode(k);
+                        if(node->IsActive()){
                                 node->DofInFSet(0); 
+                        }
                 }
         }
@@ -287 +395 @@
                         for(k=0;k<numnodes;k++){
                                 node=element->GetNode(k);
-                                node->ApplyConstraint(1,lsf[k]);
+                                if(node->IsActive()){
+                                        node->ApplyConstraint(1,lsf[k]);
+                                }
                         }
                         xDelete<IssmDouble>(lsf);
@@ -298 +408 @@
         /* Intermediaries */
         int i,k;
+        IssmDouble dmaxp,dmaxm,val;
+        Element* element;
 
         /*Initialize vector with number of vertices*/
         int numvertices=femmodel->vertices->NumberOfVertices();
-        Element* element;
 
         Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
         GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
         
+        /* set distance on elements intersected by zero levelset */
         for(i=0;i<femmodel->elements->Size();i++){
                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
-                        for(k=0;k<element->GetNumberOfVertices();k++)
-                                vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL); 
-        }
-
-        for(i=0;i<femmodel->elements->Size();i++){
-                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
+                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
                         SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
+                }
+        }
+        vec_dist_zerolevelset->Assemble();
+
+        /* Get maximum distance to interface along vertices */
+        dmaxp=0.; dmaxm=0.;
+        for(i=0;i<numvertices;i++){
+                vec_dist_zerolevelset->GetValue(&val,i); 
+                if((val>0.) && (val>dmaxp))
+                         dmaxp=val;
+                else if((val<0.) && (val<dmaxm))
+                         dmaxm=val;
+        }
+        //wait until all values are computed
+
+        /* set all none intersected vertices to max/min distance */
+        for(i=0;i<numvertices;i++){
+                vec_dist_zerolevelset->GetValue(&val,i);
+                if(val==1.) //FIXME: improve check
+                        vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
+                else if(val==-1.)
+                        vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
         }
 
@@ -327 +454 @@
         delete vec_dist_zerolevelset;
         delete dist_zerolevelset;
+
 }/*}}}*/
 void LsfReinitializationAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
@@ -334 +462 @@
 
         /* Intermediaries */
-        const int dim=3;
+        int dim=3;
         int i,d;
+        IssmDouble dist,lsf_old;
+
         int numvertices=element->GetNumberOfVertices();
-        IssmDouble s0[dim], s1[dim], v[dim];
-        IssmDouble dist,lsf_old;
 
         IssmDouble* lsf = xNew<IssmDouble>(numvertices);
         IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
         IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
+        IssmDouble* s0 = xNew<IssmDouble>(dim);
+        IssmDouble* s1 = xNew<IssmDouble>(dim);
+        IssmDouble* v = xNew<IssmDouble>(dim);
         IssmDouble* xyz_list = NULL;
         IssmDouble* xyz_list_zero = NULL;
@@ -355 +486 @@
 
         element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
-        for(d=0;d<dim;d++){
-                s0[d]=xyz_list_zero[0+d];
-                s1[d]=xyz_list_zero[3+d];
+        for(i=0;i<dim;i++){
+                s0[i]=xyz_list_zero[0+i];
+                s1[i]=xyz_list_zero[3+i];
         }
 
@@ -363 +494 @@
         for(i=0;i<numvertices;i++){
                 for(d=0;d<dim;d++)
-                        v[d]=xyz_list[3*i+d];
+                        v[d]=xyz_list[dim*i+d];
                 dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
                 signed_dist[i]=sign_lsf[i]*dist;
@@ -371 +502 @@
         for(i=0;i<numvertices;i++){
                 vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-                if(xIsNan<IssmDouble>(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
+                /* initial lsf values are +-1. Update those fields or if distance to interface smaller.*/
+                if(fabs(lsf_old)==1. || fabs(signed_dist[i])<fabs(lsf_old))
                         vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
         }
@@ -378 +510 @@
         xDelete<IssmDouble>(sign_lsf);
         xDelete<IssmDouble>(signed_dist);
+        xDelete<IssmDouble>(s0);
+        xDelete<IssmDouble>(s1);
+        xDelete<IssmDouble>(v);
+
 }/*}}}*/
 IssmDouble LsfReinitializationAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
@@ -403 +539 @@
         return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
 }/*}}}*/
-
+bool LsfReinitializationAnalysis::ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol){/*{{{*/
+
+        /*Output*/
+        bool converged = true;
+
+        /*Intermediary*/
+        Vector<IssmDouble>* dlsfg    = NULL;
+        IssmDouble          norm_dlsf,norm_lsf;
+
+        /*compute norm(du)/norm(u)*/
+        dlsfg=lsfg_old->Duplicate(); lsfg_old->Copy(dlsfg); dlsfg->AYPX(lsfg,-1.0);
+        norm_dlsf=dlsfg->Norm(NORM_TWO); norm_lsf=lsfg_old->Norm(NORM_TWO);
+        if (xIsNan<IssmDouble>(norm_dlsf) || xIsNan<IssmDouble>(norm_lsf)) _error_("convergence criterion is NaN!");
+        if((norm_dlsf/norm_lsf)<reltol){
+                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" < "<<reltol*100<<" %\n");
+        }
+        else{ 
+                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" > "<<reltol*100<<" %\n");
+                converged=false;
+        }
+
+        /*Cleanup*/
+        delete dlsfg;
+
+        return converged;
+}/*}}}*/

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

@@ -33 +33 @@
         
         /* Other */
-        void UpdateReinitSPCs(FemModel* femmodel);
+        void SetReinitSPCs(FemModel* femmodel);
         void SetDistanceOnIntersectedElements(FemModel* femmodel);
         void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);
         IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
+        bool ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol);
 };
 #endif

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

@@ -76 +76 @@
                 if(solution_enum==TransientSolutionEnum && analysis_enum==LevelsetAnalysisEnum && islevelset==false) continue;
                 if(solution_enum==TransientSolutionEnum && analysis_enum==ExtrapolationAnalysisEnum && islevelset==false) continue;
+                if(solution_enum==TransientSolutionEnum && analysis_enum==LsfReinitializationAnalysisEnum && islevelset==false) continue;