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.

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

Support (f::Fun)(::Domain) for restricting a function to a sub domain #414

Open
BoundaryValueProblems opened this issue Nov 3, 2016 · 20 comments

Comments

@BoundaryValueProblems
Copy link

BoundaryValueProblems commented Nov 3, 2016

Hello. I want to solve the Helmholtz equation on the square [-1,1] x [-1,1] with the Dirichlet or Neumann boundary conditions where the boundary values are prescribed by the boundary value of a given 2D function. For example, let's use a Gabor function of the form:

julia> d = Interval()^2
julia> g=Fun((x,y)->exp(-(x^2+3*y^2)*cos(5*pi*x+7*pi*y),d)

Now, I want to use the boundary values of this function g to solve the Helmholtz equation. Following the example in README.md, I tried:

julia> Δ = Laplacian(d)                            # Represent the Laplacian
julia> f = g(∂(d))                                 # f should represent the boundary value of g
julia> u = [dirichlet(d);Δ+500I]\[f;0.]            # Solve the PDE

However, it generates the following error at the 2nd line. What should I do?
Thanks!
PS: By the way, there is a typo in the Helmholtz solver example in README.md.
The third line above in README.md says:

julia> u = [Dirichlet(d); Δ+500I]\[f;0.]            # Solve the PDE

But Dirichlet should be replaced by dirichlet; otherwise it causes an error.

Here is the error generated by:

julia> f = g(∂(d))                                 # f should represent the boundary value of g

MethodError: no method matching start(::ApproxFun.PiecewiseInterval{Complex{Float64}})
Closest candidates are:
start(!Matched::SimpleVector) at essentials.jl:170
start(!Matched::Base.MethodList) at reflection.jl:258
start(!Matched::IntSet) at intset.jl:184
...
in append_any(::ApproxFun.PiecewiseInterval{Complex{Float64}}, ::Vararg{ApproxFun.PiecewiseInterval{Complex{Float64}},N}) at ./essentials.jl:98
in evaluate(::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2},Float64}, ::ApproxFun.PiecewiseInterval{Complex{Float64}}) at /Users/xxx/.julia/v0.5/ApproxFun/src/Fun/Fun.jl:141
in (::ApproxFun.Fun{ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2},Float64})(::ApproxFun.PiecewiseInterval{Complex{Float64}}, ::Vararg{ApproxFun.PiecewiseInterval{Complex{Float64}},N}) at...

@dlfivefifty
Copy link
Member

The README is for the development version, which hasn't been tagged yet.

The syntax f = g(∂(d)) doesn't work (though I guess it should!). Use instead

Fun((x,y)->g(x,y),(d))

@dlfivefifty dlfivefifty changed the title How to specify the boundary condition from a specified 2D function Support (f::Fun)(::Domain) for restricting a function to a sub domain Nov 3, 2016
@BoundaryValueProblems
Copy link
Author

Thanks, dlfivefifty!
However, I tried the following two possibilities, and still got errors.

julia> f=Fun((x,y)->g(x,y),∂(d))
julia> f=Fun((x,y)->exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y),∂(d))

The error messages in both cases are:
v MethodError: no method matching (::##3#4)(::Complex{Float64})
Closest candidates are:
#3(::Any, !Matched::Any) at console:1
in zerocfsFun at ApproxFun/src/Fun/constructors.jl:134
in #Fun#45 at ApproxFun/src/Fun/constructors.jl:206
in at base/
in #Fun#45 at ApproxFun/src/Fun/constructors.jl:202
in #Fun#48 at ApproxFun/src/Fun/constructors.jl:213
in ApproxFun.Fun{S,T} at ApproxFun/src/Fun/constructors.jl:213

Thanks for your help!

@dlfivefifty
Copy link
Member

dlfivefifty commented Nov 3, 2016

What version of Julia are you on?

You can you also try on development. The tagged version of PDE solving is less robust. (There was a major revamp.)

Pkg.checkout(“ApproxFun”,”development”)
d=Interval()^2
f=Fun((x,y)->exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y),(d))
u = linsolve([Dirichlet(d); Δ+500I],[f;0.];tolerance=1E-7)

Note that it is necessary to specify the tolerance as otherwise the solution time will be too long (the solver is still being optimized.)

@dlfivefifty
Copy link
Member

If you require multiple RHS, you can do the following the make subsequent solves faster:

d=Interval()^2
f=Fun((x,y)->exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y),(d))
QR=qrfact([Dirichlet(d); Δ+500I])
@time u = linsolve(QR,[f;0.];tolerance=1E-7)   # 2.392173s
@time u = linsolve(QR,[f;0.];tolerance=1E-7)   # 0.009747s

@BoundaryValueProblems
Copy link
Author

The julia version I have is: 0.5.0.
As you suggested, I just did:

Pkg.checkout("ApproxFun", "development")
d=Interval()^2
f=Fun((x,y)->exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y),∂(d))

Then, I got the same error as I reported earlier.
Do I need to use the development version of julia, i.e., 0.6.x, in order for these to work?
If possible, I want to stick with 0.5.0 if possible.
Thanks!

@dlfivefifty
Copy link
Member

No this is on 0.6. 0.5 has some outstanding bugs with conflicting packages which I suspect you are running into. You would have to delete your .julia folder to fix it.

On 4 Nov. 2016, at 8:32 am, BoundaryValueProblems [email protected] wrote:

The julia version I have is: 0.5.0.
As you suggested, I just did:

Pkg.checkout("ApproxFun", "development")
d=Interval()^2
f=Fun((x,y)->exp(-(x^2+3_y^2))_cos(5_pi_x+7_pi_y),∂(d))
Then, I got the same error as I reported earlier.
Do I need to use the development version of julia, i.e., 0.6.x, in order for these to work?
If possible, I want to stick with 0.5.0 if possible.
Thanks!


You are receiving this because you commented.
Reply to this email directly, view it on GitHub #414 (comment), or mute the thread https://github.com/notifications/unsubscribe-auth/ABLDqfgiNHaw_UIkQFvJeBnU99drsleMks5q6lL8gaJpZM4Ko5JY.

@BoundaryValueProblems
Copy link
Author

Well, I just installed julia version 0.6.0-dev.1178, and did exactly what you suggested including installing the developmental version of ApproxFun.
Then, the following line generated tons of errors:

f=Fun((x,y)->exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y),∂(d))

ERROR: MethodError: Cannot convert an object of type Int64 to an object of type ApproxFun.Infinity{Bool}
This may have arisen from a call to the constructor ApproxFun.Infinity{Bool}(...),
since type constructors fall back to convert methods.
in mapfoldl_impl(::ApproxFun.#dimension, ::Base.#*, ::ApproxFun.Infinity{Bool}, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}, ::Int64) at ./reduce.jl:43
in mapfoldl(::ApproxFun.#dimension, ::Function, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}) at ./reduce.jl:73
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:207
in (::Core.#kw#Type)(::Array{Any,1}, ::Type{ApproxFun.Fun}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at ./:0
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:206
in #Fun#84(::Array{Any,1}, ::Type{T}, ::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217
in ApproxFun.Fun{S,T}(::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217

@dlfivefifty
Copy link
Member

0.6 is not supported

On 4 Nov. 2016, at 9:13 am, BoundaryValueProblems [email protected] wrote:

Well, I just installed julia version 0.6.0-dev.1178, and did exactly what you suggested including installing the developmental version of ApproxFun.
Then, the following line generated tons of errors:

f=Fun((x,y)->exp(-(x^2+3_y^2))_cos(5_pi_x+7_pi_y),∂(d))
ERROR: MethodError: Cannot convert an object of type Int64 to an object of type ApproxFun.Infinity{Bool}
This may have arisen from a call to the constructor ApproxFun.Infinity{Bool}(...),
since type constructors fall back to convert methods.
in mapfoldl_impl(::ApproxFun.#dimension, ::Base.#*, ::ApproxFun.Infinity{Bool}, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}, ::Int64) at ./reduce.jl:43
in mapfoldl(::ApproxFun.#dimension, ::Function, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}) at ./reduce.jl:73
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:207
in (::Core.#kw#Type)(::Array{Any,1}, ::Type{ApproxFun.Fun}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at ./:0
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:206
in #Fun#84(::Array{Any,1}, ::Type{T}, ::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217
in ApproxFun.Fun{S,T}(::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217


You are receiving this because you commented.
Reply to this email directly, view it on GitHub #414 (comment), or mute the thread https://github.com/notifications/unsubscribe-auth/ABLDqaTIuLjRDUMZkRRTKfcxNuQyRZxfks5q6lyWgaJpZM4Ko5JY.

@dlfivefifty
Copy link
Member

Sorry, I made a typo. I meant to say “No, this is on 0.5"

On 4 Nov. 2016, at 9:14 am, Sheehan Olver [email protected] wrote:

0.6 is not supported

On 4 Nov. 2016, at 9:13 am, BoundaryValueProblems <[email protected] mailto:[email protected]> wrote:

Well, I just installed julia version 0.6.0-dev.1178, and did exactly what you suggested including installing the developmental version of ApproxFun.
Then, the following line generated tons of errors:

f=Fun((x,y)->exp(-(x^2+3_y^2))_cos(5_pi_x+7_pi_y),∂(d))
ERROR: MethodError: Cannot convert an object of type Int64 to an object of type ApproxFun.Infinity{Bool}
This may have arisen from a call to the constructor ApproxFun.Infinity{Bool}(...),
since type constructors fall back to convert methods.
in mapfoldl_impl(::ApproxFun.#dimension, ::Base.#*, ::ApproxFun.Infinity{Bool}, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}, ::Int64) at ./reduce.jl:43
in mapfoldl(::ApproxFun.#dimension, ::Function, ::Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}}) at ./reduce.jl:73
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:207
in (::Core.#kw#Type)(::Array{Any,1}, ::Type{ApproxFun.Fun}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at ./:0
in #Fun#81(::String, ::Type{T}, ::Function, ::ApproxFun.TensorSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Interval{Float64}},ApproxFun.Chebyshev{ApproxFun.Interval{Float64}}},ApproxFun.RealBasis,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:206
in #Fun#84(::Array{Any,1}, ::Type{T}, ::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xxx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217
in ApproxFun.Fun{S,T}(::Function, ::ApproxFun.ProductDomain{Tuple{ApproxFun.Interval{Float64},ApproxFun.Interval{Float64}},Float64,2}) at /Users/xx/.julia/v0.6/ApproxFun/src/Fun/constructors.jl:217


You are receiving this because you commented.
Reply to this email directly, view it on GitHub #414 (comment), or mute the thread https://github.com/notifications/unsubscribe-auth/ABLDqaTIuLjRDUMZkRRTKfcxNuQyRZxfks5q6lyWgaJpZM4Ko5JY.

@dlfivefifty
Copy link
Member

OK, I remembered in the last tagged version it uses complex numbers instead of x, y. So try this:

f=Fun(z->(x=real(z);y=imag(z);exp(-(x^2+3*y^2))*cos(5*pi*x+7*pi*y)),(d))
Δ = Laplacian(d) 
u = pdesolve([dirichlet(d);Δ+500I],[f;0.],200)

@dlfivefifty
Copy link
Member

Or, on development in Julia 0.5, you can try moving your .julia folder aside.

@dlfivefifty
Copy link
Member

Sorry for the typo....really annoying one.

@BoundaryValueProblems
Copy link
Author

Well, before I installed 0.6, I tried the developmental version of ApproxFun as you suggested, and as I reported earlier, I got the error.
Please clarify what setting I should use:
(1) Julia 0.5.0; tagged version of ApproxFun.jl; with complex numbers as you just suggested
or
(2) Julia 0.5.0; development version of ApproxFun.jl; without using complex numbers as you suggested.
In either case, should I rename my ~/.julia and rebuild a new ~/.julia?
I would greatly appreciate if you could clarify these once and for all.
Thanks!

@dlfivefifty
Copy link
Member

Yes. In case (1) you might be able to get away with not moving away .julia.

There's been several people who've reported weird bugs where moving away .julia fixes it. I believe it's due to the open issue JuliaLang/julia#18465 but that's just a guess.

@dlfivefifty
Copy link
Member

Note that (1) and (2) have very different solvers. The solver in (1) is faster if you are changing k. The solver in (2) is faster if yo are changing the RHS.

The solver in (1) has since been moved to https://github.com/ApproxFun/PDESchurFactorization.jl but has not been fixed up for the current development version.

@dlfivefifty
Copy link
Member

Ah, for (2) you might have to also do Pkg.checkout("BandedMatrices"). I'll tag that now to remove one step.

@BoundaryValueProblems
Copy link
Author

Thanks a lot!
I tried (1) on OS X and (2) on Windows, and both worked!
Please let me know when you plan to tag the current development version of ApproxFun so that I do not need to use two different ways to specify the boundary values from a given 2D function.
Thank you very much for your effort of creating extremely useful package!!

@dlfivefifty
Copy link
Member

I'll try to tag it by next week. I was hoping to speed up the solution to Fourier PDEs first, but I think that will have to be postponed.

@BoundaryValueProblems
Copy link
Author

Thanks! When you commit the tagged version next week, what I should do with the development version I just installed on OS X via Pkg.checkout("ApproxFun", "development")? Can I just do Pkg.update() at that point and don't need to do anything else? Or do I have to remove the development version of ApproxFun first, then upgrade it via Pkg.update()?

@dlfivefifty
Copy link
Member

help?> Pkg.free
  free(pkg)

  Free the package pkg to be managed by the package manager again. It calls
  Pkg.resolve() to determine optimal package versions after. This is an
  inverse for both Pkg.checkout and Pkg.pin.

  You can also supply an iterable collection of package names, e.g.,
  Pkg.free(("Pkg1", "Pkg2")) to free multiple packages at once.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants