I mean in solving Numerical PDE the matrix is usually quiet big but sparse. the normal “/” seems to use LU decomposition which makes a saprse matrix become two dense and will outofmemory qiuckly.
If you have a sparse matrix, then \ uses UMFPACK for LU decomposition, so it should not make it dense. Do you have an example of this happening?
2 Likes
This is the code.Most part is building a sparse matrix. It is solving a three dimensional PDE and N is the how much part we divide the PDE in different dimension. When N=30 it warks well, but when N=100
Stacktrace:
[1] umfpack_numeric!(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:258
[2] #lu#1(::Bool, ::typeof(lu), ::SparseMatrixCSC{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:160
[3] lu at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:154 [inlined]
[4] \(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SparseArrays\src\linalg.jl:1374
[5] pdesolver(::Int64, ::Float64, ::Float64, ::Float64, ::Array{Float64,2}, ::Int64) at C:\Users\lyu\Documents\DC\dachuang-again\threeDimensionalProblem.jl:349
[6] pdesolver() at C:\Users\lyu\Documents\DC\dachuang-again\threeDimensionalProblem.jl:13
[7] top-level scope at REPL[7]:1
using SparseArrays, PyPlot,NPZ
using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onehot, crossentropy, throttle,Params
using Base.Iterators: repeated
using LinearAlgebra
using Sobol
using BSON:@save
using BSON:@load
using FastGaussQuadrature,SparseGrids
## K and L is used to define the schostic layer
function pdesolver(σ=1,d=pi/2,L_y=2π,L_z=2π,ϵ=zeros(3,3),λ=-1)
N=30
K=size(ϵ,1)
L=size(ϵ,2)
dx=1/N;dy=1/N;dz=1/N
x=0:dx:1;y=0:dx:1;z=0:dx:1
Nx=length(x);Ny=length(y);Nz=length(z)
# define the stochastic layer here
h(rr1,rr2,sy,sz) = ( sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
h_y(rr1,rr2,sy,sz) = (0 + 2*rr1*pi* cos(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
h_yy(rr1,rr2,sy,sz) = (0 - (2*rr2*pi)^2* sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
h_z(rr1,rr2,sy,sz) = (0 + 2*rr2*pi* sin(2*rr1*pi*sy)*cos(2*rr2*pi*sz));
h_zz(rr1,rr2,sy,sz) = (0 - (2*rr2*pi)^2* sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
##The main finite difference code for u, u1, and u2
EPh =0;EU= zeros(Nx,Ny,Nz);
Datah = zeros(Ny,Nz); Datah_y = zeros(Ny,Nz); Datah_yy = zeros(Ny,Nz); Datah_z = zeros(Ny,Nz); Datah_zz = zeros(Ny,Nz);
sumindex=0 # size: Nz*1
for rr1=1:K
for rr2=1:L
multiplier=float(rr1)^(λ)*float(rr2)^(λ)
Datah += h.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
Datah_y += h_y.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
Datah_yy += h_yy.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
Datah_z += h_z.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
Datah_zz += h_zz.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
sumindex += multiplier
end
end
Datah=Datah/sumindex
Datah_y=Datah_y/sumindex
Datah_yy=Datah_yy/sumindex
Datah_z=Datah_z/sumindex
Datah_zz=Datah_zz/sumindex
## conputing, loop over each experimental result
## form the data structure and index for stiffness matrix and load vector
dx2=dx*dx; dy2=dy*dy; dz2=dz*dz; d4xy=4*dx*dy; d4xz=4*dx*dz; σ2=σ*σ;
DOFsp=(Nx-3)*Ny*Nz*15+Ny*Nz*10*2;
DOFu=(Nx-1)*Ny*Nz;
spI4Au=zeros(DOFsp,1); spJ4Au=zeros(DOFsp,1); spS4Au=zeros(DOFsp,1);
spI4bu=zeros(DOFu,1); spJ4bu=ones(DOFu,1);
spS4bu = zeros(DOFu,1); spS4bu1 = zeros(DOFu,1); spS4bu2 = zeros(DOFu,1);
## define the weights for trapezoidal rule
trapz = 2*ones(Nz,1); trapy = 2*ones(1,Ny); trapx= 2*ones(1,Nx);trapwt=zeros(Nx,Ny,Nz) # define the weights for trapezoidal rule
trapz[1] = 1; trapz[Nz] = 1;
trapy[1] = 1; trapy[Ny] = 1;
trapx[1] = 1; trapx[Nx] = 1;
for i=1:Nx
for j=1:Ny
for k=1:Nz
trapwt[i,j,k]=trapx[i]*trapy[j]*trapz[k];
end
end
end
wt=reshape(trapwt,DOFu+Ny*Nz,1);
##range of index for u. For x, 2<=i<=N_x; For y, 1<=j<=N_y; For z, 1<=k<=N_z
index=0; index4b=0;
for i=3:Nx-1
for j=2:Ny-1
for k=2:Nz-1
spI4Au[index+1:index+15,1]= ones(15)*((i -2)*Ny*Nz+(j -1)*Nz+k);
spJ4Au[index+1:index+15,1]= [ (i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(j+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1 ,1]= (i -2)*Ny*Nz+(j -1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
end
for j=1:Ny #front grids,i=2
for k=1:Nz
spI4Au[index+1:index+10,1]=ones(10)*((2-2)*Ny*Nz+(j-1)*Nz+k);
if j==1
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
elseif j==Ny
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
else
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
end
spI4bu[index4b+1,1]=(2-2)*Ny*Nz+(j-1)*Nz+k;
index=index+10; index4b=index4b+1;
end
end
for j=1:Ny #back grids
for k=1:Nz
if j==1
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
elseif j==Ny
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
else
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
end
spI4Au[index+1:index+10,1]=ones(10)*((Nx-2)*Ny*Nz+(j-1)*Nz+k);
spI4bu[index4b+1,1]=(Nx-2)*Ny*Nz+(j-1)*Nz+k;
index=index+10; index4b=index4b+1;
end
end
for i=3:Nx-1 #rightest grids
for j=1:Ny #k=Nz
if j==1
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
elseif j==Ny
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(1+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
else
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
end
spI4Au[index+1:index+15]=((i-2)*Ny*Nz+(j-1)*Nz+Nz)*ones(15);
spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+Nz;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #leftest grids
for j=1:Ny #k=1
if j==1
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(j+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
elseif j==Ny
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(1+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(1+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(1+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
else
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(j+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
end
spI4Au[index+1:index+15]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+1);
spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+1;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=1
j=1;
for k=2:Nz-1
spI4Au[index+1:index+15,1]=((i-2)*Ny*Nz+(1-1)*Nz+k)*ones(15);
spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(j+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=Ny
j=Ny;
for k=2:Nz-1
spI4Au[index+1:index+15,1]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+k);
spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(1+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(1+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(1+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
DataG=zeros(Nx,Ny,Nz)
for i=1:Nx
for j=1:Ny
for k=1:Nz
DataG[i,j,k]=-exp(-(d-Datah[j,k])*(1-x[i]))#-exp(-((d-Datah[j,k])*(1-x[i])).^2)
end
end
end
#-2sin(d*x[i])*(sin(L_y*y[j])+1)*(sin(L_z*z[k])+1)-sin(d*x[i])*sin(L_y*y[j])*(sin(L_z*z[k])+1)-sin(d*x[i])*(sin(L_y*y[j])+1)*sin(L_z*z[k])
## assmble and solve for u
index=0;index4b=0;
for i=3:Nx-1
for j=2:Ny-1
for k=2:Nz-1
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-2σ2*(1/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-2σ2*(1/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy; A/dx2+F/(2*dx); D/d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
end
i=2;
for j=1:Ny #front grids,i=2
for k=1:Nz
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-2σ2*(1/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-2σ2*(1L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+10,1]=[B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy;A/dx2+F/(2*dx); D/d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+10; index4b=index4b+1;
end
end
i=Nx;
for j=1:Ny #back grids
for k=1:Nz
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+10,1]=[0; 2*A/dx2; 0;0; 0;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+10; index4b=index4b+1;
end
end
k=Nz;
for i=3:Nx-1 #rightest grids
for j=1:Ny #k=Nz
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy; A/dx2+F/(2*dx); D/d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #leftest grids
k=1;
for j=1:Ny #k=1
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy; A/dx2+F/(2*dx); D/d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=1
j=1;
for k=2:Nz-1
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy;A/dx2+F/(2*dx); D/4*d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=Ny
j=Ny;
for k=2:Nz-1
A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
B=σ2/L_y^2;
C=σ2/L_z^2;
D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
-D/d4xy; A/dx2+F/(2*dx); D/d4xy; -E/d4xz; E/d4xz;];
spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
Au =sparse(vec(spI4Au),vec(spJ4Au),vec(spS4Au),DOFu,DOFu);
bu = sparse(vec(spI4bu),vec(spJ4bu),vec(spS4bu),DOFu,1);
u=Au\Matrix(bu)
u= [zeros(Ny*Nz,1); u];
U=reshape(Matrix(u),Nx,Ny,Nz);
Phtmp = (1/8)*sum(sum(sum((d*ones(Nx,Nx^2)-repeat(Datah,1,Nx)).*reshape(U,Nx,Nx^2).*reshape(trapwt,Nx,Nx^2)*dx*dy*dz)))
return Phtmp
end
A minimal example would be kind on your helpers. Have a read through Please read: make it easier to help you, in particular point 4.
3 Likes
I understand your meaning, but the problem is building a sparse matrix is the main part of the code. In order to perserve the feature of this matrix, I cannot change most part of the code. Only thing I can do is change the value in each index into a rand number like under, but it may also cause some SingularException(0) problems because of the / sometimes.
Actually, the only useful line is the last five line u = Au\bu, when the size of the matric increase it will Outofmemory
function pdesolver(σ=1,d=pi/2,L_y=2π,L_z=2π,ϵ=zeros(3,3),λ=-1)
N=20
dx=1/N;dy=1/N;dz=1/N
x=0:dx:1;y=0:dx:1;z=0:dx:1
Nx=length(x);Ny=length(y);Nz=length(z)
##The main finite difference code for u, u1, and u2
EPh =0;EU= zeros(Nx,Ny,Nz);
Datah = zeros(Ny,Nz); Datah_y = zeros(Ny,Nz); Datah_yy = zeros(Ny,Nz); Datah_z = zeros(Ny,Nz); Datah_zz = zeros(Ny,Nz);
## conputing, loop over each experimental result
## form the data structure and index for stiffness matrix and load vector
dx2=dx*dx; dy2=dy*dy; dz2=dz*dz; d4xy=4*dx*dy; d4xz=4*dx*dz; σ2=σ*σ;
DOFsp=(Nx-3)*Ny*Nz*15+Ny*Nz*10*2;
DOFu=(Nx-1)*Ny*Nz;
spI4Au=zeros(DOFsp,1); spJ4Au=zeros(DOFsp,1); spS4Au=zeros(DOFsp,1);
spI4bu=zeros(DOFu,1); spJ4bu=ones(DOFu,1);
spS4bu = zeros(DOFu,1); spS4bu1 = zeros(DOFu,1); spS4bu2 = zeros(DOFu,1);
##range of index for u. For x, 2<=i<=N_x; For y, 1<=j<=N_y; For z, 1<=k<=N_z
index=0; index4b=0;
for i=3:Nx-1
for j=2:Ny-1
for k=2:Nz-1
spI4Au[index+1:index+15,1]= ones(15)*((i -2)*Ny*Nz+(j -1)*Nz+k);
spJ4Au[index+1:index+15,1]= [ (i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(j+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1 ,1]= (i -2)*Ny*Nz+(j -1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
end
for j=1:Ny #front grids,i=2
for k=1:Nz
spI4Au[index+1:index+10,1]=ones(10)*((2-2)*Ny*Nz+(j-1)*Nz+k);
if j==1
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
elseif j==Ny
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(2-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
else
if k==1
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+2;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(2 -2)*Ny*Nz+(j-1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k; (2 -2)*Nz*Ny+(j+1-1)*Nz+k; (2 -2)*Ny*Nz+(j-1)*Nz+k-1; (2 -2)*Ny*Nz+(j-1)*Nz+k+1;
(2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
end
end
spI4bu[index4b+1,1]=(2-2)*Ny*Nz+(j-1)*Nz+k;
index=index+10; index4b=index4b+1;
end
end
for j=1:Ny #back grids
for k=1:Nz
if j==1
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
elseif j==Ny
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(2-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
else
if k==1
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
elseif k==Nz
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+2;];
else
spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(Nx -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k; (Nx -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx -2)*Ny*Nz+(j-1)*Nz+k+1;];
end
end
spI4Au[index+1:index+10,1]=ones(10)*((Nx-2)*Ny*Nz+(j-1)*Nz+k);
spI4bu[index4b+1,1]=(Nx-2)*Ny*Nz+(j-1)*Nz+k;
index=index+10; index4b=index4b+1;
end
end
for i=3:Nx-1 #rightest grids
for j=1:Ny #k=Nz
if j==1
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
elseif j==Ny
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(1+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
else
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz; (i -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
end
spI4Au[index+1:index+15]=((i-2)*Ny*Nz+(j-1)*Nz+Nz)*ones(15);
spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+Nz;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #leftest grids
for j=1:Ny #k=1
if j==1
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(j+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
elseif j==Ny
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(1+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(1+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(1+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
else
spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+1; (i -2)*Nz*Ny+(j+1-1)*Nz+1; (i -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i -2)*Ny*Nz+(j-1)*Nz+1+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
end
spI4Au[index+1:index+15]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+1);
spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+1;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=1
j=1;
for k=2:Nz-1
spI4Au[index+1:index+15,1]=((i-2)*Ny*Nz+(1-1)*Nz+k)*ones(15);
spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(Ny-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(j+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=Ny
j=Ny;
for k=2:Nz-1
spI4Au[index+1:index+15,1]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+k);
spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(1+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
(i -2)*Ny*Nz+(j-1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k; (i -2)*Nz*Ny+(1+1-1)*Nz+k; (i -2)*Ny*Nz+(j-1)*Nz+k-1; (i -2)*Ny*Nz+(j-1)*Nz+k+1;
(i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(1+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
index=index+15; index4b=index4b+1;
end
end
DataG=zeros(Nx,Ny,Nz)
for i=1:Nx
for j=1:Ny
for k=1:Nz
DataG[i,j,k]=-exp(-(d-Datah[j,k])*(1-x[i]))#-exp(-((d-Datah[j,k])*(1-x[i])).^2)
end
end
end
#-2sin(d*x[i])*(sin(L_y*y[j])+1)*(sin(L_z*z[k])+1)-sin(d*x[i])*sin(L_y*y[j])*(sin(L_z*z[k])+1)-sin(d*x[i])*(sin(L_y*y[j])+1)*sin(L_z*z[k])
## assmble and solve for u
index=0;index4b=0;
for i=3:Nx-1
for j=2:Ny-1
for k=2:Nz-1
spS4Au[index+1:index+15,1]=rand(15);
spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
end
i=2;
for j=1:Ny #front grids,i=2
for k=1:Nz
spS4Au[index+1:index+10,1]=rand(10);
spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+10; index4b=index4b+1;
end
end
i=Nx;
for j=1:Ny #back grids
for k=1:Nz
spS4Au[index+1:index+10,1]=rand(10)
spS4bu[index4b+1,1]=1#-#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+10; index4b=index4b+1;
end
end
k=Nz;
for i=3:Nx-1 #rightest grids
for j=1:Ny #k=Nz
spS4Au[index+1:index+15,1]=rand(15)
spS4bu[index4b+1,1]=1
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #leftest grids
k=1;
for j=1:Ny #k=1
spS4Au[index+1:index+15,1]=rand(15);
spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=1
j=1;
for k=2:Nz-1
spS4Au[index+1:index+15,1]=rand(15);
spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
for i=3:Nx-1 #upper grids j=Ny
j=Ny;
for k=2:Nz-1
spS4Au[index+1:index+15,1]=rand(15);
spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
index=index+15; index4b=index4b+1;
end
end
Au =sparse(vec(spI4Au),vec(spJ4Au),vec(spS4Au),DOFu,DOFu);
bu = sparse(vec(spI4bu),vec(spJ4bu),vec(spS4bu),DOFu,1);
u=Au\Matrix(bu)
u= [zeros(Ny*Nz,1); u];
U=reshape(Matrix(u),Nx,Ny,Nz);
return U
end
``I mean, for large enough systems you will run out of memory even though you use a sparse storage format. How large is the matrix? And how sparse?
A = 1020100×1020100 SparseMatrixCSC{Float64,Int64} with 15199490 stored entries
b = 1020100×1 SparseMatrixCSC{Float64,Int64} with 1020100 stored entries:
But A is a block diagonal matrix not a diagonal matrix.
A and b can be stored but canot divide, the problem seems happen in LU part
ERROR: OutOfMemoryError()
Stacktrace:
[1] umfpack_numeric!(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:258
[2] #lu#1(::Bool, ::typeof(lu), ::SparseMatrixCSC{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:160
[3] lu at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:154 [inlined]
[4] \(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SparseArrays\src\linalg.jl:1374
[5] pdesolver(::Int64, ::Float64, ::Float64, ::Float64, ::Array{Float64,2}, ::Int64) at C:\Users\lyu\Downloads\test.jl:257
b can be stored in a Vector. Dont need sparse vector for this.
You should realize that there is an intrinsic limitation of sparse-direct methods like this for matrices that come from large 3d meshes/grids — they scale poorly for 3d problems, and tend to run out of memory before the problems become large. You could try a different sparse-direct solver like MUMPS (GitHub - JuliaSmoothOptimizers/MUMPS.jl: A Julia Interface to MUMPS or GitHub - JuliaInv/MUMPSjInv.jl: MUMPS interface for Julia), which might suffice in the short term, but in the long term you will hit the same scaling issues.
Here is a nice slide on the scaling (by Per Persson) that I use in one of my classes:
Note that for an n \times n \times n grid, N=n^3 so the storage scales as O(n^4). For more information, see e.g. Tim Davis’ book.
In consequence, for large 3d PDEs, the sparse-matrix problems must ultimately be solved by iterative methods (e.g. GMRES, CG, etc.) rather than sparse-direct methods. There are a few packages for iterative solvers in Julia, e.g. IterativeSolvers.jl, but to use iterative solvers effectively you have to know a bit more, and you often want to find a problem-specific preconditioner.
9 Likes