Created
December 22, 2018 02:57
-
-
Save kailaix/809a6ea38d997a7d070a99608f3750d7 to your computer and use it in GitHub Desktop.
Assembler Programming
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # This is an example of solving the Poisson problem using the Assembler Programming method | |
| # (-Δ) u(x) = f(x) | |
| # u(x) = 0, x ∈ [0,1]^2 | |
| # The exact solution is u(x,y) = sin(pi*x)sin(pi*y) | |
| # and f(x,y) = 2pi^2 * sin(pi*x)sin(pi*y) | |
| using PyPlot | |
| import Plots | |
| using PyCall | |
| using SparseArrays | |
| @pyimport numpy as np | |
| N = 400 | |
| x = LinRange(0,1,N+1) | |
| X,Y = np.meshgrid(x,x) | |
| Ni = N+1 | |
| h = 1/N | |
| ii = Int64[] | |
| jj = Int64[] | |
| vv = Float64[] | |
| function add(i,j,k,l,v) | |
| if i<1 || i>Ni || j<1 || j>Ni || k<1 || k>Ni || l<1 || l>Ni | |
| return | |
| end | |
| idx = i + (j-1)*Ni | |
| idy = k + (l-1)*Ni | |
| push!(ii, idx) | |
| push!(jj, idy) | |
| push!(vv, v) | |
| end | |
| function assemble() | |
| return sparse(ii,jj,vv,Ni*Ni,Ni*Ni) | |
| end | |
| function add_derive(i,j,c,order) | |
| if order=="0" | |
| add(i,j,i,j,c) | |
| elseif order=="1" | |
| add(i,j,i+1,j,c/2h) | |
| add(i,j,i-1,j,-c/2h) | |
| elseif order=="2" | |
| add(i,j,i,j+1,c/2h) | |
| add(i,j,i,j-1,-c/2h) | |
| elseif order=="11" | |
| add(i,j,i+1,j,c/h^2) | |
| add(i,j,i-1,j,c/h^2) | |
| add(i,j,i,j,-2c/h^2) | |
| elseif order=="22" | |
| add(i,j,i,j+1,c/h^2) | |
| add(i,j,i,j-1,c/h^2) | |
| add(i,j,i,j,-2c/h^2) | |
| elseif order=="12" | |
| add(i,j,i-1,j+1,-c/4h^2) | |
| add(i,j,i-1,j-1,c/4h^2) | |
| add(i,j,i+1,j+1,c/4h^2) | |
| add(i,j,i+1,j-1,-c/4h^2) | |
| else | |
| error("not implemented yet") | |
| end | |
| end | |
| ################# example ################### | |
| rhs = zeros(Ni^2) | |
| I = ones(Bool, Ni^2) | |
| for i = 1:Ni | |
| for j = 1:Ni | |
| if i in [1,Ni] || j in [1,Ni] | |
| I[i + (j-1)*Ni] = false | |
| continue | |
| end | |
| add_derive(i,j,-1.0,"11") | |
| add_derive(i,j,-1.0,"22") | |
| rhs[i + (j-1)*Ni] = 2pi^2 * sin(pi*X[i,j])sin(pi*Y[i,j]) | |
| end | |
| end | |
| A = assemble() | |
| u = zeros(Ni^2) | |
| u[I] = A[I,I] \ rhs[I] | |
| uexact = vec( sin.(pi*X) .* sin.(pi*Y)) | |
| println(maximum(abs.(u-uexact))) | |
Author
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
We illustrate the idea of assembler programming idea here. In PDEs, we usually need to code the stencils or DOFs. The resulting matrix is usually sparse. The programmer only needs to focus on the local assembling, and the global assembling is taken care by the assembler.
This technique can also be seen as one derived from recursive programming.