Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

update hyperelasticity docs to show how to use NonlinearSolve #910

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

update hyperelasticity docs to show how to use NonlinearSolve #910

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
522480c
update hyperelasticity docs
oscardssmith May 13, 2024
0803edc
address code review
oscardssmith May 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 50 additions & 8 deletions docs/src/literate-tutorials/hyperelasticity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,12 +305,9 @@ function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
end
end;

# Finally, we define a main function which sets up everything and then performs Newton
# iterations until convergence.

function solve()
reset_timer!()

# Creating the grid, boundary conditions, and meshes.
function create_grid()
## Generate a grid
N = 10
L = 1.0
Expand Down Expand Up @@ -346,7 +343,6 @@ function solve()
L/2 - z + (y-L/2)*sin(θ) + (z-L/2)*cos(θ)
))
end

dbcs = ConstraintHandler(dh)
## Add a homogeneous boundary condition on the "clamped" edge
dbc = Dirichlet(:u, getfaceset(grid, "right"), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3])
Expand All @@ -364,7 +360,12 @@ function solve()
getfaceset(grid, "front"),
getfaceset(grid, "back"),
)
return dh, cv, fv, mp, ΓN, dbcs
end

# Finally, we define a main function which performs Newton iterations until convergence.
function manual_solve()
dh, cv, fv, mp, ΓN, dbcs = create_grid()
## Pre-allocation of vectors for the solution and Newton increments
_ndofs = ndofs(dh)
un = zeros(_ndofs) # previous solution vector
Expand All @@ -381,6 +382,7 @@ function solve()
newton_itr = -1
NEWTON_TOL = 1e-8
NEWTON_MAXITER = 30

prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:")

while true; newton_itr += 1
Expand All @@ -406,6 +408,7 @@ function solve()
Δu .-= ΔΔu
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
end


## Save the solution
@timeit "export" begin
vtk_grid("hyperelasticity", dh) do vtkfile
Expand All @@ -417,14 +420,53 @@ function solve()
return u
end

# Alternatively, it is possible to solve the same system with NonlinearSolve by defining a few wrappers:
# Note that as written, this is missing some of the optimizations that you would want for a real deployment.
function nonlinearsolve_solve()
dh, cv, fv, mp, ΓN, dbcs = create_grid()

# Run the simulation
## Create sparse matrix and residual vector
K = create_sparsity_pattern(dh)
g = zeros(ndofs(dh))

function nl_assemble(g, u, p)
(;K, dh, cv, fv, mp, ΓN, dbcs) = p
apply!(u, dbcs)
assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
apply_zero!(K, g, dbcs)
g
end
function nl_assemble_jac(K, u, p)
(;g, dh, cv, fv, mp, ΓN, dbcs) = p
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the reason for passing ΓN etc through p instead of just directly using them from the scope above (since this is a closure closing over it)?

apply!(u, dbcs)
assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
apply_zero!(K, g, dbcs)
K
end
Comment on lines +432 to +445
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just to be sure that I totally understand this correctly, in NonlinearSolve.jl there is currently no support to evaluate both at the same time, right? If so, then we can separate the functions above into jacobian and residual only:

function assemble_element_jac!(ke, cell, cv, fv, mp, ue, ΓN)
    ## Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    b = Vec{3}((0.0, -0.5, 0.0)) # Body force
    tn = 0.1 # Traction (to be scaled with surface normal)
    ndofs = getnbasefunctions(cv)

    for qp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, qp)
        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇u = function_gradient(cv, qp, ue)
        F = one(∇u) + ∇u
        C = tdot(F) # F' ⋅ F
        ## Compute stress and tangent
        S, ∂S∂C = constitutive_driver(C, mp)
        P = F ⋅ S
        I = one(S)
        ∂P∂F =  otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I)

        ## Loop over test functions
        for i in 1:ndofs
            ## Test function and gradient
            δui = shape_value(cv, qp, i)
            ∇δui = shape_gradient(cv, qp, i)
            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
            for j in 1:ndofs
                ∇δuj = shape_gradient(cv, qp, j)
                ## Add contribution to the tangent
                ke[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
            end
        end
    end
end;

function assemble_element_residual!(ge, cell, cv, fv, mp, ue, ΓN)
    ## Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    b = Vec{3}((0.0, -0.5, 0.0)) # Body force
    tn = 0.1 # Traction (to be scaled with surface normal)
    ndofs = getnbasefunctions(cv)

    for qp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, qp)
        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇u = function_gradient(cv, qp, ue)
        F = one(∇u) + ∇u
        C = tdot(F) # F' ⋅ F
        ## Compute stress and tangent
        S, ∂S∂C = constitutive_driver(C, mp)
        P = F ⋅ S

        ## Loop over test functions
        for i in 1:ndofs
            ## Test function and gradient
            δui = shape_value(cv, qp, i)
            ∇δui = shape_gradient(cv, qp, i)
            ## Add contribution to the residual from this test function
            ge[i] += ( ∇δui ⊡ P - δui ⋅ b ) * dΩ
        end
    end

    ## Surface integral for the traction
    for face in 1:nfaces(cell)
        if (cellid(cell), face) in ΓN
            reinit!(fv, cell, face)
            for q_point in 1:getnquadpoints(fv)
                t = tn * getnormal(fv, q_point)
                dΓ = getdetJdV(fv, q_point)
                for i in 1:ndofs
                    δui = shape_value(fv, q_point, i)
                    ge[i] -= (δui ⋅ t) * dΓ
                end
            end
        end
    end
end;

+- typos and we can do the same with with assemble_global!

I will benchmark this after I finish the prioritized tasks.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't it be better to keep them in the same function with a flag to chose what to evaluate. That would make it easier to evaluate both for solvers that support it, to share some computations?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure about the performance hit associated with splitting up the Jacobian and residual computation, especially for element formulations which come with some kind of condensation for the Jacobian.


nlf = NonlinearFunction(nl_assemble, jac=nl_assemble_jac, jac_prototype=K)
p = (;K, g, dh, cv, fv, mp, ΓN, dbcs)
nlp = NonlinearProblem(nlf, u, p)
u = solve(nlp, NewtonRaphson(), abstol=1e-8, maxiter=30).u

## Save the solution
@timeit "export" begin
vtk_grid("hyperelasticity", dh) do vtkfile
vtk_point_data(vtkfile, dh, u)
end
end

u = solve();
print_timer(title = "Analysis with $(getncells(grid)) elements", linechars = :ascii)
return u
end

# Run the simulation
u = manual_solve();

## test the result #src
using Test #src
@test norm(u) ≈ 4.761404305083876 #src
@test u ≈ nonlinearsolve_solve() #src

#md # ## Plain program
#md #
Expand Down