Newest questions tagged advection-diffusion - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2019-06-27T10:09:07Z https://scicomp.stackexchange.com/feeds/tag?tagnames=advection-diffusion&sort=newest http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://scicomp.stackexchange.com/q/31375 0 Error for the finite differences scheme -- Advection equation Smilia https://scicomp.stackexchange.com/users/29701 2019-04-05T06:36:25Z 2019-04-06T18:56:32Z <p>Consider the advection equation (1D in space) <span class="math-container">$$\frac{\partial u}{\partial t} + V\, \frac{\partial u}{\partial x}=0$$</span> and we solve it numerically on <span class="math-container">$[0,1]\times [0,1]\ni (t,x)$</span> using a forward scheme in time and a backward decentered in space. This scheme is of order <span class="math-container">$1$</span> in space and of order <span class="math-container">$1$</span> in time and is stable iff <span class="math-container">$V\frac{dt}{dx}\le 1$</span>.</p> <p>Let <span class="math-container">$N$</span> be the number of points in the subdivision of <span class="math-container">$[0,1]$</span> and <span class="math-container">$dt$</span> be the time step. Numerically, to observe the order, if I fix <span class="math-container">$N$</span> and <span class="math-container">$dt$</span> is varying the error between the analyitc solution and the numerical solution is represented in the log-plot figure:</p> <p>The red curve is the linear regression between the points, but useless here because we don't observe a straight line so the order can be checked. More importantly when <span class="math-container">$dt$</span> is decreasing, the error is increasing ... <strong>What is the reason for that ?</strong></p> <ul> <li>problem in the implementation ? maybe</li> <li>the fact is when we say of order <span class="math-container">$1$</span> in space and time it means we approximate the equation by an error which can be written: <span class="math-container">$$C_1\, dt + C_2\, dx$$</span> where <span class="math-container">$C_1$</span> and <span class="math-container">$C_2$</span> are constant. But <span class="math-container">$C_1$</span> may depend on <span class="math-container">$dx$</span> and <span class="math-container">$C_2$</span> may depend on <span class="math-container">$dt$</span>. So when <span class="math-container">$dt$</span> is decreasing, <span class="math-container">$C_2$</span> is not constant and the error can increase, which could explain the plot.</li> </ul> <p><a href="https://i.stack.imgur.com/GdtGP.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/GdtGP.png" alt="(blue) Error when dt varies and N is fixed"></a></p> <p>Note that when the ratio <span class="math-container">$V\frac{dt}{dx}$</span> is constant equal to <span class="math-container">$1$</span> and <span class="math-container">$dt,dx \to 0$</span>, the error is constant and very small (no diffusion)</p> https://scicomp.stackexchange.com/q/31317 0 Advection-Diffusion by using Lattice Boltzmann Method, Is it practical for engineering applications? Alone Programmer https://scicomp.stackexchange.com/users/28872 2019-03-26T19:44:07Z 2019-03-26T19:53:06Z <p>I want to use lattice Boltzmann method to solve advection-diffusion in three-dimensional space. In fact, my problem is related to drug release in human blood vessels and as a results, I'm interested to use lattice Boltzmann to solve this equation:</p> <p><span class="math-container">$$\frac{\partial \phi}{\partial t} + \vec{u} \cdot \nabla \phi = \nabla \cdot (D \nabla \phi)$$</span></p> <p>Where <span class="math-container">$\phi$</span> is concentration of drug, <span class="math-container">$\vec{u}$</span> is blood velocity derived from other lattice Boltzmann flow solver (for a moment let's forgot about the coupling of advection-diffusion and Navier-Stokes and assume that <span class="math-container">$\vec{u}$</span> is already a known vector field), and <span class="math-container">$D$</span> is diffusion coefficient of drug in blood.</p> <p>Due to that <span class="math-container">$\phi$</span> is just a scalar variable, I use lattice Boltzmann equation with single-relaxation time collision operator like this:</p> <p><span class="math-container">$$f_{i}(\vec{r}+\vec{c}_{i}\Delta t,t+\Delta t) = f_{i}(\vec{r},t) - \frac{f_{i}-f^{eq}_{i}}{\tau}$$</span></p> <p>Where <span class="math-container">$f_{i}$</span> is distribution function and <span class="math-container">$f^{eq}_{i}$</span> is the equilibrium distribution function derived from Maxwell-Boltzmann distribution and expanded up to first order like this:</p> <p><span class="math-container">$$f^{eq}_{i} = \omega_{i}\phi(1+\frac{\vec{u} \cdot \vec{c}_{i}}{c_{s}^{2}})$$</span></p> <p>Where <span class="math-container">$\omega_{i}$</span> is the ith weight in the lattice and <span class="math-container">$c_{s}$</span> is the speed of sound of the lattice.</p> <p>For my practical applications, <span class="math-container">$D$</span> (diffusion coefficient of drug in blood) is in the order of ~<span class="math-container">$10^{-11}$</span> <span class="math-container">$\frac{\mathrm{m}^{2}}{\mathrm{s}}$</span>. The blood velocity is in the order of 0.1 <span class="math-container">$\frac{\mathrm{m}}{\mathrm{s}}$</span>. The characteristic length of the blood vessel (L), which is its diameter, is 3 mm. As a result, the Peclet number could be calculated as:</p> <p><span class="math-container">$$Pe = \frac{u L}{D} = \frac{0.1 \times 3 \times 10^{-3}}{10^{-11}} = 3\times 10^{7}$$</span></p> <p>Typically in lattice Boltzmann simulations, people use <span class="math-container">$Mach &lt; 0.02$</span> to ensure stability but due to my extremely complicated case, I want to push the borders and say Mach number could have some freedom to be as high as 0.1. So, <span class="math-container">$Mach = 0.1$</span>.</p> <p>On the other hand, we know that <span class="math-container">$\tau$</span> the relxation time of lattice Boltzmann scheme is related to diffusion coefficient as:</p> <p><span class="math-container">$$\tau = \frac{\Delta t}{2} + \frac{1}{c_{s}} \frac{D}{c_{s}} = \frac{\Delta t}{2} + \frac{L}{c_{s}} \frac{Mach}{Pe}$$</span></p> <p>Due to calculated Pe number, assumed Mach number, and characteristic length of the blood vessel and the fact that at least for stability of this lattice Boltzmann scheme, we need to have <span class="math-container">$\tau &gt; 0.5 \Delta t$</span>, we would take <span class="math-container">$\tau$</span> as <span class="math-container">$\tau = 0.501 \Delta t$</span> (which is really small but still acceptable due to linear equilibrium distribution function used in this study). Finally, we have:</p> <p><span class="math-container">$$\tau = 0.501 \Delta t = 0.5 \Delta t + \frac{L}{c_{s}} \frac{Mach}{Pe}$$</span></p> <p><span class="math-container">$$0.001 = \frac{3\times 10^{-3}\times 0.1}{\frac{1}{\sqrt{3}}\Delta x \times 3 \times 10^{7}}$$</span></p> <p><span class="math-container">$$\Delta x = 1.732 \times 10^{-8}\mathrm{m}$$</span></p> <p>Where <span class="math-container">$\Delta x$</span> is the mesh size.</p> <p>Now, if we assume that the length of the cylinder which represent the blood vessel is 30 mm, we could estimate the total number of meshes that are needed to fill this cylinder as:</p> <p><span class="math-container">$$N = \frac{V}{\Delta x^{3}} = \frac{\frac{\pi}{4} L^{2} H}{\Delta x^{3}} = 4.081 \times 10^{16}$$</span></p> <p>It means, I need <span class="math-container">$4.081 \times 10^{16}$</span> meshes to fill this cylinder and make sure my simulation by using this scheme will remain stable (still there is a doubt that my simulation will remain stable or not with even this unrealistic number of meshes). So, finally my question: <strong>Am I missing something here or something is wrong with my calculations and theoretical scheme, or really lattice Boltzmann method is unsuitable for this purpose (simulating drug diffusion in blood vessels in the presence of blood perfusion)?</strong></p> <p>I appreciate any insight or suggestion or useful direction.</p> https://scicomp.stackexchange.com/q/31268 0 How to make a less diffusive code to solve 2D advection equation? Peanutlex https://scicomp.stackexchange.com/users/26454 2019-03-20T10:51:38Z 2019-04-20T17:02:55Z <p>I would like to solve the following differential equation numerically in 2D, <span class="math-container">$$\frac{\partial z^-}{\partial t}+(\vec{B}\cdot\vec{\nabla})z^-=0,$$</span> see <a href="https://en.wikipedia.org/wiki/Magnetohydrodynamic_turbulence" rel="nofollow noreferrer">Wikipedia</a> if you are curious about what the symbols mean, where <span class="math-container">$$\vec{B}=(B_x,B_y,0).$$</span> I have attached some code which can solve this to first order accuracy in python by following the corner transport upstream (CTU) algorithm outlined <a href="https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-950-atmospheric-and-oceanic-modeling-spring-2004/lecture-notes/lec10.pdf" rel="nofollow noreferrer">here</a>. The problem is the code is too diffusive for my purposes. Does anyone know of an algorithm I could follow to make the code a higher order of accuracy while still remaining stable? Also, it would nice if the algorithm was an upstream algorithm i.e. each cell update only required information from cells in the upstream direction as it is easier to impose open boundary conditions if this is the case.</p> <p>Here is the python code:</p> <pre><code>import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import os nx = 64 x_min = -2 x_max = 2 dx = (x_max - x_min) / (nx - 1) x = np.linspace(x_min, x_max, nx) ny = 64 y_min = -2 y_max = 2 dy = (y_max - y_min) / (ny - 1) y = np.linspace(y_min, y_max, ny) bx = np.zeros((nx,ny)) + 1 by = np.zeros((nx,ny)) + 1 t_max = 1 dt = dx * dy / np.max(abs(bx) * dy + abs(by) * dx) # CFL condition nt = int(t_max / dt) t = np.linspace(0, t_max, nt) zm = np.zeros((nx, ny, nt)) for i in range(nx): for j in range(ny): r = np.sqrt(x[i] ** 2 + y[j] ** 2) if r &lt; 1: zm[i,j,0] = 2 * np.cos(np.pi * r / 2) ** 2 print(nt) zm1 = np.zeros((nx, ny)) for n in range(nt - 1): if n % 100 == 0: print('n =', n) cx = bx[1:,1:] * dt / dx cy = by[1:,1:] * dt / dy zm1[1:,1:] = (1 - cx) * zm[1:,1:,n] + cx * zm[0:-1,1:,n] zm[1:,1:,n+1] = (1 - cy) * zm1[1:,1:] + cy * zm1[1:,0:-1] zm[0,:,n+1] = 0 zm[:,0,n+1] = 0 print('Computation finished, now generating figures...') os.makedirs('Figures/2d/zm', exist_ok = True) X = np.zeros((nx,ny)) Y = np.zeros((nx,ny)) for j in range(ny): X[:,j] = x for i in range(nx): Y[i,:] = y i = -1 for n in range(nt - 1): # if n % 10 == 0: i = i + 1 fig = plt.figure() ax = fig.gca(projection = '3d') ax.plot_surface(X, Y, zm[:,:,n]) plt.title('t = ' + str(t[n])) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('zp') ax.set_zlim(-1, 1) plt.savefig('Figures/2d/zm/' + "{0:0=4d}".format(i) + '.png') plt.close(fig) </code></pre> <p><strong>My attempt at Beam-Warming (see comments):</strong></p> <pre><code>import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import os nx = 68 x_min = 0 x_max = 2 dx = (x_max - x_min) / (nx - 1) x = np.linspace(x_min - dx, x_max + dx, nx) # Add ghost cells ny = 68 y_min = -2 y_max = 0 dy = (y_max - y_min) / (ny - 1) y = np.linspace(y_min - dy, y_max + dy, ny) bx = np.zeros((nx,ny)) + 1 by = np.zeros((nx,ny)) + 1 t_max = 10 dt = dx * dy / np.max(abs(bx) * dy + abs(by) * dx) # CFL condition nt = int(t_max / dt) t = np.linspace(0, t_max, nt) zm = np.zeros((nx, ny, nt)) zp = np.zeros((nx, ny, nt)) for i in range(nx): for j in range(ny): r = np.sqrt((x[i] - 1) ** 2 + (y[j] + 1) ** 2) if r &lt; 1: zm[i,j,0] = 2 * np.cos(np.pi * r / 2) ** 2 print(nt) for n in range(nt - 1): if n % 100 == 0: print('n =', n) zm[2:,2:,n+1] = zm[2:,2:,n] \ - 0.5 * dt * bx[2:,2:] / dx * \ (3 * zm[2:,2:,n] - 4 * zm[1:-1,2:,n] + zm[0:-2,2:,n]) \ + 0.5 * (dt * bx[2:,2:] / dx) ** 2 * \ (zm[2:,2:,n] - 2 * zm[1:-1,2:,n] + zm[0:-2,2:,n]) \ - 0.5 * dt * by[2:,2:] / dy * \ (3 * zm[2:,2:,n] - 4 * zm[2:,1:-1,n] + zm[2:,0:-2,n]) \ + 0.5 * (dt * by[2:,2:] / dy) ** 2 * \ (zm[2:,2:,n] - 2 * zm[2:,1:-1,n] + zm[2:,0:-2,n]) zm[0,:,n+1] = 0 zm[:,0,n+1] = 0 print('Computation finished, now generating figures...') os.makedirs('Figures/2d/zm', exist_ok = True) X = np.zeros((nx,ny)) Y = np.zeros((nx,ny)) for j in range(ny): X[:,j] = x for i in range(nx): Y[i,:] = y i = -1 for n in range(nt - 1): if n % 100 == 0: i = i + 1 fig = plt.figure() ax = fig.gca(projection = '3d') ax.plot_surface(X, Y, zm[:,:,n]) plt.title('t = ' + str(t[n])) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('zp') ax.set_zlim(-1, 1) plt.savefig('Figures/2d/zm/' + "{0:0=4d}".format(i) + '.png') plt.close(fig) </code></pre> <p>Pictures of the instability: <a href="https://i.stack.imgur.com/T900m.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/T900m.png" alt=""></a> <a href="https://i.stack.imgur.com/sVc0A.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/sVc0A.png" alt="enter image description here"></a></p> https://scicomp.stackexchange.com/q/30818 2 Analytical Solution of Transport Equation Natasha https://scicomp.stackexchange.com/users/29087 2019-01-01T13:04:09Z 2019-01-02T06:20:59Z <p>I'm looking at the analytical solution of the convection-diffusion equation</p> <p><span class="math-container">$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x}$$</span> with initial condition <span class="math-container">$$c(x,0) = 4000$$</span> and with Dirichlet boundary condition</p> <p><span class="math-container">$$C(x=0,t&gt;0) = 4100$$</span></p> <p>Neumann boundary condition <span class="math-container">$$\frac{\partial C}{\partial x}=0\text{ at } x=L ; t&gt;0.$$</span></p> <p><a href="https://i.stack.imgur.com/8xQzi.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8xQzi.png" alt="enter image description here"></a></p> <p>However, when the above analytical solution is coded I obtain negative values of C which is unrealistic.</p> <pre><code>function AnalyticalSoln2() format long R = 1; D = 900; v = 200; L = 60; co = 4100; ci = 4000; X = linspace(0,60,10); t = 0:0.001:2; sol=[]; for pos = 1:length(X) x = X(pos); A1 = 0.5*erfc((x.*R-t.*v)./(2*(t.*D*R).^0.5)); A2 = 0.5*exp(x.*v/D).*erfc((x.*R+t.*v)./(2*(D*t.*R).^0.5)); A31 = 0.5*(2+v*(2*L-x)/D + (t.*v^2)/(D*R))*exp(v*L/D); A32 = erfc((R*(2*L-x)+t.*v)./(2*(D*t.*R).^0.5)); A41 = -((t.*v^2)./(pi*D*R)).^0.5; A42 = exp((v*L/D)-(R./(4*t.*D)).*(2*L -x + t.*v/R).^2); A = A1 + A2 + A31.*A32 + A41.*A42; if t==0 C = ci + (co - ci)*A' else C = ci + (co - ci)*A' - co*A'; end sol = horzcat(sol,C); end sol(1:100:end,:) end </code></pre> <p>Whereas, the numerical solution obtained from the pdepe solver is non-negative.</p> <p>Here's the solution obtained using pdepe</p> <pre><code>function DiffusionConvectionMATLAB format short global D m = 0; x = linspace(0,60,10); t = 0.1:0.1:2; D = 900; sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t) function [g,f,s] = pdefun(x,t,c,DcDx) v = 200; g = 1; f = D*DcDx; s = -v*DcDx; end function c0 = icfun(x) c0 = 4000; end function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t) pl = cl-4100; ql = 0; pr = 0; qr = 1; end end </code></pre> <p>Could someone suggest if there is any issue with the way I'm computing the solution form the analytical expression?</p> https://scicomp.stackexchange.com/q/30817 0 Simulating Brownian motion in 3-D for first hitting time? Userhanu https://scicomp.stackexchange.com/users/23559 2019-01-01T10:41:30Z 2019-01-02T10:40:04Z <p>I want to simulate Brownian motion in 3-D for the following conditions: <span class="math-container">$$p(x=0,y=0,z=0,t=0)=1$$</span> <span class="math-container">$$p(x,y,z=c,t)=0$$</span> where <span class="math-container">$p$</span> is the probability of finding molecules in the 3-D environment. I want to find the number of particles getting absorbed on a certain region on the absorbing wall in 5 secs. As can be inferred the absorbing boundary is at <span class="math-container">$z=c$</span>, and I want the probability of molecules hitting a circular area centered at <span class="math-container">$(0,0,2.5\times10^{-4})$</span> with a radius <span class="math-container">$5\times10^{-6}$</span>m. The fluid environment has drift,i.e. <span class="math-container">$u,v,w$</span> in three directions. The particles in <span class="math-container">$x,y,z$</span> directions are moving as <span class="math-container">\begin{equation} \text{d}X_j(t)=u\text{d}t+\sqrt{2D}\text{d}B_j(t),\\ \text{d}Y_j(t)=v\text{d}t+\sqrt{2D}\text{d}B_j(t),\\ \text{d}Z_j(t)=w\text{d}t+\sqrt{2D}\text{d}B_j(t), \end{equation}</span> where <span class="math-container">$\text{d}B$</span> is the derivative of Brownian process which I am simulating through normal random varaible.</p> <p>I am using following MATLAB code to simulate particle diffusion:</p> <pre><code>close all; clear all; clc; TotalNumberReceivedUp=0; Ts=5; TotalSimulations1=100000; flg=0; for particlesNum=1:TotalSimulations1 flg=flg+1 D=4*10^-9; c=250*10^-6;u=5*10^-6;v=2*10^-6;w=3*10^-6; e=5*10^-6; X=zeros;Y=zeros;Z=zeros; X(1)=0;Y(1)=0;Z(1)=0; j=1; for i=5*10^-4:5*10^-4:Ts t2=i; t1=t2-5*10^-4; j=j+1; r = normrnd(0,1); q= normrnd(0,1); p= normrnd(0,1); X(j)=X(j-1)+(sqrt(2*D))*sqrt(t2-t1)*p+u*(t2-t1); Y(j)=Y(j-1)+(sqrt(2*D))*sqrt(t2-t1)*q+v*(t2-t1); Z(j)=Z(j-1)+(sqrt(2*D))*sqrt(t2-t1)*r+w*(t2-t1); if ((Z(j)&gt;=c) &amp;&amp; ... (sqrt((X(j))^2+(Y(j))^2)&lt;=e)) TotalNumberReceivedUp=TotalNumberReceivedUp+1; break end end end ProbabilityUp=TotalNumberReceivedUp/TotalSimulations1 </code></pre> <p>But my results are not coming good, the probable particles hitting are more than those I am getting through analysis. I think it may be because of the excess over the boundary problem which has been discussed in "A comparison of four methods for simulating the diffusion process" for 1-D movement with two absorbing barriers. I have tried to use the same factor in 3-D but it doesn't work well. </p> https://scicomp.stackexchange.com/q/30744 0 Analytical solution of 1D advection -diffusion equation Natasha https://scicomp.stackexchange.com/users/29087 2018-12-18T06:23:52Z 2018-12-19T18:50:59Z <p>I am looking for the analytical solution of 1-dimensional advection-diffusion equation with Neumann boundary condition at both the inlet and outlet of a cylinder through which the fluid flow occurs. <span class="math-container">$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x}$$</span> with initial condition <span class="math-container">$$c(x,0) = C_i$$</span> and with Neumann boundary condition <span class="math-container">$$\frac{\partial C}{\partial x}=0\text{ at }t&gt;0.$$</span></p> <p>Could someone suggest a reference?</p> <p>I had a chance to look at the answer posted <a href="https://scicomp.stackexchange.com/questions/7873/does-the-time-dependent-advection-diffusion-equation-have-an-analytical-solution">here</a>. Out of the solutions listed, I couldn't find the analytical solution for the transport equation with Neumann boundary condition at both the ends.</p> https://scicomp.stackexchange.com/q/30605 0 Impose Neumann Boundary Condition in advection-diffusion equation 1D tnt235711 https://scicomp.stackexchange.com/users/28594 2018-11-28T09:02:32Z 2018-11-30T04:14:12Z <p>when solving the advection equation in 1D that is: </p> <p><span class="math-container">$$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$$</span> with <span class="math-container">$u'(t,0) = 0$</span> and <span class="math-container">$u(t,L) = 0$</span> , <span class="math-container">$u(0,x) = u_{0}$</span></p> <p>one numerical scheme is the FTCS (Forward time-centered space), but this numerical scheme is unstable. </p> <p><span class="math-container">$$\frac{u_{j}^{n+1}-u_{j}^{n}}{ h_{t}} = c \frac{u_{j+1}^{n}-u_{j-1}^{n}}{ 2h_{x}}$$</span></p> <p>But when solving </p> <p><span class="math-container">$$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \alpha \frac{\partial^2 u}{\partial x^2}$$</span> the advection-diffusion equation in 1D with <span class="math-container">$u'(t,0) = 0$</span> and <span class="math-container">$u(t,L) = 0$</span> , <span class="math-container">$u(0,x) = u_{0}$</span></p> <p>Since the advection-diffusion equation is a second order equation I'd like to use a second order approximation. </p> <p>if we define <span class="math-container">$u_{k}^{n} := u(t_{n},x_{k})$</span>; <span class="math-container">$\ \ x_{k} = kh$</span> and <span class="math-container">$\ \ k = 0,1,2,...,N$</span>. <span class="math-container">$h$</span> is known as the mesh size or step size. </p> <p>For the second derivative:</p> <p><span class="math-container">$$\frac{\partial^2u_{k}^{n}}{\partial x^2} \approx \frac{u_{k+1}^{n}-2u_{k}^{n}+u_{k-1}^{n}}{h^2} = \frac{u_{k-1}^{n}-2u_{k}^{n}+u_{k+1}^{n}}{h^2}$$</span> for <span class="math-container">$k=0,1,...,N-1$</span> </p> <p>Since <span class="math-container">$u'(t,0) = 0$</span> and <span class="math-container">$u(t_{n},L) = u(t_{n},x_{N}) = u_{N}^{n} = 0$</span> we get the following matrix representation of the second derivative operator</p> <p><span class="math-container">\begin{equation} \frac{\partial^2}{\partial x^2} \approx L_{2} = \frac{1}{h^2}\left(\begin{matrix} -2 &amp; 1 &amp; &amp; 0\\ 1 &amp; \ddots &amp; \ddots &amp; \\ &amp; \ddots &amp; \ddots &amp; 1 \\ 0 &amp; &amp; 1 &amp; -2 \end{matrix} \right) \end{equation}</span></p> <p>for <span class="math-container">$k=0$</span> , we get </p> <p><span class="math-container">$$\frac{\partial u_{0}^{n}}{\partial x} = \frac{u_{0+1}^{n}-u_{0-1}^{n}}{ 2h} = 0$$</span> this implies that <span class="math-container">$u_{1}^{n}=u_{-1}^{n}$</span> and </p> <p><span class="math-container">$$\frac{\partial^2u_{0}^{n}}{\partial x^2} \approx \frac{u_{0+1}^{n}-2u_{0}^{n}+u_{0-1}^{n}}{h^2} = \frac{u_{0-1}^{n}-2u_{0}^{n}+u_{0+1}^{n}}{h^2} = \frac{-2u_{0}^{n}+2u_{1}^{n}}{h^2}$$</span></p> <p>thus we have to modify the entry <span class="math-container">$1,2$</span> of <span class="math-container">$L_{2}$</span></p> <p><span class="math-container">\begin{equation} L_{2} = \frac{1}{h^2}\left(\begin{matrix} -2 &amp; 2 &amp; &amp; 0\\ 1 &amp; \ddots &amp; \ddots &amp; \\ &amp; \ddots &amp; \ddots &amp; 1 \\ 0 &amp; &amp; 1 &amp; -2 \end{matrix} \right) \end{equation}</span></p> <p>What I have done, is <span class="math-container">$\mathbf{impose \ the \ Neumann \ boundary \ condition}$</span> in <span class="math-container">$L_{2}$</span> . </p> <p>I want to approximate the first derivative using central difference(Second order approximation):</p> <p><span class="math-container">$$\frac{\partial u_{k}^{n}}{\partial x} = \frac{u_{k+1}^{n}-u_{k-1}^{n}}{ 2h}$$</span></p> <p>The matrix representation is: <span class="math-container">\begin{equation} \frac{\partial}{\partial x} \approx L_{1} = \frac{1}{2h}\left(\begin{matrix} 0 &amp; 1 &amp; &amp; 0\\ -1 &amp; \ddots &amp; \ddots &amp; \\ &amp; \ddots &amp; \ddots &amp; 1 \\ 0 &amp; &amp; -1 &amp; 0 \end{matrix} \right) \end{equation}</span></p> <p>for <span class="math-container">$k=0$</span> , we get </p> <p><span class="math-container">$$\frac{\partial u_{0}^{n}}{\partial x} = \frac{u_{0+1}^{n}-u_{0-1}^{n}}{ 2h} = 0$$</span> this implies that <span class="math-container">$u_{1}^{n}=u_{-1}^{n}$</span> </p> <p>But I'm stuck when I try to <span class="math-container">$\mathbf{impose \ the \ neumann \ boundary \ condition \ in}$</span> <span class="math-container">$L_{1}$</span>. I don't know how to do that.</p> <p>If we solve that problem, we can solve the differential equation</p> <p><span class="math-container">$$\frac{ \partial u }{\partial t} = \Big(-cL_{1} +\alpha L_{2} \Big)u$$</span> integrating in time( C-N, Back-Euler, RK4 )</p> <p><span class="math-container">$\mathbf{please \ help!\ \ How \ do \ you \ impose \ the \ Neumann \ Boundary \ Condition \ in \ L_{1}?}$</span></p> https://scicomp.stackexchange.com/q/30152 1 Inflow and outflow boundary conditions for advection-diffusion equation Mehrdad Yousefi https://scicomp.stackexchange.com/users/28407 2018-09-04T20:47:54Z 2018-09-04T20:47:54Z <p>I'm trying to solve this advection-diffusion equation (ADE):</p> <p>$$\frac{\partial \phi}{\partial t} + \nabla \cdot (-D \nabla \phi + \mathbf{u} \phi) = 0$$</p> <p>In fact, this ADE framework is coupled to a Navier-Stokes (NS) solver and I get velocity field ($\mathbf{u}$) from NS solver. My question is about ADE boundary conditions on inflow and outflow planes. I postulate that because the flow is directed from inflow to outflow obviously, I put zero concentration ($\phi = 0$) boundary condition on inlet plane because my source of species is located far from inflow boundary. Also, I assumed that the dominant transport mechanism at the outlet plane is convection and as a result of that normal diffusive flux will be vanished at the outlet plane ($-D \frac{\partial \phi}{\partial \vec{n}} = 0$).</p> <p>I wanted to know how much these assumptions are true and could be justified physically? I should say that the outflow boundary condition is verified by comparing the simulation results with analytical solution of the Graetz problem. Any reference, suggestion, or idea is appreciated.</p> https://scicomp.stackexchange.com/q/30099 1 How write a integration loop in fortran, leapfrog scheme to solvind PDE (advection)? PCat27 https://scicomp.stackexchange.com/users/28581 2018-08-28T00:57:59Z 2018-09-28T20:01:06Z <p>I want to resolve numerically this equation using of difference finite method with Leapfrog Scheme $$\frac{\partial{u}}{\partial t}+ v \frac{\partial{u}}{\partial x}= 0$$</p> <p>I'm trying to write a code that resolves that PDE wit $v=1$ in leapfrog scheme (<a href="http://www4.ncsu.edu/~zhilin/TEACHING/MA584/Chapter5_fd_fem.pdf" rel="nofollow noreferrer">http://www4.ncsu.edu/~zhilin/TEACHING/MA584/Chapter5_fd_fem.pdf</a> page 7), this scheme uses two level of time, but, I don't know how to write it in my code</p> <p>$$u_j^{n+1}=u_j^{n-1}-\frac{\Delta t}{\Delta x}\left[u_{j+1}^{n}-u_{j-1}^{n}\right]+\mathcal O(\Delta t^3,\Delta x^2)$$ So I guess this is how could be.</p> <pre><code>do l=1,Nt t = t + dt u_n_p_o=u_n_o u_n_p=u_n </code></pre> <p>I calculated $u^{n-1}$ to using FTCS(centering difference) Scheme</p> <pre><code>do i=1,Nx-1 u_n_o(i)= u_n_p_o(i)- (dt/dx)*(u_n_p_o(i+1)- u_n_p_o(i-1)) end do </code></pre> <p>Actually <code>u_n_o</code> is correct ?</p> <pre><code>u_n_p_o(i)=u_n_o(i-1) !2do reciclado de variable do i=1,Nx-1 u_n(i)= u_n_p_o(i)- (dt/dx)*(u_n_p(i+1)-u_n_p(i-1)) end do u_e = amp * exp( - ( x - x0 - t)**2 / sigma**2 ) + 1e-20 call save1Ddata(Nx,t,x,u_e,'u_e',1) call save1Ddata(Nx,t,x,u_n,'u_n',1) call save1Ddata(Nx,t,x,u_n,'u_n_o',1) end do </code></pre> <p>What are the features of Leapfrog? I know that in this case, it's no matter the Courant factor (because it is stable); however, the solution numerically dissipates in its amplitude. So how can you test your numeric solution? </p> https://scicomp.stackexchange.com/q/29931 5 Finite Differencing schemes for Convection-Diffusion equation Scotty1- https://scicomp.stackexchange.com/users/27868 2018-07-26T18:50:15Z 2018-08-11T04:36:43Z <p>I'm using the Convection(/advection)-Diffusion(-Reaction) equation to calculate the temperature over time in different hydraulic parts like a pipe or a heat exchanger.<br> The flow/convection is <strong>always 1D</strong>, while the diffusion, in this case, heat conduction, <strong>can be 1D, 2D or 3D</strong>.<br> Let's consider an example fully in 1D for the sake of easiness:<br> $$\frac{\delta T}{\delta t} = \frac{\lambda}{c_p \rho}\frac{\delta^2 T}{\delta x^2} + v\frac{\delta T}{\delta x}$$ $$\alpha = \frac{\lambda}{c_p \rho}$$ where $\lambda$ is the heat conductivity, $c_p$ the specific heat capacity, $\rho$ the density and $v$ the velocity.</p> <p>To discretize this equation, I use the <strong>forward-time central-space scheme for the diffusion</strong> (second order central) and the <strong>forward-time forward-space scheme for the convection</strong> (first order upwind), as shown here: $$\frac{T^{n+1} - T^n}{\Delta t} \approx \alpha \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}$$ where the sign $\pm$ in $T^n_{i\pm1}$ is depending on the sign of $v$.</p> <p>Now there are two common cases:</p> <ul> <li>$v \neq 0$ and the resulting Peclet-Number is (in most cases) $|Pe| &gt; 2$</li> <li>$v = 0$ and the resulting Peclet-Number is thus $Pe = 0$ (and the PDE becomes fully parabolic)</li> </ul> <p>Having a relatively small timestep is ok in most cases since time continuous systems like PID controllers interact with the simulation results in between each step. These, for example, control the massflows, internal heat gains etc. Thus I've been using an embedded Heun and/or RK method with adaptive stepsize control to get the explicit time results.</p> <h2>From this multiple questions arise:</h2> <ul> <li>Is using this kind of mixed discretization ok?</li> <li>Should I use central differencing scheme for convection for $|Pe| &lt; 2$</li> <li>Could I implement a QUICK-scheme for the convective term to reduce numeric errors?</li> <li>Is implementing higher accuracy schemes for these derivative orders, as shown in this <a href="https://en.wikipedia.org/wiki/Finite_difference_coefficient" rel="nofollow noreferrer">table</a>, useful?</li> <li>When having mixed kinds of "parts" or parts with a big differences in the size of the grid spacing (and thus of the resulting mass/volume of the cells around each node) this results in the parts with the smallest cells never (or really slow) reaching the temperature of the inflowing 1D massflow (neglecting cell heat loss by conduction). Is this called numeric diffusion or dissipation? And can I reduce this by using other differencing schemes like QUICK?</li> <li>I'd like to solve a few parts with implicit schemes, Crank-Nicholson is preferred. Is Crank-Nicholson space discretization centered scheme for <em>both</em> convection AND diffusion <em>or</em> can I use centered scheme for diffusion and upwind/QUICK for convection?</li> </ul> <p><strong>edit:</strong> Is there any additional information I could provide to make it easier to answer this question?</p> https://scicomp.stackexchange.com/q/29847 2 Finite difference Neumann boundary conditions: uneven weighting of edge nodes? WilliamMorris https://scicomp.stackexchange.com/users/28191 2018-07-13T03:00:19Z 2018-07-15T15:34:59Z <p>Originally asked this on <a href="https://math.stackexchange.com/q/2848196/332683">math.stackexchange</a>, but I figure it's also appropriate here. I'm reading through some finite difference code for a diffusion equation and came across something odd for the boundary conditions, and I was wondering about the validity of the method. I'm familiar with the finite difference method, but I've never seen this certain method before.</p> <p>The equation I am solving is </p> <p>$$\frac{\partial f}{\partial t}=C(x)\frac{\partial^2 f}{\partial x^2}$$</p> <p>with boundary condition </p> <p>$$\frac{\partial f}{\partial x}=0\ \text{on}\ x=0,1$$</p> <p>across time, $t$, and position, $x$, with position-dependent diffusion constant, $C(x)$. It's being solved in an explicit way using the forward difference for time and the second-order central difference for space</p> <p>$$\frac{f^\text{new}_i-f^\text{old}_i}{\Delta t}=C_i\frac{f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}}{(\Delta x)^2}$$</p> <p>along discretised space $\{x_1,\dots,x_n\}$. Rearranging the above equation gives </p> <p>\begin{align} f^\text{new}_i=&amp;f^\text{old}_i + \frac{C_i\Delta t}{(\Delta x)^2}\Big(f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}\Big)\\ =&amp;f^\text{old}_{i-1}(S_i)+f^\text{old}_{i}(1-2S_i)+f^\text{old}_{i+1}(S_i) \end{align}</p> <p>where $S_i=\frac{C_i\Delta t}{(\Delta x)^2}$.</p> <p>This is simply represented in the code. For the boundary at $x_1$ (and at $x_n$, but I'll only talk about the first as they are similar), there are some strange things going on. Usually I would do a Taylor series expansion of function $f$ at points $x_2$ and $x_3$</p> <p>$$f_2=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i \Delta x + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$ $$f_3=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i (2\Delta x) + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (2\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$</p> <p>which, due to boundary condition $\frac{\partial f}{\partial x}=0$ can be rearranged to </p> <p>$$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{2}{(\Delta x)^2}\Big(f_2 - f_1\Big)\tag{A}$$ $$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{2(\Delta x)^2}\Big(f_3 - f_1\Big)\tag{B}.$$</p> <p>Usually, I'd average $A$ and $B$ to determine how the boundary conditions are applied in the eventual matrix solution</p> <p>$$\frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{5}{4}f_1 + f_2 + \frac{1}{4}f_3 \Big).$$</p> <p>However, the code that I'm reading does something different and takes the boundary as a weighted sum of $A/5+4B/5$, as shown below</p> <p>$$\frac{A+4B}{5} \Rightarrow \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{4}{5}f_1 + \frac{2}{5}f_2 + \frac{2}{5}f_3 \Big).$$</p> <p>I don't think this as being inherently wrong, as the LHS is still $\frac{\partial^2 f}{\partial x^2}\big\rvert_i$, but I'm not sure. Is this a valid step to make? Is it ok to "weight" results from different nodes, and if so, in what situations would weighting or not weighting be useful? </p> https://scicomp.stackexchange.com/q/29245 1 Jacobian Elements for Coupled Drift-Diffusion System using Vertex-Centered Finite Volume Chronum https://scicomp.stackexchange.com/users/17947 2018-04-06T01:44:53Z 2018-04-06T01:44:53Z <p>I'm trying to solve the fully coupled drift-diffusion system using Newton's Method. Although I eventually plan to potentially use a Jacobian-Free Newton-Krylov approach, this is still something that I want to implement and look into.</p> <p>A bit of notation for the equations presented below:</p> <p>$d_{ij}$ and $l_{ij}$ are control area/volume constants, not physical entities. The $i$ and $j$ indices are for the central point being considered ($i$), and all neighbouring points that are connected ($j$'s). $A_{\Omega_i}$ is the control area/volume associated to each vertex/node, and is also constant.</p> <p>All quantities that are explicitly required during calculation are at node points. In the equations $B(\Delta_{ij}) \equiv B(\phi_i - \phi_j)$, the Bernoulli function, used as a smoothing method for the differences. This difference should <em>technically</em> exist between node points, but that's for another time.</p> <p>Here are the equations: $\phi$ is electrostatic potential, $n$ and $p$ are carrier concentrations. These are the independent variables to solve for. Everything else is either a constant, or a function of one of these variables. </p> <p>$$\sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) = -\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i}$$</p> <p>$$\frac{\partial{n_i}}{\partial t}A_{\Omega_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$</p> <p>$$\frac{\partial{p_i}}{\partial t}A_{\Omega_i} = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$</p> <p>In residual form, and assuming equilibrium conditions so all time derivatives go to 0, these equations are:</p> <p>$$\Phi = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) +\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i} = 0$$</p> <p>$$N = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$</p> <p>$$P = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$</p> <p>Now, my question is whether I have these Jacobian elements correct. Here is the Jacobian element formula that I've tried to come up with. Only a few equations have been derived, since the remainining elements are just trivial variations of these 3 formulations. This post is mostly about getting more experienced eyes to look at these equations, and potentially see something that I haven't quite seen yet. Are these formulations for the Jacobian elements correct, or is there something wrong with them that I'm missing?</p> <p>$$\frac{\partial\Phi}{\partial\phi_i} = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}} +\frac{q}{\epsilon}(\frac{\partial p_i}{\partial \phi_i} - \frac{\partial n_i}{\partial \phi_i})A_{\Omega_i}$$</p> <p>$$\frac{\partial{N}}{\partial \phi_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \left\{\frac{\partial{n_i}}{\partial \phi_i}B(\Delta_{ij}) + n_i \frac{\partial B(\Delta_{ij})}{\partial \phi_i})\right\}\right] + \frac{\partial(G_i - R_i)}{\partial \phi_i}A_{\Omega_i}$$</p> <p>$$\frac{\partial N}{\partial n_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} B(\Delta_{ij})\right] + \frac{\partial (G_i - R_i)}{\partial n_i}A_{\Omega_i}$$</p> https://scicomp.stackexchange.com/q/29240 2 Closed boundary conditions in finite difference method for diffusive-advective equation Gabriele https://scicomp.stackexchange.com/users/27377 2018-04-04T18:21:26Z 2018-04-05T14:47:44Z <p>I am implementing a finite difference method in solving the diffusive-advective equation: $$u_t + v \cdot u_x = D\cdot u_{xx}$$ (v, D are constants). Planning to use the operator splitting method (see below), I am now <strong>focusing on the advective part</strong>.</p> <p>I chose staggered leapfrog method (from <a href="http://numerical.recipes/" rel="nofollow noreferrer">Numerical Recipes book</a>) $$u^{n+1}_j = u^{n-1}_j - \frac{v\cdot\Delta t }{\Delta x}(u^n_{j+1} - u^n_{j-1}) \;\;\;\;$$ because it should be Courant-stable in case of conserved flux; this turns out to be the case, since it works perfectly using periodic conditions. The problems show when trying to change the boundaries. My question, in brief, is: how to explicitly implement these conditions in such a method? Below, I am going to provide more details on my failings. </p> <h2>No-flux conditions</h2> <p>Let us start from a more general case: as it is pointed out <a href="https://scicomp.stackexchange.com/a/5440/27377">here</a> for diffusion-advection equation, Robin conditions are required to achieve closed boundaries $$v\cdot u - Du_x = 0$$ I tried to discretize them: $$v\cdot u_j - \frac{D}{\Delta x} (c_j - c_{j-1}) = 0$$ However, when D=0 the only condition remaining is that $$v\cdot u = 0$$ thus it is to me unclear how to interpret this: imposing u=0 would result in the (failing) strategy for absorbing conditions (see below), while v is constant (and I don't understand where I could impose it equals to zero). </p> <h2>Absorbing conditions</h2> <p>Leapfrog method is meant for flux-conserving situations. In fact, imposing</p> <pre><code>u(j=0,n) = 0 u(jmax,n) = 0 </code></pre> <p>creates artificial instabilities. The leapfrog updating has got three parts: one for the previous (two steps before) value in the point, one increasing this value proportionally to the uphill point, the other reducing proportionally to the downhill one. This is also what one could expect from a physical point of view: some quantity enters, some goes out. The reason why the simple imposition of null boundaries fails is now evident: if the downhill point is null, there is no reducing term, so the point before the last explodes; the point before it has now a big reducing term, so goes to zero; the point before then has a small reducing term and so on, resulting in an alternance of increasing and decreasing terms (see figure <a href="https://i.stack.imgur.com/WaDj7.png" rel="nofollow noreferrer">1</a>). <a href="https://i.stack.imgur.com/WaDj7.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WaDj7.png" alt="fig.1"></a> </p> <h2>Operator splitting</h2> <p>As explained in the Numerical Recipes book (cited above), an equation like $$u_t = \mathcal{L} u$$ where a decomposition in linear operators $$\mathcal{L}=\sum^m_{i}\mathcal{L}_i$$ is possible can be resolved applying, for each step, the sequence $$\mathcal{U}_1(\mathcal{U}_2(\dots\mathcal{U}_m(u^n)\dots)) = u^{n+1}$$ provided that each operator U solves for a term in the sum and is stable.</p> <hr> <p>Even if I know alternatives are possible, which contemplate diffusion-advection function in general, before digging into them I am willing to solve this doubt of mine.</p> https://scicomp.stackexchange.com/q/29227 3 Finite Element Stabilization for Drift-Diffusion/Advection-Diffusion Equations Chronum https://scicomp.stackexchange.com/users/17947 2018-04-03T02:46:22Z 2018-04-03T02:46:22Z <p>I've tried my best to look through the relevant suggested similar questions when posting this, and hopefully this contains enough new material to not be considered a duplicate.</p> <p>I'm currently trying to draw up a semiconductor drift-diffusion simulation framework, and the work that I'm supposed to be working off of, even though it's called 'finite element' in the section heading, is control-volume finite element, or vertex-centered finite volume (depending on which field you come from, one may be more common notation than the other).</p> <p>As best as I understand it right now, the finite volume method preserve fluxes by construction, and is therefore preferred for implementation of conservation laws. It also allows the Scharfetter-Gummel Finite Box method to be implemented for current densities across the finite volume elements. The SG-FB method is a generalization of the original SG formulation, where the finite differences between carrier densities (the $\nabla n$ term below) are calculated using a special exponential smoothing scheme. This is done because $n$ can often change discontinuously between two points in space, and the SG method is a stable difference scheme.</p> <p>However, I would like to implement this in finite elements, not volume, almost entirely because of the libraries available for FEM (deal.II, libMesh, FEniCS, ...). The relevant current density equations for electrons are below (these, combined with an identical equation for holes, and a nonlinear Poisson equation, make up the system)</p> <p>$$\frac{\partial n}{\partial t} = \frac{1}{q} \nabla \cdot \mathbf{J_n} + G_n - R_n$$</p> <p>$$\mathbf{J_n} = qD_n\nabla n - q\mu_ nn\nabla \phi$$</p> <p>For the questions:</p> <ol> <li><p>Are there any flux/mass-conserving schemes in FEM that have properties similar/identical to the flux/mass-conserving properties of FVM? I keep seeing in a few places that discontinuous-Galerkin conserves flux/mass, but I can't quite get a consensus on it.</p></li> <li><p>What would be appropriate stabilization scheme(s) to apply to this system, given the rapid fluctuations in $n$ (and $p$, identically, for the hole continuity equation)? I keep reading about streamline-upwind Petrov-Galerkin, Galerkin least Squares, and shock-capturing terms, but this seems to be a generally open question.</p></li> <li><p>Given that I keep seeing shock capturing terms in my research, is it reasonable to claim that rapid, discontinuous changes in carrier density are identical to 'static' shock waves/fronts?</p></li> </ol> <p>I'm still quite new to these equations, having dealt with them only for the past few months. I'd be happy to provide any more information, if required.</p> https://scicomp.stackexchange.com/q/29109 1 finite differences on a slanted grid --- advection diffusion equation shamalaia https://scicomp.stackexchange.com/users/27064 2018-03-21T23:33:22Z 2018-03-23T17:26:58Z <p>I used pretty much all my expertise with finite differences to solve an advection-diffusion equation with space-dependent coefficient with a grid in the $x$,$z$ domain with regular spacings. Something like this:</p> <p><a href="https://i.stack.imgur.com/l3TWu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/l3TWu.png" alt="enter image description here"></a></p> <p>This is nice because it is easy to initialize with some experimental measurements that I have. However, near the sloping boundary, I have to carefully choose the grid spacings or sometimes my code behaves funny. This also limits the vertical resolution that I can have very close to the boundary.</p> <p>To solve this problem I was thinking to change reference system by rotating it by an angle $\theta$. Something like this:</p> <p><a href="https://i.stack.imgur.com/qxcSz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qxcSz.png" alt="enter image description here"></a></p> <p>so that I could easily increase the resolution near the boundary. What I do not know is if this imposes some limitations on the physics. Is it sufficient to multiply my equations by a rotation matrix in order to be consistent with the new reference system?</p> <p>Alternatively, is it possible (and fairly easy to do) to add triangular cells near the slope and still use finite difference approximations?</p> https://scicomp.stackexchange.com/q/29005 1 Don't we care about the numerical diffusion in the diffusion term? toliveira https://scicomp.stackexchange.com/users/18932 2018-03-10T04:20:12Z 2018-03-14T12:49:08Z <p>In the context of the solution of advection-diffusion equations by finite volume method, many numerical schemes, papers and book chapters are dedicated to address the numerical diffusion and/or numerical dispersion that comes from the discretization of the advection term.</p> <p>If I understand it correctly, the discretization of the diffusion term also creates numerical diffusion and/or dispersion. However, given the lack of literature about it, it seems not to be a problem.</p> <p>Why is that so?</p> https://scicomp.stackexchange.com/q/27970 2 CFL condition in Stokes equation nabber https://scicomp.stackexchange.com/users/25093 2017-10-02T13:33:07Z 2017-10-25T03:50:07Z <p>Does the CFL condition play any role in a pure Stokes flow, i.e. convective term is neglibile, or vanishing? If not, what is the "equivalent" condition for stability? I have read something about the diffusional time scale but it's quite vague. Anyone with in-depth insight?</p> <p>And a quotation from a paper :<br> <em>"The time step is $3.10^-3\dot\gamma^-1$, corresponding to a CFL number based on the diffusional time scale $CFL_d=v\Delta t/\Delta^2$ of 50."</em> </p> <p> Gallier, S., Lemaire, E., Lobry, L., &amp; Peters, F. (2014). A fictitious domain approach for the simulation of dense suspensions. Journal of Computational Physics, 256, 367-387.</p> https://scicomp.stackexchange.com/q/27633 1 How to simulate 3D diffusion in python? Kama https://scicomp.stackexchange.com/users/25151 2017-08-14T19:09:39Z 2017-08-21T12:17:39Z <p>I want to simulate a simple 3D diffusion (e.g., an ink released from one side of a vessel) using <code>SciPy</code>. There are some tutorials for one-dimensional diffusion. Although the ink goes in one direction, it will not be straight. I appreciate if you can help me to consider the spread of the ink along with other dimensions too.</p> <p>Solving the Fick's second law in one-dimension gives a straightforward equation to draw the concentration profiles against <em>x</em> at different times. How to compute the profile of <em>c(x,y,z,t)</em> at different times?</p> <p>If introducing the code in any other programming language, I can transform it into python. I just want to understand the algorithm structure.</p> https://scicomp.stackexchange.com/q/27625 1 Jacobian of the electron and hole drift diffusion equations with respect to potential in semiconductors P. Biswas https://scicomp.stackexchange.com/users/25146 2017-08-14T08:47:21Z 2017-08-14T17:45:28Z <p>These are the discretized drift diffusion equations as taken from the book "Analysis and Simulation of Semiconductor Devices". The electron continuity equation is:</p> <p>$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T).n_{i+1,j} + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^n.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T).n_{i-1,j} - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}).D_{i-1/2,j}^n.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T))].n_{i,j} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T).n_{i,j+1} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T).n_{i,j-1} = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2)$$</p> <p>The hole continuity equation is:</p> <p>$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T).p_{i+1,j} + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T).p_{i-1,j} - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T)) + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T))].p_{i,j} + \\ ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T).p_{i,j+1} + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T).p_{i,j-1} = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2\ ,$$</p> <p>where $\Delta_i$ and $\Delta_j$ are the mesh spacing for finite difference grid along X and Y directions. $D_n$,$D_p$ are electron and hole diffusion coefficients.$B(x) = x/(e^x-1)$ is Bernoulli's function and $V_T$ is the thermal voltage.</p> <p>Also, $$R^{EFFECTIVE} = R^{SRH} + R^{AU} + R^{OPT} + R^{II}$$ where $R^{EFFECTIVE}$ is the effective generation/recombination term. $R^{SRH}$, $R^{AU}$, $R^{OPT}$ and $R^{II}$ are SRH, Auger, Radiative and Impact Ionization Recombination terms.</p> <p>Now, I am applying the coupled Newton's iterative scheme to solve the above equations where I need to differentiate the equations. When differentiating with respect to potential '$\varphi$', I noticed that $\varphi$ appears in three places namely,</p> <ol> <li><p>High field Electron/Hole Mobility, $\mu_{n,p}^{high}=\mu_{n,p}^{low}⁄ [1 + (\mu_{n,p}^{LICN}.E_{n,p}/v_{n,p}^{sat})^{\beta_{n,p}}]^{1/\beta_{n,p}}$ and $D_{n,p} = \mu_{n,p}.V_T$;</p></li> <li><p>Bernoulli's function, $B(x)$ where 'x' here means difference between potentials at neighboring grid points and</p></li> <li><p>Impact Ionization $R^{II} = -\alpha_n.(|J_{i,j}^n|/q) - \alpha_p.(|J_{i,j}^p|/q)$, $\alpha_n = a_n.e^{-b_n/|E_{i,j}|}$, $\alpha_p = a_p.e^{-b_p/|E_{i,j}|}$ where $J^{n,p}$ and $E$ represent electron/hole currents and electric field respectively. Current terms are made of potential terms as in, $J_{i+1/2,j}^{nx}=D_{i+1/2,j}^n.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).n_{i,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).n_{i+1,j}]/\Delta_i^x \rightarrow (31)$, $J_{i,j+1/2}^{ny}=D_{i,j+1/2}^n.[B(\frac{\varphi_{i,j}-\varphi_{i,j+1}}{V_T}).n_{i,j} - B(\frac{\varphi_{i,j+1}-\varphi_{i,j}}{V_T}).n_{i,j+1}]/\Delta_j^y \rightarrow (32)$, $J_{i+1/2,j}^{px}=D_{i+1/2,j}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i+1,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_i^x \rightarrow (33)$, $J_{i,j+1/2}^{py}=D_{i,j+1/2}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i,j+1} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_j^y \rightarrow (34)$</p></li> </ol> <p>So, basically I need to differentiate all of these to find the Jacobian with respect to $\varphi$. But, the problem arises when the difference between potentials become zero. This means that the entire electron/hole equation becomes a constant when the potential difference is zero. Hence, the Jacobian is also zero. </p> <p>This would mean that the final Jacobian matrix would have entire columns and rows made out of zeros thus making the matrix singular and it is not possible to solve a singular matrix.</p> <p>Is this the correct approach to solving using the Newton's method. What should I do when the potential difference is zero?</p> <h2>References</h2> <ol> <li>Selberherr, Siegfried. Analysis and simulation of semiconductor devices. Springer Science &amp; Business Media, 2012.</li> </ol> https://scicomp.stackexchange.com/q/27389 1 How to formulate Poisson's equation into flux eqution cbcoutinho https://scicomp.stackexchange.com/users/19791 2017-07-12T04:12:27Z 2017-07-12T04:12:27Z <p>I have a small 2D system I'm trying to model using a non-linear extension of Darcy's law for fluid flow in porous media. I'm primarily interested in the local flow velocity, not necessarily the pressure, that's later used in the real model I'm working on:</p> <p>The original equation known as Darcy's Law used in ground water modeling is essentially Poisson's equation, where $\vec{q}$ is the fluid velocity, $\mu$ and $\kappa$ are constants that together describe the hydraulic conductivity of the porous medium, $p$ is local hydrostatic head and f is a source:</p> <p>$$-\frac{\mu}{\kappa} \nabla p = \vec{q}$$</p> <p>$$\nabla \cdot ( \vec{q} ) = f = -\frac{\mu}{\kappa} \nabla^2 p$$</p> <p>The non-linear equation I want to solve instead is the following, known as the Darcy–Forchheimer law. As you can see it's essentially a polynomial instead of a linear model like Darcy's law:</p> <p>$$\frac {\partial p}{\partial x} = \frac {\mu }{\kappa} \vec{q} - \frac {\rho }{\kappa _{1}}\vec{q}^{2}$$</p> <p>I've picked up on some authors  that go through the calculus and have figured out a way to represent $\vec{q}$ as a function of $p$ by calculating the inverse of the flow/pressure equation, and I see somewhat how they arrive at their methodology, but what I'm interested in is actually calculating the flux ($\vec{q}$) itself, which I use elsewhere in the model I'm developing - I don't have much of a use for $p$ besides in the initial problem setup with some simple Dirichlet boundary conditions.</p> <p>Somewhat related, I've also read that for some numerical methods such as DG FEM, diffusion dominated problems can very unstable due to the non-directional flux. In this paper , the authors investigate some 'equation splitting' methods where they solve one equation for the state variable, and another for the flux, essentially the first two equations I showed above. This apparently relieves some of the instability issues caused by diffusion</p> <hr> <p>Looking at these two, it looks like I have two options:</p> <ol> <li>Use the methodology of  to calculate the pressure throughout the system, and then calculate the gradient of the hydrostatic head to generate the local velocity, Or...</li> <li>Split my equation into two equations and solve them together. I end up having the solve two equations at once, but I need to have the velocity regardless so that doesn't seem so bad. Other than that I know nothing about this method.</li> </ol> <p>I don't know how to implement the equation splitting method as I've never done it before, so I have two questions:</p> <ol> <li>Is this even a good idea? Is splitting the equations up into double the number of equations disadvantageous for any other reason than doubling the number of state equations to solve?</li> <li>Is there a source that goes into the actual implementation of equation splitting for (preferably) the finite element or finite difference methods. I am thoroughly intimidated by any DG theory I come across, and I'm hoping there is an easier way to solve this problem</li> </ol> <p> <a href="https://arxiv.org/pdf/1508.00294.pdf" rel="nofollow noreferrer">https://arxiv.org/pdf/1508.00294.pdf</a></p> <p> <a href="https://www.cs.utah.edu/~kirby/Publications/Kirby-11.pdf" rel="nofollow noreferrer">https://www.cs.utah.edu/~kirby/Publications/Kirby-11.pdf</a></p> https://scicomp.stackexchange.com/q/27282 1 Solve ODE with two unknown functions Alessandro Simon https://scicomp.stackexchange.com/users/24727 2017-06-30T08:31:36Z 2017-06-30T08:31:36Z <p>I want to solve a diffusion-convection problem numerically. I start from a PDE of the form $$\partial_t f(x,t) = \partial_x(K(x)\;f(x,t))+\partial_{xx}\;f(x,t)^{m}\;.$$ It is possible to calculate a stationary solution $f_{st}(x)$, which looks something like this $$f_{st}(x) = (a \cosh(x-b))^c\;.$$ Now I try to simplify the original PDE be assuming that the solution is similar to the stationary one, the only difference being that the coefficients $a$ and $b$ are now time dependent, i.e. $a(t)$ and $b(t)$. This allows me to calculate the RHS of the PDE, reducing it formally to a first order ODE in time. My problem is that I'm not sure how to feed it to the computer, because due to the chain rule the spatial derivatives 'extract' additional instances of $a(t)$ and $b(t)$ out of $f$. In other words, I end up with a ODE like $$\partial_t(g(t)\;h(u(t)) = f(g(t),u(t))$$ Intuitively I'd expect there to be a problem because I have two unknowns and just one equation but I'm trying to retrace a paper that seems to be doing just that. Anyways, thanks for reading, hope you can understand what I mean. </p> https://scicomp.stackexchange.com/q/27178 15 BDF vs implicit Runge Kutta time stepping user107904 https://scicomp.stackexchange.com/users/11604 2017-06-16T22:38:13Z 2017-06-17T16:39:52Z <p>Are there any reasons for why one should choose high order implicit Runge Kutta (IMRK) over BDF time stepping? BDF seems much easier to me as $q$ stage IMRK needs $q$ linear solves per time step. Stability for BDF and IMRK looks to be a moot point. I can't find any resources comparing/contrasting implicit time steppers.</p> <p>If it helps, the end goal is for me to select a high order implicit time stepper for advection-diffusion PDE.</p> https://scicomp.stackexchange.com/q/27097 1 Can this nonlinear advection-diffusion equation be discretized as to only have to solve SPD systems? eepperly16 https://scicomp.stackexchange.com/users/16757 2017-06-11T01:16:12Z 2017-07-11T17:39:23Z <p>Consider a nonlinear advection-diffusion equation of the form</p> <p>$$\frac{\partial u}{\partial t} = \nabla \cdot (a(u) \nabla b(u) - \vec{c}(u)u) \tag{1}$$</p> <p>on a rectangular domain with Dirichlet boundary conditions. We have $a(u) &gt; 0$ for all $u \in \mathbb{R}$ and $b(u)$ is a monotonically increasing function of $u$. I wish to solve (1) by finite differences.</p> <p>Can (1) always be discretized in such a way such that all matrix solves involve only symmetric positive definite matrices? If not, are there conditions on $a$, $b$, and $\vec{c}$ such that this is the case?</p> https://scicomp.stackexchange.com/q/26192 1 When is it safe to ignore the diffusion term in an advection-diffusion equation? balborian https://scicomp.stackexchange.com/users/7023 2017-02-16T03:28:45Z 2017-02-16T18:29:15Z <p>Given the one dimensional equation:</p> <p>$\epsilon\frac{\partial^2u}{\partial x^2} +\frac{\partial u}{\partial x} = 0$ </p> <p>with $0\le\epsilon \ll1$ with boundary conditions $u(0) = 0$ and $u(1) = 2$, we can't neglect the diffusive term because of the boundary conditions. In which situations (with very small $\epsilon$ could we? What's the mathematical theory behind that could explain it?</p> <p>Also, if we had a time-dependent equation:</p> <p>$\frac{\partial u}{\partial t} = \epsilon\frac{\partial^2u}{\partial x^2} +\frac{\partial u}{\partial x}$ </p> <p>how would the situation change? what if the convective term were nonlinear, such as in the Burgers equation? I haven't been able to find references on this topic, I would appreciate any. </p> https://scicomp.stackexchange.com/q/26016 1 How can I numericaly solve a convection-diffusion equation with a large diffusion term? fabian https://scicomp.stackexchange.com/users/14651 2017-01-25T09:57:37Z 2017-01-25T17:43:28Z <p>I want to numerically solve the advection-diffusion equation: \begin{equation} u_t(x,t) + cu_x(x,t) = v u_{xx}(x,t) \end{equation} for $x \in [0,1]$ and $t \geq 0$ subject to the boundary conditions $u(0,t) = u(1,t) = 0$ and $u(x,0) = f(x)$. To compare the accuracy of the numerical solution, I will at first derive the analytical solution.</p> <p><strong>Analytical Solution</strong> The solution satisfies \begin{equation} u(x,t) = h(x,t) \exp \left( \alpha x + \beta t \right) \end{equation} If one sets $\alpha = c /2v$ and $\beta = - c^{2} / 4v$ the function $h(x,t)$ adheres to the conventional heat equation: $h_t(x,t) = vh_{xx}(x,t)$ suject to $h(0,t) = h(1,t) = 0$ and $h(x,0) = g(x)$. The solution is particularly inexpensive to calculate if one sets $g(x) = \sin ( 2\pi x)$. Then, the function $u$ is: \begin{equation} u(x,t) = \sin (2 \pi x) \exp \left( -4 \pi^2 v t - \dfrac{c^2 t}{4v} + \dfrac{cx}{2v} \right ) \end{equation}</p> <p><strong>Numerical Solution with BCTS and Crank-Nicolson</strong> To solve the advection-diffusion equation numerically, I use the BTCS and Crank-Nicolson algorithms. Define first \begin{align*} r &amp;= \dfrac{v \Delta t}{\Delta x^{2}} \\ R &amp;= \dfrac{c \Delta t}{\Delta x} \end{align*} The discretized version of the advection-diffusion equation adherring to the BTCS algorithm is: \begin{equation} u_{k}^{n}(1 + 2r) + u_{k-1}^{n} ( -R/2 - r) + u_{k+1}^{n} (R/2 -r) = u_{k}^{n-1} \end{equation} The differential equation for the Crank-Nicolson algorithm can be written as: \begin{align*} &amp;u_{k}^{n} ( 1 + r) + u_{k+1}^{n} ( R/4 - r/2) + u_{k-1}^{n} ( -R/4 - r/2) = \\ &amp; u_{k}^{n-1} ( 1 - r) + u_{k+1}^{n} ( -R/4 + r/2) + u_{k-1}^{n} ( R/4 + r/2) \end{align*}</p> <p><strong>Simulations</strong> When I solve the advection-diffusion algorithm with $c=1, v=2$ and $\Delta t = 0.001$ and $\Delta x \approx 0.0416$ the solution looks as follows: <a href="https://i.stack.imgur.com/clwJb.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/clwJb.png" alt="enter image description here"></a></p> <p>The analytical solution is both positive and negative, the numerical solution however is an inverted parabola. With a reduced strength of the diffusion term $v=1/6$ the solution is very accurate:</p> <p><a href="https://i.stack.imgur.com/SVe5Z.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/SVe5Z.png" alt="enter image description here"></a></p> <p>How can I compute the numerical solution accurately over a large range of values for $v$?</p> <p>I am particularly concerned about this question, as I solved the advection-diffusion equation to understand the algorithms and want to apply them to a non-linear PDE of which I don't know the analytical solution. The computer code used for the example is:</p> <pre><code>import numpy as np from scipy import linalg import matplotlib.pyplot as plt class ConvectionDiffusion(object): """ Class to construct solutions to the convection-diffusion equation """ def __init__(self, DELTA_T, M, V, C): """ Parameters: ----------- DELTA_T: scalar(float): The time step Delta_t M: scalar(float): The number of grid-points for x V: scalar(float): The constant for the diffusion term C: scalar(float): The constant for the advection term """ self.DELTA_T = DELTA_T self.M = M self.DELTA_X = 1 / (M - 1) self.xVec = np.linspace(0, 1., num=M) self.C = C self.V = V self.R = C * DELTA_T / self.DELTA_X self.r = V * DELTA_T / self.DELTA_X**2 self.u0 = np.sin(2 * np.pi * self.xVec) * \ np.exp(self.C * self.xVec / (2 * self.V)) def Implicit(self, stencil, RHS, Sol, l_and_u, T): """ Solving a PDE with implicit algorithms through (banded) matrix inversion Parameters: ----------- stencil: fun Function defining for each value u^{n-1} the stencil such that the discretized PDE can be iterated forward RHS: fun Constructs the RHS of A x = b. l_and_u: tupel(int) The number of lower and upper diagonal elements used in stencil to use linalg.solve_banded T: scalar(float): The time at which the solution to the PDE is expressed. Sol: fun Function generating next periods u^{n} from the solution in the interior x and the boundary values Return: ------- initial: array_like(float): The solution to the PDE at time T """ DELTA_T, initial = self.DELTA_T, self.u0 t = 0.0 while t &lt; T: A = stencil(initial) b = RHS(initial) x = linalg.solve_banded(l_and_u, A, b) initial = Sol(x) t = t + DELTA_T return initial def stencilCN(self, initial): """ Stencil for the Crank-Nicolson algorithm; Stencil * u^n = b Parameters: ---------- initial: array_like(float): The previous solution Returns: -------- A: banded_matrix(float): The matrx in baned for to pass it to scipy.linalg.banded_solve """ A = np.zeros((3, len(initial) - 2)) R, r = self.R, self.r A[0, 1:] = R / 4 - r / 2 A[1, :] = 1 + r A[2, :-1] = -R / 4 - r / 2 return A def stencil(self, initial): """ Stencil for the BCTS algorithm, ; Stencil * u^n = u^{n-1} Parameters: ---------- initial: array_like(float): The previous solution Returns: -------- A: banded_matrix(float): The matrx in baned for to pass it to scipy.linalg.banded_solve """ A = np.zeros((3, len(initial) - 2)) R, r = self.R, self.r A[0, 1:] = R / 2 - r # uper diagonal A[1, :] = 1 + 2 * r A[2, :-1] = -R / 2 - r # lower diagonal return A def RHS(self, initial): """ Pepare the RHS for the BCTS algorithm. Parameters: ----------- initial: array_like(float): The previous solution u^{n-1} Returns: -------- b: array_like(float): Only the interior values of u^{n-1} """ return initial[1: -1] def RHSCN(self, initial): """ Pepare the RHS for the Crank-Nicolson algorithm. Parameters: ----------- initial: array_like(float): The previous solution u^{n-1} Returns: -------- b: array_like(float): A weighted sum of previous u^{n-1} values """ R, r = self.R, self.r present = initial[1:-1] * (1 - r) past = initial[:-2] * (R / 4 + r / 2) future = initial[2:] *(-R / 4 + r / 2) b = present + past + future return b def Sol(self, x): """ Add the boundary condition Parameters: x: array_like(float): The solution in the interior Returns: -------- uNew: array_like(float): The entire solution """ uNew = np.zeros(len(x) + 2) uNew[1:-1] = x return uNew def analyticalSol(self, x, T): """ Analytical Solution of the advection-difussion equation Parameters: ----------- x array_like(float): The state space T: scalar(float): The time at which the solution is evaluated Returns: v: array_like(float): The solution for the PDE """ C, V = self.C, self.V exponential = (- 4 * np.pi**2 * V * T - C**2 * T / (4 * V) + C * x[1:-1] / (2 * V) ) interior = np.sin(2 * np.pi * x[1:-1]) * np.exp(exponential) v = self.Sol(interior) return v def ComparisonSol(self, T): """ Calculates the analytical solution and the its numercial approximation according to the BCTS, C-N and Explicit algorithm at time T Parameters: T: scalar(float): The time value Returns: -------- sol: array_like(float): The array with all solutions e: array_like(float): The normalized second error norm of each approximation """ xVec, DELTA_T, DELTA_X = self.xVec, self.DELTA_T, self.DELTA_X sol, E = np.empty((len(self.u0), 3)), np.zeros(2) ## == Anaytical == ## sol[:, 0] = self.analyticalSol(xVec, T) ## == BTCS == ## sol[:, 1] = self.Implicit(self.stencil, self.RHS, self.Sol, (1, 1), T) E = np.linalg.norm(sol[:, 0] - sol[:, 1]) / np.linalg.norm(sol[:, 0]) ## == CN == ## sol[:, 2] = self.Implicit( self.stencilCN, self.RHSCN, self.Sol, (1, 1), T) E = np.linalg.norm(sol[:, 0] - sol[:, 2]) /np.linalg.norm(sol[:, 0]) return sol, E if __name__ == "__main__": Simulation = ConvectionDiffusion(DELTA_T=0.001, M=50, V=1/6, C=1.) y, e = Simulation.ComparisonSol(0.5) x = Simulation.xVec fig, ax1 = plt.subplots() ax1.plot(x, y[:, 0], 'b-', label='Anaytical Solution') ax1.set_xlabel('x') ax1.set_ylabel('Axis Anaytical Solution', color='b') ax1.tick_params('y', colors='b') ax1.legend(loc=2) ax2 = ax1.twinx() ax2.plot(x, y[:, 1], 'r--', label='BTCS ') ax2.plot(x, y[:, 2], 'r:', label='C-N ') ax2.set_ylabel('Axis Numerical Solutions', color='r') ax2.tick_params('y', colors='r') ax2.legend() fig.tight_layout() plt.show() </code></pre> https://scicomp.stackexchange.com/q/25525 1 Implementing initial conditions into the solution domain of a 1-D advection-diffusion equation user22328 https://scicomp.stackexchange.com/users/22328 2016-11-13T02:33:51Z 2016-11-13T19:33:08Z <p>I have the following PDE.</p> <p>$\frac{\partial c}{\partial t} + v\frac{\partial c}{\partial x} - D\frac{\partial^2 c}{\partial x^2} = 0 \\$.</p> <p>I have discretized it such that i now have $\frac{dC}{dt} = aC_{i+1} +bC_{i} +cC_{i-1} \\$ (1)</p> <p>I have been given the following initial conditions:</p> <p>$\\C=C_0\quad\forall \ 0&lt;t&lt;\tau\\$</p> <p>$\\C=0\quad\forall \ t&gt;\tau\\$</p> <p>i have been asked how i would go about: ' describe how to implement the initial conditions in terms of the solution domain and in terms of any extra vectors'. </p> <p>This question is from a degree course and is about injecting a trace of dye.</p> <p>I have discretised the convective component using simple upwind and the diffusive component using central differences. I'm struggling to understand what the question is asking me. </p> <p>I understand that evaluation at i=1 poses a problem in that it requires C at $\ C_0\\$. This is outside of the solution domain and would require a ghost node, but this is outside the solution domain.</p> <p>What are peoples thoughts? </p> <p>For the 1st initial condition i was thinking of ignoring the diffusive component and for the 2nd initial condition ignoring both the diffusive and convective terms since there is no dye left. Then i would have two equations that describe the initial conditions.</p> https://scicomp.stackexchange.com/q/25462 1 Finite-difference form of the reaction-term in the solute transport equation ToNoY https://scicomp.stackexchange.com/users/22244 2016-11-08T06:41:28Z 2016-12-13T21:17:54Z <p>The partial differential equation is a combination of the diffusion plus convective trans­port equations and an adsorption sink. The equation for one-dimensional solute transport model is:</p> <p>$$\frac{\partial C}{\partial t} = D\frac{\partial^ 2C}{\partial^ 2x}-v\frac{\partial C}{\partial x} - \frac{\rho}{\theta} \frac{\partial S}{\partial t}$$</p> <p>where, C = solute concentration, D = dispersion coefficient, v = average pore-water velocity, x = distance from the inflow position, and t = time. Assuming the adsorption process is a first order reversible reaction, the rate of mass transfer to the adsorbed phase, $\frac{\partial S}{\partial t} = \frac{k_{A}\theta C}{\rho}-k_{D}S$; where, $k_{A}$ and $k_{D}$ are the adsorption (forward) and desorption (backward) rate coefficients (unit: 1/time), $\theta$ is the soil-water content by volume, and $\rho$ is the bulk density of the soil system.</p> <p>The fully explicit finite-difference approximation for all except for the first order reversible reaction term can be written simply as (also, tested to work fine against <a href="https://scicomp.stackexchange.com/questions/25441/numerical-approximation-for-a-known-exact-solution-of-advection-dispersion-equat">exact solution</a>):</p> <p>$$C_{x,t} = C_{x,t-\Delta t} + \frac{D \Delta t}{\Delta x^2}(C|_{x+\Delta x, t-\Delta t}- 2C|_{x,t-\Delta t} + C|_{x-\Delta x,t-\Delta t}) - \frac{v\Delta t}{2 \Delta x}(C|_{x+\Delta x,t-\Delta t} - C|_{x-\Delta x, t-\Delta t})$$</p> <p>I cannot seem to figure out how the above finite-difference approximation could be modified to incorporate the reaction-term defined above. <strong>Hint:</strong> <em>Page#96-99 of the <a href="https://books.google.com/books?id=IYThyxERqMgC&amp;pg=PA94&amp;lpg=PA94&amp;dq=finite%20difference%20form%20of%20first%20order%20reactions%20soil%20column&amp;source=bl&amp;ots=M_0iOGLH_y&amp;sig=-EIXYY2lfyilWddtkvx2QENhkpg&amp;hl=en&amp;sa=X&amp;ved=0ahUKEwjLh-z1ipjQAhVHxGMKHXjJAgc4ChDoAQgdMAE#v=onepage&amp;q=finite%20difference%20form%20of%20first%20order%20reactions%20soil%20column&amp;f=false" rel="nofollow noreferrer">this book</a> does provide a solution</em> but I just cannot get my head around it. I'm supplying the best known articles for the <a href="https://drive.google.com/file/d/0B6GUNg-8d30vOE12SURKb0FMUE0/view?usp=sharing" rel="nofollow noreferrer">analytical solution</a> and <a href="https://drive.google.com/file/d/0B6GUNg-8d30vWlRDLXRJS0ZLb1U/view?usp=sharing" rel="nofollow noreferrer">numerical solution</a> that I could find. Any help with reproducible example codes would be hihgly appreciated.</p> https://scicomp.stackexchange.com/q/25441 3 Numerical approximation for a known exact solution of advection-dispersion equation ToNoY https://scicomp.stackexchange.com/users/22244 2016-11-06T04:24:27Z 2016-11-08T13:22:10Z <p>My goal is to create a numerical solution of 1D-solute transport (Convective-dispersion equation, CDE) to match it's analytical solution based on experimental data. The CDE can be written as (where, C= solute concentration, D= dispersion coefficient, v=velocity, t=time):</p> <p>$$-v\frac{\partial C}{\partial x} + D\frac{\partial^ 2C}{\partial^ 2x} = \frac{\partial C}{\partial t}$$</p> <p>The fully explicit finite-difference approximation for the CDE is:</p> <p>$$C_{x,t} = C_{x,t-\Delta t} + \frac{D \Delta t}{\Delta x^2}(C|_{x+\Delta x, t-\Delta t}- 2C|_{x,t-\Delta t} + C|_{x-\Delta x,t-\Delta t}) - \frac{v\Delta t}{2 \Delta x}(C|_{x+\Delta x,t-\Delta t} - C|_{x-\Delta x, t-\Delta t})$$</p> <p>when, the initial condition: $$C(x,0) = 0$$ Left boundary condition: $$C(0,t) = C_{0}$$ Right boundary condition at infinite distance: $$\frac{\partial C}{\partial x} (\infty, t) = 0$$ </p> <p>the analytical solution can be given by: $$\frac{C}{C_{0}} = \frac{1}{2}\operatorname{erfc}\Bigg[\frac{Rx-vt}{2(DRt)^\frac{1}{2}}\Bigg]+ \frac{1}{2}\exp\Bigg(\frac{vx}{D}\Bigg)\operatorname{erfc}\Bigg[\frac{Rx+vt}{2(DRt)^\frac{1}{2}}\Bigg]$$</p> <p>For the analytical expression, breakthrough $(C=C_{0})$ occurs at a distance of $x=0.075$ $m$ at $0.61$ $days$ (see black line in the attcahed figure) given a $D = 0.0033$ $m^2/day$, $v =0.66$ $m/day$, a retardation, $R = 1$ (dimensionless), and a constant input concentration $C_{0} = 320$ $ppt$. The numerical solution plotted for the same distance, $x$, was simulated with the same $D$ and $v$ values with a $\Delta t = 0.005$ $days$ and $\Delta x = 0.009375$ $m$ (magenta). The parameters used for numerical solution were well within the prescribed stability limits: $\frac{v\Delta t}{2D}=0.00001 &lt;&lt;1$; $\frac{D\Delta t}{\Delta x^2}=0.19 &lt;0.5$; $P_{e}=$$\frac{v\Delta x}{2D}=0.94 &lt;=1; and 1-\Delta t\Big(\frac{2D}{\Delta x^2}+\frac{v}{2\Delta x}\Big)=0.45&gt;=0 along with a Reynold's number, R_{L}=15. The numerical solution was set up in an excel spreadsheet with the same boundary conditions described above. Since it is not possible to set up the right boundary at "\infty" distance, I set it up at an arbitrary distance of x=1.74 m, but extending this distance futhur doesn't change the outcome (e.g. early breakthrough) in any significant way.</p> <p><a href="https://i.stack.imgur.com/alPwb.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/alPwb.jpg" alt="enter image description here"></a></p> <p>My question is: is the numerical approximation good enough? Is there a better way of implementing it especially when retardation and reaction terms needs to be incorporated?</p> https://scicomp.stackexchange.com/q/25029 1 Finite Difference advection-reaction-diffusion in spherical coordinates, problem with Diffusion Gesetzt https://scicomp.stackexchange.com/users/21748 2016-09-21T17:31:53Z 2016-09-26T07:11:49Z <p>I have a problem with the following equation: </p> <p>$$ V_r\frac{\partial C_A}{\partial r} + \frac{V_\theta}{r}\frac{\partial C_A}{\partial \theta} = \frac{2}{Pe_A} \left[ \frac{\partial ^2 C_A}{\partial r ^2} + \frac{2}{r} \frac{\partial C_A}{\partial r} - k C_A\right] $$</p> <p>I'm trying to solve this with finite difference technique. My discretization :</p> <p>$$ V_r \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} + \frac{V_\theta}{r} \frac{A_{i,j+1}-A_{i,j-1}}{2\Delta \theta} - \frac{2}{Pe} \left [ \frac{A_{i+1,j}-2A_{i,j}+A_{i-1,j}}{\Delta r ^2}+\frac{2}{r} \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} - kA_{i,j}\right] = 0$$</p> <p>I first tried it with this, but because I used central-difference in doesn't work for$Pe(cell)&gt;2$, I had strange oscillations at the boundary. So I thought I use a upwind scheme now. </p> <p>$V_r$and$V_\theta$are given through polynoms as a function ($r, \theta$).</p> <p>So for$\theta&lt;90 °$I Use second order upwind for the convection term$V_r &gt; 0$.</p> <p>And for$\theta &gt;90 °$in the other direction.$V_\theta$is always$&lt;0$, so I only use this form for this.</p> <p>My problem is, for high Peclet number and low reaction constants, the concentration is rising up where it shouldn't. This happens in a region, where the convection term gets stronger.</p> <p>What am I supposed to do ? Is there a difference about using upwind in polar coordinates ?</p> <p>Edit :</p> <p>I want go more in detail. This is my coordinate system :</p> <p>The problem appears at the zero degress line. I also have seen that the problem is strongly dependent on$V_\theta$. If I multiply$V_\theta$with 0.9 for example there is no problem. The concentration is not rising. I've looked at the speed vectors in more detail. The problem is vanishing when the speed vectors near the zero degree line are horizontal or show away from it. I have Neumann boundary condition there ( also at 180°). ($\frac{\partial C_A}{\partial \theta}=0 $). Now I also have seen that for the regions where this "pseudo-source" exist the boundary condition isn't met perfectly. There is a gradient. So I am thinking maybe this acts like a "pseudo-source ". I also have drichilet boundary condition for r=1 and r-->infinite. ( My calculation begins at r=1).</p> <p>My grid number is the following ( for example I have three radial and five angular steps ). My notation is c(i,j), i is the radial direction and j the angular direction. So say for i=1, it is c(1,j) with j=0...4 from phi=0...180°. Let's have a look for zero degress with constants l_i :$c10*l1+...c12*l2+c14*l4=0$.</p> <p>I need no ghoistpoints here. I just sad c12=c14 and rearranged the equation to$c10*l1...+c12(l2+l4)\$. For 180 ° I need ghoistpoints because of second order upwind scheme I have c(i,j-2) which is outside my domain, but at 180 ° there is no problem.</p> <p>Should I implement Neumann on another way ? Or do you have any other suggestions ?</p> <p><a href="https://i.stack.imgur.com/3VCYE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3VCYE.png" alt="enter image description here"></a></p> https://scicomp.stackexchange.com/q/24982 3 How numerical diffusion is related to advection term? Ather Cheema https://scicomp.stackexchange.com/users/20089 2016-09-15T09:36:56Z 2016-09-15T10:03:17Z <p>I have crude idea that numerical diffusion arises while using upwind scheme and causes solution to deviate from its original one. But I am unable to understand how numerical diffusion phenomenon is (directly) related to advection phenomenon?</p>