Loading 1d_diffusion.ipynb +19 −19 Original line number Diff line number Diff line %% Cell type:code id: tags: ``` julia include("src/fp.jl") include("src/flux_array.jl") ``` %% Output Main.FP Main.ModFluxArray %% Cell type:code id: tags: ``` julia using .FP using .ModFluxArray ``` %% Cell type:code id: tags: ``` julia N = 10 D = 1.; h = 1e-3; k = 1.; ``` %% Cell type:code id: tags: ``` julia n = DistributionArray{Float64}(0., 1., N); f = FluxArray{Float64}(N); ``` %% Cell type:code id: tags: ``` julia n.p[Integer(size(n)/2)] = 1; n.p ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_fr!(f); fill_lr!(f, n.p); ``` %% Cell type:code id: tags: ``` julia f.left ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.0022742765607139948 0.1936697622933632 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia f.right ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.997725723439286 0.8063302377066368 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia f.left .+ f.right ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_df!(f, D, h); ``` %% Cell type:code id: tags: ``` julia f.diffusive ``` %% Output 11-element Array{Float64,1}: 0.0 0.0 0.0 0.0 2.2742765607139948 -997.725723439286 193.6697622933632 -806.3302377066368 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_sf!(f, D, h, k); ``` %% Cell type:code id: tags: ``` julia f.stochastic ``` %% Output 11-element Array{Float64,1}: -0.8640646960498399 0.0 1.662956306086986 0.0 -0.0 0.3672085709169771 -2.924360797089483 0.0 0.0 -0.0 4.240011808025856 2.710972258785979 -0.0 0.0 0.3838915655069409 -0.0 -0.0 -0.3416384648516665 %% Cell type:code id: tags: ``` julia apply_bc!(f); ``` %% Cell type:code id: tags: ``` julia f.stochastic ``` %% Output 11-element Array{Float64,1}: 0.0 0.0 0.0 -0.0 0.3672085709169771 -2.924360797089483 0.0 0.0 -0.0 4.240011808025856 2.710972258785979 -0.0 0.0 -0.0 -0.0 0.0 %% Cell type:code id: tags: ``` julia ``` src/fp.jl→src/flux_array.jl +2 −52 Original line number Diff line number Diff line module FP module StochFluxArray export DistributionArray export size Loading Loading @@ -172,7 +172,7 @@ function fill_ft!(flx::FluxArray{T}, h::T, k::T) where T for i=2:size_cc(flx) # Total particle number P::T = flx.right[i-1] + flx.left[i] P::T = flx.right[i-1] + flx.left[i]; # Number in right box (i.e. number of particles going left) x0::T = flx.left[i]; # Valid range of fluxes Loading @@ -182,55 +182,6 @@ function fill_ft!(flx::FluxArray{T}, h::T, k::T) where T end #function getFlux(n, h, N, k, D, FP, time) # #get the diffusive flux for FV solver # # #create flux vectors # Ftot = zeros(1,N+1); #total flux # Fd = zeros(1,N+1); #diffusive flux # Fs = zeros(1,N+1); #stochastic flux # # #get the fraction of particles going in each direction # fr = rand(1,N); #fraction going right, 1-fr is fraction going left # left = n.*(1-fr);right = n.*fr; #define left and right particle numbers # # #compute the diffusive flux # for i = 2:N # Fd(i) = D/h * (left(i) - right(i-1)); #random fractions # #Fd(i) = D/h * (n(i) - n(i-1)); #usual fv # end # # #compute the stochastic flux - assume gaussian # R = randn(1,N+1); #standard normals # for i = 2:N # avg = left(i) + right(i-1); #average over cell, ignore factor of 2 # Fs(i) = R(i)*sqrt(D*avg/(k*sqrt(h))); #set correct variance, 2 cancels # end # # #check if each box has a physical flux # nu = k/h; #ratio of time to space step to propogate soln # Ftot = Fd + Fs; #set total flux # Q = -Ftot*nu; #take time step to get real fluxes # for i=1:N-1 # P = right(i)+left(i+1); #total particle number # x0 = left(i+1); #number in right box # lower = -x0; upper = P-x0; #valid range of fluxes # if (Q(i+1)>upper || Q(i+1) < lower) # #if particle num is low enough, sample uniform - approx FP # if (P < 1e-3) # Ftot(i+1) = -1/nu*(P*rand()-x0); # else # #do a fokker planck solve - sample flux # #divide by nu to remain consistent. # flux = FPsample(P,x0,D,h,k); # Ftot(i+1) = -1/nu*flux; # FP{i}(time) = 1; # end # end # end # # return Ftot, FP #end end No newline at end of file Loading
1d_diffusion.ipynb +19 −19 Original line number Diff line number Diff line %% Cell type:code id: tags: ``` julia include("src/fp.jl") include("src/flux_array.jl") ``` %% Output Main.FP Main.ModFluxArray %% Cell type:code id: tags: ``` julia using .FP using .ModFluxArray ``` %% Cell type:code id: tags: ``` julia N = 10 D = 1.; h = 1e-3; k = 1.; ``` %% Cell type:code id: tags: ``` julia n = DistributionArray{Float64}(0., 1., N); f = FluxArray{Float64}(N); ``` %% Cell type:code id: tags: ``` julia n.p[Integer(size(n)/2)] = 1; n.p ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_fr!(f); fill_lr!(f, n.p); ``` %% Cell type:code id: tags: ``` julia f.left ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.0022742765607139948 0.1936697622933632 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia f.right ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.997725723439286 0.8063302377066368 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia f.left .+ f.right ``` %% Output 10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_df!(f, D, h); ``` %% Cell type:code id: tags: ``` julia f.diffusive ``` %% Output 11-element Array{Float64,1}: 0.0 0.0 0.0 0.0 2.2742765607139948 -997.725723439286 193.6697622933632 -806.3302377066368 0.0 0.0 0.0 0.0 0.0 %% Cell type:code id: tags: ``` julia fill_sf!(f, D, h, k); ``` %% Cell type:code id: tags: ``` julia f.stochastic ``` %% Output 11-element Array{Float64,1}: -0.8640646960498399 0.0 1.662956306086986 0.0 -0.0 0.3672085709169771 -2.924360797089483 0.0 0.0 -0.0 4.240011808025856 2.710972258785979 -0.0 0.0 0.3838915655069409 -0.0 -0.0 -0.3416384648516665 %% Cell type:code id: tags: ``` julia apply_bc!(f); ``` %% Cell type:code id: tags: ``` julia f.stochastic ``` %% Output 11-element Array{Float64,1}: 0.0 0.0 0.0 -0.0 0.3672085709169771 -2.924360797089483 0.0 0.0 -0.0 4.240011808025856 2.710972258785979 -0.0 0.0 -0.0 -0.0 0.0 %% Cell type:code id: tags: ``` julia ```
src/fp.jl→src/flux_array.jl +2 −52 Original line number Diff line number Diff line module FP module StochFluxArray export DistributionArray export size Loading Loading @@ -172,7 +172,7 @@ function fill_ft!(flx::FluxArray{T}, h::T, k::T) where T for i=2:size_cc(flx) # Total particle number P::T = flx.right[i-1] + flx.left[i] P::T = flx.right[i-1] + flx.left[i]; # Number in right box (i.e. number of particles going left) x0::T = flx.left[i]; # Valid range of fluxes Loading @@ -182,55 +182,6 @@ function fill_ft!(flx::FluxArray{T}, h::T, k::T) where T end #function getFlux(n, h, N, k, D, FP, time) # #get the diffusive flux for FV solver # # #create flux vectors # Ftot = zeros(1,N+1); #total flux # Fd = zeros(1,N+1); #diffusive flux # Fs = zeros(1,N+1); #stochastic flux # # #get the fraction of particles going in each direction # fr = rand(1,N); #fraction going right, 1-fr is fraction going left # left = n.*(1-fr);right = n.*fr; #define left and right particle numbers # # #compute the diffusive flux # for i = 2:N # Fd(i) = D/h * (left(i) - right(i-1)); #random fractions # #Fd(i) = D/h * (n(i) - n(i-1)); #usual fv # end # # #compute the stochastic flux - assume gaussian # R = randn(1,N+1); #standard normals # for i = 2:N # avg = left(i) + right(i-1); #average over cell, ignore factor of 2 # Fs(i) = R(i)*sqrt(D*avg/(k*sqrt(h))); #set correct variance, 2 cancels # end # # #check if each box has a physical flux # nu = k/h; #ratio of time to space step to propogate soln # Ftot = Fd + Fs; #set total flux # Q = -Ftot*nu; #take time step to get real fluxes # for i=1:N-1 # P = right(i)+left(i+1); #total particle number # x0 = left(i+1); #number in right box # lower = -x0; upper = P-x0; #valid range of fluxes # if (Q(i+1)>upper || Q(i+1) < lower) # #if particle num is low enough, sample uniform - approx FP # if (P < 1e-3) # Ftot(i+1) = -1/nu*(P*rand()-x0); # else # #do a fokker planck solve - sample flux # #divide by nu to remain consistent. # flux = FPsample(P,x0,D,h,k); # Ftot(i+1) = -1/nu*flux; # FP{i}(time) = 1; # end # end # end # # return Ftot, FP #end end No newline at end of file