v0.8.1

evolve! for time dependent solution

Project.toml

@@ -1,7 +1,7 @@
 name = "VoronoiFVM"
 uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
 authors = ["Juergen Fuhrmann <[email protected]>"]
-version = "0.8.0"
+version = "0.8.1"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
@@ -19,7 +19,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 DiffResults = "^1.0.0"
 DocStringExtensions = "^0.8.0"
-ExtendableGrids = "^0.2.0"
+ExtendableGrids = "^0.2.2"
 ExtendableSparse = "^0.2.6"
 ForwardDiff = "^0.10.0"
 IterativeSolvers = "^0.8.0"

docs/src/changes.md

@@ -1,6 +1,9 @@
 # Changes
 
-## v0.8 Apr 28, 2019
+## v0.8.1 May 2 2020
+- Introduce evolve! : time solver with automatic timestep control
+
+## v0.8 Apr 28, 2020
 - Replaced VoronoiFVM grid module by  [ExtendableGrids.jl](https://github.com/j-fu/ExtendableGrids.jl)
 - Moved grid generation, modification, plotting  over to ExtendableGrids
 - Necessary changes in codes using VoronoiFVM:
@@ -12,7 +15,7 @@
   - For using any methods on grids like cellmask! one nees to use `ExtendableGrids`
   - Subgrids now are of the same type `ExtendableGrids`,  views are currently defined for vectors only.
 
-## v0.7 Feb 28 2019
+## v0.7 Feb 28 2020
 - API modification:
   - __Breaking__:
     - `data` parameter passed to physics callbacks only if `Physics` object is created with `data` parameter.
@@ -36,7 +39,7 @@
   - Return `nothing` from mutating methods to avoid some allocations 
   - Indexing in `formfactors.jl` with `Int`
 
-## v0.6.5 Jan 25 2019
+## v0.6.5 Jan 25 2020
 - use updateindex! for matrix, depend on ExtendableSparse 0.2.6
 
 ## v0.6.4 2020-01-20

src/VoronoiFVM.jl

@@ -1,3 +1,8 @@
+"""
+$(README)
+
+$(EXPORTS)
+"""
 module VoronoiFVM
 
 
@@ -72,6 +77,7 @@ export getdof
 export setdof!
 export value
 export solve!
+export evolve!
 export embed!
 export freqdomain_impedance
 export measurement_derivative

src/vfvm_newtoncontrol.jl

@@ -85,6 +85,11 @@ mutable struct NewtonControl
     """
     edge_cutoff::Float64
 
+    Δt::Float64
+    Δt_max::Float64
+    Δt_min::Float64
+    Δt_grow::Float64
+    Δu_opt::Float64
 
     NewtonControl()=NewtonControl(new())
 end
@@ -109,6 +114,11 @@ function NewtonControl(this)
     this.Δp=1
     this.Δp_max=1
     this.Δp_min=1.0e-3
+    this.Δt=0.1
+    this.Δt_max=1
+    this.Δt_min=1.0e-3
+    this.Δt_grow=1.2
+    this.Δu_opt=0.1
     this.handle_exceptions=false
     this.edge_cutoff=0.0
 

src/vfvm_solver.jl

@@ -723,6 +723,81 @@ function embed!(
     return solution
 end
 
+################################################################
+"""
+$(TYPEDSIGNATURES)
+
+Time dependent solver for abstract system
+
+Use iimplicit Euler method with + damped   Newton's method  to 
+solve time dependent problem.
+"""
+function evolve!(
+    solution::AbstractMatrix{Tv}, # Solution
+    xinival::AbstractMatrix{Tv},   # Initial value 
+    this::AbstractSystem{Tv},    # Finite volume system
+    times::AbstractVector;
+    control=NewtonControl(),      # Newton solver control information
+    pre=function(sol,t) end,       # Function for preparing step
+    post=function(sol,oldsol, t, Δt) end,      # Function for postprocessing successful step
+    delta=(u,v)->norm(u-v,Inf)
+) where Tv
+    inival=copy(xinival)
+    Δt=control.Δt
+    if control.verbose
+        @printf("  Evolution: start\n")
+    end
+    for i=1:length(times)-1
+        Δt=max(Δt,control.Δt_min)
+        tstart=times[i]
+        tend=times[i+1]
+        t=tstart
+        istep=0
+
+        while t<tend
+            solved=false
+            t0=t
+            while !solved
+                solved=true
+                try
+                    t=t0+Δt
+                    pre(solution,t)
+                    solve!(solution,inival, this ,control=control,tstep=Δt)
+                catch err
+                    if (control.handle_exceptions)
+                        _print_error(err,stacktrace(catch_backtrace()))
+                    else
+                        rethrow(err)
+                    end
+                    # here we reduce the embedding step and retry the solution
+                    Δt=Δt*0.5
+                    if tp<control.tp_min
+                        throw(EmbeddingError())
+                    end
+                    solved=false
+                    if control.verbose
+                        @printf("  Evolution: retry: Δt=%.3e\n",Δt)
+                    end
+                end
+            end
+            istep=istep+1
+            Δu=delta(solution, inival)
+            if control.verbose
+                @printf("  Evolution: istep=%d t=%.3e Δu=%.3e\n",istep,t,Δu)
+            end
+            post(solution,inival,t, Δt)
+            inival.=solution
+            if t<tend
+                Δt=min(control.Δt_max,Δt*control.Δt_grow,Δt*control.Δu_opt/(Δu+1.0e-14),tend-t)
+            end
+        end
+    end
+    if control.verbose
+        @printf("  Evolution: success\n")
+    end
+    return solution
+end
+
 
 
 ################################################################