Merge branch 'main' into l1-gs-smoother

Author: Abdelrahman912

Project.toml

@@ -10,6 +10,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -44,6 +45,8 @@ BlockArrays = "1"
 DiffEqBase = "6"
 FastBroadcast = "0.3.5"
 Ferrite = "1"
+ForwardDiff = "0.10.38"
+JET = "0.9"
 LinearSolve = "2"
 Logging = "1.10"
 ModelingToolkit = "9"

bak/examples/Manifest.toml

@@ -63,8 +63,6 @@ odeform = semidiscretize(
 u₀ = zeros(Float32, OS.function_size(odeform))
 steady_state_initializer!(u₀, odeform)
 
-# io = ParaViewWriter("spiral-wave-test")
-
 timestepper = #Thunderbolt.ReactionTangentController(
     OS.LieTrotterGodunov((
         BackwardEulerSolver(
@@ -87,19 +85,20 @@ integrator = OS.init(problem, timestepper, dt=dt₀, verbose=true)
 TimerOutputs.enable_debug_timings(Thunderbolt)
 step!(integrator) # precompile for benchmark below
 
+io = ParaViewWriter("Niederer-Test")
 TimerOutputs.reset_timer!()
 for (u, t) in OS.TimeChoiceIterator(integrator, tspan[1]:dtvis:tspan[2])
-    # dh = odeform.functions[1].dh
-    # φ = u[odeform.dof_ranges[1]]
-    # @info t,norm(u)
+    dh = odeform.functions[1].dh
+    φ = u[odeform.solution_indices[1]]
+    @info t,norm(u)
     # sflat = ....?
-    # store_timestep!(io, t, dh.grid) do file
-    #     Thunderbolt.store_timestep_field!(file, t, dh, φ, :φₘ)
-    #     # s = reshape(sflat, (Thunderbolt.num_states(ionic_model),length(φ)))
-    #     # for sidx in 1:Thunderbolt.num_states(ionic_model)
-    #     #    Thunderbolt.store_timestep_field!(io, t, dh, s[sidx,:], state_symbol(ionic_model, sidx))
-    #     # end
-    # end
+    store_timestep!(io, t, dh.grid) do file
+        Thunderbolt.store_timestep_field!(file, t, dh, φ, :φₘ)
+        # s = reshape(sflat, (Thunderbolt.num_states(ionic_model),length(φ)))
+        # for sidx in 1:Thunderbolt.num_states(ionic_model)
+        #    Thunderbolt.store_timestep_field!(io, t, dh, s[sidx,:], state_symbol(ionic_model, sidx))
+        # end
+    end
 end
 TimerOutputs.print_timer()
 # TimerOutputs.disable_debug_timings(Thunderbolt)

bak/examples/conduction-velocity-benchmark.jl

@@ -0,0 +1,58 @@
+using Thunderbolt
+using Thunderbolt.TimerOutputs
+
+TimerOutputs.enable_debug_timings(Thunderbolt)
+
+mesh = generate_mesh(Hexahedron, (4,4,4))
+material = Thunderbolt.LinearMaxwellMaterial(
+    E₀ = 70e3,
+    E₁ = 20e3,
+    μ  = 1e3,
+    η₁ = 1e3,
+    ν  = 0.3,
+)
+tspan = (0.0,1.0)
+Δt = 0.1
+
+# Clamp three sides
+dbcs = [
+    Dirichlet(:d, getfacetset(mesh, "left"), (x,t) -> (0.0, 0.0, 0.0), [1,2,3]),
+    Dirichlet(:d, getfacetset(mesh, "right"), (x,t) -> (0.1, 0.0, 0.0), [1,2,3]),
+]
+
+quasistaticform = semidiscretize(
+    QuasiStaticModel(:d, material, ()),
+    FiniteElementDiscretization(
+        Dict(:d => (LagrangeCollection{1}()^3 => QuadratureRuleCollection(2))),
+        dbcs,
+    ),
+    mesh
+)
+problem = QuasiStaticProblem(quasistaticform, tspan)
+
+# Create sparse matrix and residual vector
+timestepper = BackwardEulerSolver(;
+    inner_solver=Thunderbolt.MultiLevelNewtonRaphsonSolver(;
+        # global_newton=NewtonRaphsonSolver(),
+        # local_newton=NewtonRaphsonSolver(),
+    )
+)
+integrator = init(problem, timestepper, dt=Δt, verbose=true)
+
+step!(integrator) # Precompile
+
+TimerOutputs.reset_timer!()
+solve!(integrator)
+TimerOutputs.print_timer()
+
+cc = first(CellIterator(mesh.grid))
+F = one(Tensor{2,3})
+coefficient = nothing
+qr = QuadratureRule{RefHexahedron}(2)
+element_cache = Thunderbolt.setup_element_cache(integrator.cache.stage.nlsolver.global_solver_cache.op.integrator.volume_model, qr, integrator.cache.stage.nlsolver.global_solver_cache.op.dh.subdofhandlers[1])
+state_cache   = element_cache.internal_cache
+qp = QuadraturePoint(1, qr.points[1])
+
+using BenchmarkTools
+@btime Thunderbolt.solve_local_constraint(F, coefficient, material, state_cache, cc, qp, 0.0)
+@btime Thunderbolt._query_local_state(state_cache, cc, qp)