adding a permutation field to the sparse matrix pattern to alleviate …

Project.toml

@@ -1,7 +1,7 @@
 name = "FiniteElementContainers"
 uuid = "d08262e4-672f-4e7f-a976-f2cea5767631"
 authors = ["Craig M. Hamel <[email protected]> and contributors"]
-version = "0.10.1"
+version = "0.10.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -27,16 +27,16 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 
 [extensions]
-AbaqusReaderExt = "AbaqusReader"
 AMDGPUExt = "AMDGPU"
+AbaqusReaderExt = "AbaqusReader"
 BlockArraysExt = "BlockArrays"
 CUDAExt = "CUDA"
 MPIExt = "MPI"
 PartitionedArraysExt = "PartitionedArrays"
 
 [compat]
-AbaqusReader = "0.2.7"
 AMDGPU = "2"
+AbaqusReader = "0.2.7"
 Adapt = "4"
 Aqua = "0.8"
 Atomix = "1"
@@ -60,7 +60,6 @@ julia = "1"
 
 [extras]
 AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"

examples/moose-tutorials-clone/ex01/script.jl

@@ -0,0 +1,76 @@
+using FiniteElementContainers
+using LinearAlgebra
+
+# physics code we actually need to implement
+struct Poisson{F <: Function} <: AbstractPhysics{1, 0, 0}
+    func::F
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u_q = unpack_field(∇u_q, 1)
+    R_q = ∇N_X * ∇u_q - N * physics.func(X_q, 0.0)
+    return JxW * R_q[:]
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::Poisson, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    K_q = ∇N_X * ∇N_X'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+f(_, _) = 0.0
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Poisson(f)
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, one_func)
+        DirichletBC(:u, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    # solver = NewtonSolver(DirectLinearSolver(asm))
+    # integrator = QuasiStaticIntegrator(solver)
+    # evolve!(integrator, p)
+
+    FiniteElementContainers.update_bc_values!(p)
+    residual_int = VectorIntegral(asm, residual)
+    stiffness_int = MatrixIntegral(asm, stiffness)
+
+    K = stiffness_int(U, p)
+    remove_fixed_dofs!(stiffness_int)
+
+    for n in 1:2
+        R = residual_int(U, p)
+        remove_fixed_dofs!(residual_int)
+        ΔU = -K \ R.data
+        U.data .+= ΔU
+        @show n norm(R) norm(ΔU)
+    end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u",), U)
+    close(pp)
+end
+
+simulate()

examples/moose-tutorials-clone/ex02/script.jl

@@ -0,0 +1,78 @@
+using FiniteElementContainers
+using LinearAlgebra
+using StaticArrays
+
+# physics code we actually need to implement
+struct Advection{N} <: AbstractPhysics{1, 0, 0}
+    v::SVector{N, Float64}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u_q = unpack_field(∇u_q, 1)
+    term = dot(∇u_q, physics.v)
+    R_q = ∇N_X * ∇u_q + term * N
+    return JxW * R_q[:]
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::Advection, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    term = ∇N_X * physics.v
+    K_q = ∇N_X * ∇N_X' + N * term'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Advection(SVector{3, Float64}(0., 0., 1.))
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, one_func)
+        DirichletBC(:u, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    # solver = NewtonSolver(DirectLinearSolver(asm))
+    # integrator = QuasiStaticIntegrator(solver)
+    # evolve!(integrator, p)
+
+    FiniteElementContainers.update_bc_values!(p)
+    residual_int = VectorIntegral(asm, residual)
+    stiffness_int = MatrixIntegral(asm, stiffness)
+
+    K = stiffness_int(U, p)
+    remove_fixed_dofs!(stiffness_int)
+
+    for n in 1:2
+        R = residual_int(U, p)
+        remove_fixed_dofs!(residual_int)
+        ΔU = -K \ R.data
+        U.data .+= ΔU
+        @show n norm(R) norm(ΔU)
+    end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u",), U)
+    close(pp)
+end
+
+simulate()

examples/moose-tutorials-clone/ex03/script.jl

@@ -0,0 +1,107 @@
+using FiniteElementContainers
+using ForwardDiff
+using LinearAlgebra
+using StaticArrays
+
+struct CoupledPhysics <: AbstractPhysics{2, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    ∇u = unpack_field(∇u_q, 1)
+    ∇v = unpack_field(∇u_q, 2)
+    term = dot(∇u, ∇v)
+    R_u = ∇N_X * ∇u + term * N
+    R_v = ∇N_X * ∇v
+    R = vcat(R_u', R_v')
+    return JxW * R[:]
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::CoupledPhysics, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    # interps = map_interpolants(interps, x_el)
+    # (; X_q, N, ∇N_X, JxW) = interps
+    # ∇u_q = interpolate_field_gradients(physics, interps, u_el)
+    # ∇u = unpack_field(∇u_q, 1)
+    # ∇v = unpack_field(∇u_q, 2)
+
+    # # term = ∇N_X * physics.v
+    # term_1 = ∇N_X * ∇v
+    # term_2 = ∇N_X * ∇v
+    # K_uu = ∇N_X * ∇N_X' + N * term_1'
+    # K_uv = N * term_2'
+    # K_vu = zeros(typeof(K_uv))
+    # K_vv = ∇N_X * ∇N_X'
+
+    # # K_q = ∇N_X * ∇N_X' + N * term'
+    # K = vcat(
+    #     hcat(K_uu, K_uv),
+    #     hcat(K_vu, K_vv)
+    # )
+    # return JxW * K
+    ForwardDiff.jacobian(
+        z -> FiniteElementContainers.residual(
+            physics, interps, x_el, t, dt, z, u_el_old, state_old_q, state_new_q, props_el
+        ), u_el
+    )
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+two_func(_, _) = 2.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("mug.e")
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = CoupledPhysics()
+    props = create_properties(physics)
+    u = FiniteElementContainers.GeneralFunction(
+        ScalarFunction(V, :u),
+        ScalarFunction(V, :v)
+    )
+
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, two_func)
+        DirichletBC(:u, :top, zero_func)
+        DirichletBC(:v, :bottom, one_func)
+        DirichletBC(:v, :top, zero_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(mesh, asm, physics, props; dirichlet_bcs = dbcs)
+
+    solver = NewtonSolver(DirectLinearSolver(asm))
+    integrator = QuasiStaticIntegrator(solver)
+    evolve!(integrator, p)
+
+    # FiniteElementContainers.update_bc_values!(p)
+    # residual_int = VectorIntegral(asm, residual)
+    # stiffness_int = MatrixIntegral(asm, stiffness)
+
+    # K = stiffness_int(U, p)
+    # remove_fixed_dofs!(stiffness_int)
+
+    # for n in 1:4
+    #     R = residual_int(U, p)
+    #     remove_fixed_dofs!(residual_int)
+    #     ΔU = -K \ R.data
+    #     U.data .+= ΔU
+    #     @show n norm(R) norm(ΔU)
+    # end
+
+    U = p.h1_field
+
+    pp = PostProcessor(mesh, "output.e", u)
+    write_times(pp, 1, 0.0)
+    write_field(pp, 1, ("u", "v"), U)
+    close(pp)
+end
+
+simulate()

examples/moose-tutorials-clone/ex06/script.jl

@@ -0,0 +1,69 @@
+using FiniteElementContainers
+using LinearAlgebra
+
+# physics code we actually need to implement
+struct Transient <: AbstractPhysics{1, 0, 0}
+end
+
+@inline function FiniteElementContainers.residual(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    u_q, ∇u_q = interpolate_field_values_and_gradients(physics, interps, u_el)
+    u_q_old = interpolate_field_values(physics, interps, u_el_old)
+    ∇u = unpack_field(∇u_q, 1)
+    dudt = 20. * (u_q[1] - u_q_old[1]) / dt
+    R = dudt * N + ∇N_X * ∇u
+    return JxW * R[:]
+end
+
+@inline function FiniteElementContainers.stiffness(
+    physics::Transient, interps, x_el, t, dt, u_el, u_el_old, state_old_q, state_new_q, props_el
+)
+    interps = map_interpolants(interps, x_el)
+    (; X_q, N, ∇N_X, JxW) = interps
+    K_q = (20. / dt) * N * N' + ∇N_X * ∇N_X'
+    return JxW * K_q
+end
+
+# define some helper funcs upfront
+one_func(_, _) = 1.0
+zero_func(_, _) = 0.0
+
+function simulate()
+    # load a mesh
+    mesh = UnstructuredMesh("cyl-tet.e")
+    times = TimeStepper(0., 75., 75)
+    V = FunctionSpace(mesh, H1Field, Lagrange)
+    physics = Transient()
+    props = create_properties(physics)
+    u = ScalarFunction(V, :u)
+    asm = SparseMatrixAssembler(u; use_condensed=true)
+    dbcs = [
+        DirichletBC(:u, :bottom, zero_func)
+        DirichletBC(:u, :top, one_func)
+    ]
+    U = create_field(asm)
+    p = create_parameters(
+        mesh, asm, physics, props; 
+        dirichlet_bcs = dbcs,
+        times = times
+    )
+    solver = NewtonSolver(DirectLinearSolver(asm))
+    integrator = QuasiStaticIntegrator(solver)
+    pp = PostProcessor(mesh, "output.e", u)
+
+    n = 1
+    while times.time_current[1] < 75.0
+        evolve!(integrator, p)
+        U = p.h1_field
+        write_times(pp, n, times.time_current[1])
+        write_field(pp, n, ("u",), U)
+        n = n + 1
+    end
+
+    close(pp)
+end
+
+simulate()