Changeset 26699

Timestamp:

12/02/21 12:35:35 (3 years ago)

Author:

Mathieu Morlighem

Message:

NEW: done with convergence criterion

Location:

issm/trunk-jpl/src/jl/solve

Files:

• inputs.jl (modified) (1 diff)
• modules.jl (modified) (1 diff)
• solutionsequences.jl (modified) (2 diffs)
• toolkits.jl (modified) (1 diff)

issm/trunk-jpl/src/jl/solve/inputs.jl

@@ -112 +112 @@
 
         if input.interp==P0Enum
-                return input.element_values
+                return input.element_values[1]
         elseif input.interp==P1Enum
                 value = input.element_values[1]*gauss.coords1[i] +  input.element_values[2]*gauss.coords2[i] +  input.element_values[3]*gauss.coords3[i]

issm/trunk-jpl/src/jl/solve/modules.jl

@@ -158 +158 @@
         end
 
-        return ug
+        return uf
 
 end#}}}

issm/trunk-jpl/src/jl/solve/solutionsequences.jl

@@ -29 +29 @@
                 Mergesolutionfromftogx(ug, uf, ys, femmodel.nodes)
 
-                print(ug)
-                error("compare with ISSM...")
+                #Check for convergence
+                converged = convergence(Kff,pf,uf,old_uf,restol,reltol,abstol)
+                InputUpdateFromSolutionx(analysis,ug,femmodel)
 
                 #Increase count
@@ -40 +41 @@
         end
 
+        print("\n   total number of iterations: ",  count,  "\n")
+
         error("STOP")
 
 end# }}}
+function convergence(Kff::IssmMatrix, pf::IssmVector, uf::IssmVector, old_uf::IssmVector, restol::Float64, reltol::Float64, abstol::Float64)#{{{
+
+        print("   checking convergence\n");
+
+        #If solution vector is empty, return true
+        if(IsEmpty(uf))
+                return true
+        end
+
+        #Convergence criterion #1: force equilibrium (Mandatory)
+        #compute K[n]U[n-1] - F
+        KUold  = Duplicate(uf);    MatMult!(Kff,old_uf,KUold)
+        KUoldF = Duplicate(KUold); VecCopy!(KUold, KUoldF); AXPY!(KUoldF, -1.0, pf)
+        nKUoldF = Norm(KUoldF,2)
+        nF      = Norm(pf,2)
+        res = nKUoldF/nF
+        if ~isfinite(res)
+                println("norm nf = ", nF, " and norm kuold = ",nKUoldF)
+                error("mechanical equilibrium convergence criterion is not finite!")
+        end
+        if(res<restol)
+                print("   mechanical equilibrium convergence criterion ", res*100, " < ", restol*100, " %\n")
+                converged=true
+        else
+                print("   mechanical equilibrium convergence criterion ", res*100, " > ", restol*100, " %\n")
+                converged=false;
+        end
+
+        #Convergence criterion #2: norm(du)/norm(u)
+        if ~isnan(reltol)
+                duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf)
+                ndu = Norm(duf, 2); nu = Norm(old_uf, 2)
+                if ~isfinite(ndu) | ~isfinite(nu) 
+                        error("convergence criterion is not finite!")
+                end
+                if((ndu/nu)<reltol)
+                        print("   Convergence criterion: norm(du)/norm(u)      ", ndu/nu*100, " < ", reltol*100, " %\n")
+                else
+                        print("   Convergence criterion: norm(du)/norm(u)      ", ndu/nu*100, " > ", reltol*100, " %\n")
+                        converged=false;
+                end
+        end
+
+        #Convergence criterion #3: max(du)
+        if ~isnan(abstol)
+                duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf)
+                nduinf= Norm(duf, 3)
+                if ~isfinite(nduinf) 
+                        error("convergence criterion is not finite!")
+                end
+                if(nduinf<abstol)
+                        print("   Convergence criterion: max(du)               ", nduinf, " < ", abstol, "\n")
+                else
+                        print("   Convergence criterion: max(du)               ", nduinf, " > ", abstol, "\n")
+                        converged=false;
+                end
+        end
+
+        return converged
+
+end#}}}

issm/trunk-jpl/src/jl/solve/toolkits.jl

@@ -58 +58 @@
 
 end#}}}
+function IsEmpty(vector::IssmVector)#{{{
+
+        return GetSize(vector)==0
+
+end#}}}
+function Duplicate(vector::IssmVector)#{{{
+
+        #Copy data structure
+        M=GetSize(vector)
+        return IssmVector(M)
+
+end#}}}
+function VecCopy!(x::IssmVector,y::IssmVector)#{{{
+
+        y.vector = x.vector
+
+end#}}}
+function Norm(x::IssmVector,type::Int64)#{{{
+
+        norm = 0
+
+        if type==2
+                for i in 1:length(x.vector)
+                        norm += x.vector[i]^2
+                end
+                norm = sqrt(norm)
+        elseif type==3
+                #Infinite norm
+                for i in 1:length(x.vector)
+                        if(abs(x.vector[i])>norm) norm = abs(x.vector[i]) end
+                end
+        else
+                error("type ",type," not supported yet")
+        end
+
+        return norm
+
+end#}}}
 
 #Operations