Automatic difference of complex constitutive law

New to Julia
, ,
1

Hello, these days I am striving to realize a Jacobin of constitutive law in finite element simulation. I use the SparseDiffTools.jl, which accurately works at the linear elastic stage. However, when the displacements of nodes (finite element method) reaches the softening stage of the constitutive, I always encounter the error as below. I also marked the place where the error occurs in my function code. Any ideas are appreciated.

ERROR: UndefVarError: effectivestress not defined
Stacktrace:
 [1] constitutive42(nodalforce1::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, noddis::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}})
   @ Main c:\Users\DOJ14\OneDrive - University of Pittsburgh\PITT\Code\Julia\test_Mechanical\55 
nonlinear LDPM_explicit working now 2\2\2.27 nonlinear LDPM_explicit working copy 4.jl:561      
 [2] forwarddiff_color_jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::typeof(constitutive42), x::Vector{Float64}, jac_cache::ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}})
   @ SparseDiffTools C:\Users\DOJ14\.julia\packages\SparseDiffTools\zGdIo\src\differentiation\compute_jacobian_ad.jl:377
 [3] forwarddiff_color_jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::Function, x::Vector{Float64}; dx::Vector{Float64}, colorvec::Vector{Int64}, sparsity::SparseArrays.SparseMatrixCSC{Float64, Int64})
   @ SparseDiffTools C:\Users\DOJ14\.julia\packages\SparseDiffTools\zGdIo\src\differentiation\compute_jacobian_ad.jl:322
 [4] top-level scope
   @ c:\Users\DOJ14\OneDrive - University of Pittsburgh\PITT\Code\Julia\test_Mechanical\55 nonlinear LDPM_explicit working now 2\2\2.27 nonlinear LDPM_explicit working copy 4.jl:1025

My constitutive function is

function constitutive42(nodalforce1, noddis)

    #get volumetric strain for every element
    V_epsi = zeros(length1)
    division = zeros(length1)

    newverticecoor = verticeoftetrahedron1 + [noddis[(jj*6-5):(jj*6-3)] for jj = 1:length_nod][tttt1]
    for tnumber = 1:length2
        epsiVforte = (PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(newverticecoor[tnumber, 1]), PVector(newverticecoor[tnumber, 2]), PVector(newverticecoor[tnumber, 3]), PVector(newverticecoor[tnumber, 4]))) -
                      PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))) / PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))

        for numele = 1:length1
            if elements41[numele][1] ∈ tetrahedron1[tnumber, :] && elements41[numele][2] ∈ tetrahedron1[tnumber, :]
                V_epsi = V_epsi + vcat(vec(zeros(numele - 1, 1)), epsiVforte, vec(zeros(length1 - numele, 1)))
                division = division + vcat(vec(zeros(numele - 1, 1)), 1, vec(zeros(length1 - numele, 1)))
            end
        end
    end
    V_epsi_1 = deepcopy(@. V_epsi / division)
    #get other strains for every element
    for jii = 1:length_nod*6
        nodalforce1[jii] = 0.0
    end
    for kjh in element_num_inter

        epsiNk = BNjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BNik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiMk = BMjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BMik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiLk = BLjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BLik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiTk = sqrt(epsiMk^2 + epsiLk^2)
        ##epsiNk
        if epsiNk >= 0
            effectiveepsi = sqrt(epsiTk^2 * alpha + epsiNk^2)
            # tensile-shear boundary
            omig = atan(epsiNk / sqrt(alpha) / epsiTk)
            s = sin(omig)
            c = cos(omig)
            if omig == atan(sqrt(alpha) / mu)
                effecstress0b = 0.5 * (sigmt + 2 * sigma) * sigmt / ((sigmt + sigma) * sin(omig))
            else
                effecstress0b = (-(sigmt + sigma) * s + sqrt(((sigmt + sigma) * s)^2 + (alpha * (c / mu)^2 - s^2) * (sigmt + 2 * sigma) * sigmt)) / (alpha * (c / mu)^2 - s^2)
            end
            Hs = 2 * alpha * EN / ((2 * alpha * EN * Gs) / (sigms^2 * le2[kjh]) - 1)
            Ht = 2 * EN / ((2 * EN * Gt) / (sigmt^2 * le2[kjh]) - 1)
            H = Hs + (Ht - Hs) * (2 * omig / pi)^nt
            if maximum(effectivestrain_his[:, kjh]) == effectiveepsi
            effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([effectiveepsi - effecstress0b / EN, 0]))
            else
            effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([maximum(effectivestrain_his[:, kjh]) - effecstress0b / EN, 0]))
            end
   
            if mode[t_step::Int64, kjh] == "continuous load"
                effectivestress = minimum([effectiveepsi * EN, effectstressb])
            elseif mode[t_step::Int64, kjh] == "unloading"
                effectivestress = minimum([maximum([effectivestress_unlo[t_step::Int64, kjh] - (effectivestrain_unlo[t_step::Int64, kjh] - effectiveepsi) * EN, 0+3*effectiveepsi]), effectstressb])
            elseif mode[t_step::Int64, kjh] == "reloading"
                effectivestress = minimum([maximum([(effectiveepsi - kt * (effectivestrain_unlo[t_step::Int64, kjh] - effectivestress_unlo[t_step::Int64, kjh] / EN)) * EN, 0+3*effectiveepsi]), effectstressb])
            end


           
                 #***!!!!This is exactly where the error happens!***
                sigmNk = effectivestress / effectiveepsi * epsiNk
            

           
                sigmTk = effectivestress / effectiveepsi * epsiTk * alpha
            

            
                sigmLk = sigmTk * epsiLk / epsiTk
            
            
                sigmMk = sigmTk * epsiMk / epsiTk

        else #epsiNk < 0
            #compressive

            # rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
            # H = Hc / (1 + kapha2 * maximum([rdv-kapha1,0]))
            # sigmcb = sigmc0 * exp(-H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))
            

            # if mode_compre[t_step::Int64, kjh] == "continuous load"
            rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
            H = Hc / (1 + kapha2 * maximum([rdv-kapha1,0]))
            sigmcb = sigmc0 * exp(-H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))
            

            if mode_compre[t_step::Int64, kjh] == "continuous load"
                sigmNk = maximum([epsiNk * EN, sigmcb])
            else
                sigmNk = maximum([minimum([comprestress_unlo[t_step::Int64, kjh] - (comprestrain_unlo[t_step::Int64, kjh] - epsiNk) * Ed, sigmc0 - (comprestrain_unlo[t_step::Int64, kjh] - (comprestress_unlo[t_step::Int64, kjh] - sigmc0) / Ed - epsiNk) * EN, 0+3*epsiNk]), sigmcb])
            end
            # else
            #     sigmNk = maximum([minimum([comprestress_unlo[t_step::Int64, kjh] - (comprestrain_unlo[t_step::Int64, kjh] - epsiNk) * Ed, sigmc0 - (comprestrain_unlo[t_step::Int64, kjh] - (comprestress_unlo[t_step::Int64, kjh] - sigmc0) / Ed - epsiNk) * EN, 0]), sigmcb])
            # end##epsiNk = -0.0019340612033711502
            
            #shear stress under compression
      
            if sigmNk <= sigmN0
                if sigmNk> sigmc0
                sigmTb = sqrt((1-(sigmNk-sigmN0)^2/(sigmN0-sigmc0)^2)*sigmT0^2)
                else
                sigmTb =0.0
                end
            else
                sigmTb = mu * sqrt((sigmNk - sigmt - sigma)^2 - sigma^2)
            end
         
            if t_step::Int64 == 1
                sigmTk_f = epsiTk * ET
            elseif fricstrain_his[t_step::Int64-1, kjh] == 300.0
                sigmTk_f = shearsigm_his[t_step::Int64-1, kjh] + (epsiTk - sqrt(shearMstrain_his[t_step::Int64-1, kjh]^2 + shearLstrain_his[t_step::Int64-1, kjh]^2)) * ET
            else
                sigmTk_f = fricsigm_his[t_step::Int64-1, kjh] + (epsiTk - fricstrain_his[t_step::Int64-1, kjh]) * ET
            end
            if sigmTk_f > 0
                sigmTk = minimum([sigmTk_f, sigmTb+3*epsiTk])
            else
                sigmTk = maximum([sigmTk_f, -sigmTb+3*epsiTk])
                sigmTk = -sigmTk
            end

                if fricsigm_his[t_step::Int64-1, kjh] * fricsigm_his[t_step::Int64, kjh] < 0

                    epsiMk_1 = epsiMk - fricMstrain_his[t_step::Int64-1, kjh]# - (0 - fricsigm_his[t_step::Int64-1, kjh])/ET / sqrt(fricMstrain_his[t_step::Int64-1, kjh]^2 + fricLstrain_his[t_step::Int64-1, kjh]^2) * fricMstrain_his[t_step::Int64-1, kjh]
                    epsiLk_1 = epsiLk - fricLstrain_his[t_step::Int64-1, kjh]# - (0 - fricsigm_his[t_step::Int64-1, kjh])/ET / sqrt(fricMstrain_his[t_step::Int64-1, kjh]^2 + fricLstrain_his[t_step::Int64-1, kjh]^2) * fricLstrain_his[t_step::Int64-1, kjh]
                else
                    epsiMk_1 = fricMstrain_reset_his[t_step::Int64-1, kjh] + epsiMk - fricMstrain_his[t_step::Int64-1, kjh]
                    epsiLk_1 = fricLstrain_reset_his[t_step::Int64-1, kjh] + epsiLk - fricLstrain_his[t_step::Int64-1, kjh]
                end


            
                sigmLk = sigmTk * epsiLk_1 / sqrt(epsiMk_1^2 + epsiLk_1^2) ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk
            
           
                sigmMk = sigmTk * epsiMk_1 / sqrt(epsiMk_1^2 + epsiLk_1^2)
    

        end
        #println([sigmNk, sigmMk, sigmLk])
        # if isnan(sigmNk) || isnan(sigmMk) || isnan(sigmLk)
        #     println(kjh)
        # end
        # nodalforce1 = nodalforce1.*0.0 + vcat(vec(zeros(6*elements41[kjh][1]-6,1)), (-le2[kjh] * Ak2[kjh] * (sigmNk * BNik2[kjh] + sigmMk * BMik2[kjh] + sigmLk * BLik2[kjh]))',vec(zeros(length_nod*6 - (6*elements41[kjh][1]),1)))
        # nodalforce1 = nodalforce1.*0.0 +[vec(zeros(1,6*elements41[kjh][2]-6))' (le2[kjh] * Ak2[kjh] * (sigmNk * BNjk2[kjh] + sigmMk * BMjk2[kjh] + sigmLk * BLjk2[kjh])) vec(zeros(1,length_nod*6 - (6*elements41[kjh][2])))']'
        nodalforce1[(elements41[kjh][1]*6-5):(elements41[kjh][1]*6)] += (-le2[kjh] * Ak2[kjh] * (sigmNk * BNik2[kjh] + sigmMk * BMik2[kjh] + sigmLk * BLik2[kjh]))'
        nodalforce1[(elements41[kjh][2]*6-5):(elements41[kjh][2]*6)] += (le2[kjh] * Ak2[kjh] * (sigmNk * BNjk2[kjh] + sigmMk * BMjk2[kjh] + sigmLk * BLjk2[kjh]))'
  
    end
    nothing

end
2

I suspect none of the conditions in the if-else block matched. To test this, you can rewrite this as

effectivestress = if mode[t_step::Int64, kjh] == "continuous load"
    minimum([effectiveepsi * EN, effectstressb])
elseif mode[t_step::Int64, kjh] == "unloading"
    minimum([maximum([effectivestress_unlo[t_step::Int64, kjh] - (effectivestrain_unlo[t_step::Int64, kjh] - effectiveepsi) * EN, 0+3*effectiveepsi]), effectstressb])
elseif mode[t_step::Int64, kjh] == "reloading"
    minimum([maximum([(effectiveepsi - kt * (effectivestrain_unlo[t_step::Int64, kjh] - effectivestress_unlo[t_step::Int64, kjh] / EN)) * EN, 0+3*effectiveepsi]), effectstressb])
end

if isnothing(effectivestress)
    @show mode[t_step::Int64, kjh]
end

Then, if none of the conditions match effectivestress will be nothing and you can see what unexpected value is present.

3

You are right! I just forget to check it. The funny thins is when I face this yesterday I checked and find the bug. Thank you so much!