BrianChung10 / SOP

2 stars 1 forks source link

Problem with the Allen-Cahn equation #12

Open fz920 opened 2 years ago

fz920 commented 2 years ago

Hi @dlfivefifty @ioannisPApapadopoulos, I am now working with the Allen-Cahn equation and trying to find the solution using the finite difference method in 2D. But I am unsure how to construct the residual function effectively, and my construction of residual F causes some problems when applying the jacobian function. Could you have a look at it? Thanks!

The codes are in the file Finite_Difference_2D / allen.jl.

fz920 commented 2 years ago

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}) Closest candidates are: (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/rounding.jl:200 (::Type{T})(::T) where T<:Number at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/boot.jl:770 (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/char.jl:50 ... Stacktrace: [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}) @ Base ./number.jl:7 [2] setindex!(::Matrix{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}, ::Int64, ::Int64) @ Base ./array.jl:905 [3] F(u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}}) @ Main ~/Desktop/group_project_deflation/SOP/Finite_Differences_2D/allen.jl:43 [4] chunk_mode_jacobian(f::typeof(F), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}}}) @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:227 [5] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}}}, ::Val{true}) @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:23 [6] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(F), Float64}, Float64, 12}}}) (repeats 2 times) @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:19 [7] newton(f::typeof(F), x0::Vector{Float64}, max_iter::Int64, eps::Float64) @ Main ~/Desktop/group_project_deflation/SOP/Finite_Differences_2D/allen.jl:58 [8] newton(f::Function, x0::Vector{Float64}) @ Main ~/Desktop/group_project_deflation/SOP/Finite_Differences_2D/allen.jl:52 [9] top-level scope @ ~/Desktop/group_project_deflation/SOP/Finite_Differences_2D/allen.jl:70

dlfivefifty commented 2 years ago

The error is saying you are trying to set an entry in a Matrix{Float64} to a Dual. The issue is here:

https://github.com/BrianChung10/SOP/blob/db2323263793699a594db8ae8a37db870e9b8d95/Finite_Differences_2D/allen.jl#L12

where you are creating a Matrix{Float64}. You should take in T and create a A = zeros(T,n+1, n+1)... hopefully that fixes it

fz920 commented 2 years ago

Thank you for your suggestion. That indeed solves the problem. Now I have an error called the LAPACKException(1) error which I believe is caused by the fact the Jacobian of my function F is not invertible. I print out the Jacobian matrix and find that most entries are 0 which makes the Jacobian non-invertible. Do you think there is any problem with my function F in Finite_Difference_2D / allen.jl? Thanks.

ioannisPApapadopoulos commented 2 years ago

probably. A few things I noticed

  1. In your definition of the residual F, I think it will always spit out a vector that always has some entries of +1 or -1. Therefore the norm of F will never go to zero and your Newton method will never converge.
  2. There are as many zero rows as there are discretization points on the boundary. Since these are fixed to be constant, differentiating should give you a zero row as expected. You need to pick the correct entries to sub in to make the Jacobian invertible (just like in your previous examples for ODEs/PDEs).