Closed ChrisRackauckas closed 4 years ago
@YingboMa @HarrisonGrodin solidifying the discussion here.
Thinking about some ways to automatically handle boundary condition operator generation for polar circle domains down in DiffEqOperators.jl - the boundary condition on theta can simply be periodic, but the condition on r is a bit more difficult, The higher index end can be arbitrarily chosen, but it is the lower index end i'm concerned with here.
An important question is what the lower index end of r should correspond to: r=0 or r=dr. In the first case it would be necessary to ensure somehow that the lower index points of r for all theta contained the same value. Would it be sufficient to place in the ghost nodes along this boundary the value at r[2] at the diametrically opposed value of theta? In this case it would be necessary to have a uniformly sampled theta with an even number of points. In the second case, where r[1] = dr[1], the value at r=0 could be interpolated from the rest of r, and this value placed in the ghost nodes along this boundary - this may be preferable because then theta could be non uniform.
That's a question for DiffEqOperators. It has the information to do the discretization, and it has to find out how it wants to do that. However, this is all about making sure there's general enough information such that it can make the correct choice, and I believe that is the case.
I'm anxiously waiting that this gets finished. I practically only solve PDEs in my work, so I usually discretize them using some Collocation Method, transforming it into a DAE and solve it with the IDA solver. The problem is that the only linear solver that does not give me IC errors is the Dense solver, which is not appropriate at all for my problems. When this is done will I be able to solve my problems more efficiently?
Hopefully it'll set them up in a very efficient manner haha. That's the goal
That would be amazing! I actually wanted to use Modia to do process simulations but unfortunately they cannot handle, yet, sparse jacobians and event handling, things that I absolutely need. Looking forward to see this development!
https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/342 https://github.com/JuliaDiffEq/DiffEqOperators.jl/pull/179 https://github.com/JuliaDiffEq/ModelingToolkit.jl/pull/180
Makes this work:
using ModelingToolkit, DiffEqOperators, DiffEqBase, LinearAlgebra
# Define some variables
@parameters t x
@variables u(..)
@derivatives Dt'~t
@derivatives Dxx''~x
eq = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ - x * (x-1) * sin(x),
u(t,0) ~ 0, u(t,1) ~ 0]
domains = [t ∈ IntervalDomain(0.0,1.0),
x ∈ IntervalDomain(0.0,1.0)]
pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
discretization = MOLFiniteDifference(0.1)
prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent
using OrdinaryDiffEq
sol = solve(prob,Tsit5(),saveat=0.1)
using Plots
plot(prob.space,Array(prob.extrapolation*sol[1]))
plot!(prob.space,Array(prob.extrapolation*sol[2]))
plot!(prob.space,Array(prob.extrapolation*sol[3]))
plot!(prob.space,Array(prob.extrapolation*sol[4]))
savefig("pde_solve.png")
Implement a smart discretize function in DiffEqOperators.jl for automatic finite difference discretizations. This is where a ton of work will have to go, since it will need to parse the equation and find out what CenteredDifference and UpwindDifference operators to create. This work will continue after this closes.
This is now https://github.com/JuliaDiffEq/DiffEqOperators.jl/issues/180
Hi Chris I am trying to run the above equation to learn about the PDE solution using NN. I am running the code in JupyterLab I am getting the following error. Can you kindly help to understand the error and how to resolve it.
**MethodError: no method matching one(::Type{Array{Float64,1}}) Closest candidates are: one(!Matched::Type{Missing}) at missing.jl:103 one(!Matched::BitArray{2}) at bitarray.jl:422 one(!Matched::Missing) at missing.jl:100 ...
Stacktrace: [1] DirichletBC(::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\Administrator.julia\packages\DiffEqOperators\QiB3a\src\derivative_operators\BC_operators.jl:140
Best Regards Irfan
You might need to have to work at it to make this work again. That was only a preliminary prototype, not a robust result expected to be used by anyone yet. The issue will be closed when this is done.
There are still things to do.
Add another interesting test case: [ ] Implementation of a Saint Venant test case, see DOI:10.3384/ecp15119237.
Boundary conditions: [ ] Support Neumann boundary conditions [ ] Support Robin boundary conditions Note: a discussion about PDESystem boundary specification is taking place in SciML/ModelingToolkit.jl#526
Stationary cases and ODE cases: [ ] Auto-Discretization of stationary cases. E.g. 0 = d2u/dx2 + d2u/dy2 + f(x,y) [ ] Auto-Discretization of ODE cases, where the ODESystem is packed inside a PDESystem.
More dimensions: [ ] Auto-Discretization of 3D PDE problems. [ ] Auto-Discretization of 4 or more dimension PDE problems.
Some ToDo's were defined within src/MOL_discretization.jl. I will solve them in the following months, so it is not worth creating new Issues.
With https://nextjournal.com/kirill_zubov/physics-informed-neural-networks-pinns-solver-on-julia-gsoc-2020-final-report and https://github.com/SciML/DiffEqOperators.jl/pull/250, I think we can mark this, which is about establishing a common interface, as completed. We can have specific follow-up issues and discussions in the repos for the specific discretizations though. And we should probably find a nice package-independent way of documenting the PDE interface (maybe in MTK?)
Sorry for my question.
How is the number of spatial discretization points determined in a PDESystem
?
Thanks!
Depends on the method. Some don't even need to discretize.
Sir, I am interested to contribute to these issues but I am new to this. I am familiar with numerical differential equations, linear algebra, infinite series & PDE. Can you please guide me how to get started and contribute to these issues and could you assign me some issues? Thanks.
Hey, this work lives on in https://github.com/SciML/MethodOfLines.jl. See the issues and start making contributions. Join the Slack chat https://julialang.org/slack/ #diffeq-bridged for dev discussions
Building off of the story in https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260 , DiffEqOperators getting auto finite difference discretizations likely done this summer. Thus the next step is automating the PDE discretization from a high level description. Since ModelingToolkit.jl is well-developed, it makes sense to utilize that to give the high level description. Thus I put a prototype together:
Description of the code
The equation and boundary conditions are just given by ModelingToolkit.jl operations. The ConstrainedEquation allows for specifying equations only on specific parts. Then each independent variable needs a domain. Domains can be multi-dimensional, and variables can be multi-dimensional, so this lets us smartly define variables in domains (and in manifolds!) and equations using them.
The totality of the PDE is in the PDESystem, just a ModelingToolkit.AbstractSystem object. The workflow for solving a PDE is thus:
So (1) is all ModelingToolkit.jl. (2) is just a
discretize
function and some problem types. We will need to add LinearProblem and NonlinearProblem with some solve dispatches for this. After (3), we will want to wrap the solution so it can get a nicer plot recipe.Thus the steps are:
\
, GMRES (LinSolveGMRES()
),nlsolve
discretize
stub function in DiffEqBasediscretize
function in DiffEqOperators.jl for automatic finite difference discretizations. This is where a ton of work will have to go, since it will need to parse the equation and find out whatCenteredDifference
andUpwindDifference
operators to create. This work will continue after this closes.