# Effect of mesh size on solution curves for a 1D problem

I'm interested in studying the effect of mesh size on the behavior of the solution curves of 1D convection-diffusion problem.

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2} - v\frac{\partial C}{\partial x}$$

I could use MATLAB's pdepe solver to solve the above system.

function sol=so()
format short
global D nnode init_co find_index v
m = 0;
delx = 0.25;
xend = 5;
D = 500;
v = 200;
x = 0:delx:xend; % mesh
find_index  = 0:delx:xend;
t = 0:0.00001:0.5;
init_co = [1 ; zeros(length(x)-1,1)];
nnode = length(x);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
plot(t,sol)
xlabel('time')
ylabel('c')
grid on
function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = -v*DcDx;
end

function c0 = icfun(x)
c0 = init_co(find(find_index==x));
end

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
pl = cl - 3;
ql = 0;
pr = 0;
qr = 1;
end
end


I changed delx = 1, the following is observed

By just observing the above plots, I find the behavior (i.e the steepness of the curves) is the same and doesn't vary much with mesh size.

I think it's a crude way to observe just the steepness of the curves. I'd like to ask for suggestions on other quantitative measures that can be computed from solution curves to study the effect of the mesh size.

EDIT: I've compared one quantity(absolute error) from the suggestions received below.

I understand the error increases with an increase in mesh size. I'd like to know how to strike a balance -- fine mesh will increase computational cost and coarse mesh will produce errors and even lead to oscillations in solution curves for high delx. For instance, in the literature, I came across the CFL number, which should be less than 1 to obtain stable solutions. However, I'm not sure how to compute the CFL number for my system since I don't know the step size that the solver(ode15s) takes. Also, I'd like to know what is the acceptable range of error (what percent of deviation from true solution is acceptable? is there a thumb rule?).

• What do the different curves mean on each plot? I couldn't glance it from the code. Feb 6, 2020 at 18:35
• The usual way to study the effect of the mesh size is of course to compute the error in the solution. Feb 6, 2020 at 19:30
• @Javier-Acuna Thanks a lot for the response. The curves denote the transient change in the quantity c observed positions defined in x = 0:delx:xend; % mesh. Feb 7, 2020 at 12:18
• @WolfgangBangerth Thank a lot. I've been looking at absolute and relative errors. Are there instances in which one shouldn't look at absolute/ relative error? For example, if I have 0.5 and 0.2e-25 at point x=1 and time = 0.02, computing(0.5-0.2e-25)/0.2e-25 gives a really high relative error. Feb 7, 2020 at 12:23
• @Natasha: I would have computed the relative error relative to some "typical size" of the solution -- for example, the mean absolute value of the solution. If you make it relative to the exact value at some point, you're always going to get meaningless numbers anywhere the solution is zero or close to zero. Feb 7, 2020 at 17:26

I'm a bit surprised that nobody gave you an analytical solution for comparing your simulation results with it and measure your mesh convergence. It's always good to specify the boundary conditions explicitly, because it's really difficult to understand it from your code unless you dig deep into the MATLAB documentation of pdepe. Your boundary conditions are:

$$C(0,t) = C_{b}$$

and

$$-D \frac{\partial C}{\partial x}|_{x=L} = 0$$

and you have this advection-diffusion equation:

$$\frac{\partial C}{\partial t} = D \frac{\partial^{2} C}{\partial x^{2}} - v \frac{\partial C}{\partial x}$$

Let's define these dimensionless parameters:

$$\theta = \frac{C - C_{b}}{C_{b}}$$

$$\hat{x} = \frac{x}{L}$$

$$\hat{t} = \frac{Dt}{L^{2}}$$

$$Pe = \frac{vL}{D}$$

$$\frac{\partial \theta}{\partial \hat{t}} = \frac{\partial^{2} \theta}{\partial \hat{x}^{2}} - Pe \frac{\partial \theta}{\partial \hat{x}}$$

And your boundary conditions:

$$\theta(0,\hat{t}) = 0$$

$$-\frac{\partial \theta}{\partial \hat{x}}|_{\hat{x} = 1} = 0$$

Your initial condition is:

$$\theta(\hat{x},0) = -1$$

By using separation of variable technique you have:

$$\theta(\hat{x},\hat{t}) = X(\hat{x})T(\hat{t})$$

So:

$$XT^{'} = X^{''}T - PeX^{'}T$$

or:

$$\frac{T^{'}}{T} = \frac{X^{''}}{X} - Pe \frac{X^{'}}{X} = -\lambda^{2}$$

For spatial part:

$$X^{''} - Pe X^{'} + \lambda^{2} X = 0$$

The general solution is:

$$X(\hat{x}) = \exp{(\frac{Pe \cdot \hat{x}}{2})} (A\cos{(\eta \hat{x})} + B \sin{(\eta \hat{x})})$$

If put $$\hat{x} = 0$$:

$$A = 0$$

So:

$$X(\hat{x}) = \exp{(\frac{Pe \cdot \hat{x}}{2})} (B \sin{(\eta \hat{x})})$$

if we put $$\hat{x} = 1$$:

$$X^{'}(\hat{x}) = \frac{Pe}{2} \exp{(\frac{Pe \cdot \hat{x}}{2})} (B \sin{(\eta \hat{x})}) + \exp{(\frac{Pe \cdot \hat{x}}{2})} B \eta \cos{(\eta \hat{x})}$$

$$X^{'}(\hat{x}) = B \exp{(\frac{Pe \cdot \hat{x}}{2})} (\frac{Pe}{2} \sin{(\eta \hat{x})} + \eta \cos{(\eta \hat{x})})$$

So:

$$X^{'}(1) = B \exp{(\frac{Pe}{2})} (\frac{Pe}{2} \sin{(\eta)} + \eta \cos{(\eta)}) = 0$$

So the eigenvalue equation is:

$$\frac{Pe}{2} \sin{(\eta)} + \eta \cos{(\eta)} = 0$$

or:

$$\tan{(\eta_{n})} = -\frac{2 \eta_{n}}{Pe}$$

We know that $$\eta$$ and $$\lambda$$ should be related as:

$$\lambda^{2}_{n} = \frac{Pe^{2}}{4} +\eta^{2}_{n}$$

The temporal solution is:

$$T(\hat{t}) = \exp{(-\lambda^{2}\hat{t})}$$

So:

$$\theta_{n}(\hat{x},\hat{t}) = B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} \exp{\Bigg(-\Big(\frac{Pe^{2}}{4}+\eta_{n}^{2}\Big) \hat{t}\Bigg)}$$

Final solution is:

$$\theta(\hat{x},\hat{t}) = \sum_{n=0}^{\infty} B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} \exp{\Bigg(-\Big(\frac{Pe^{2}}{4}+\eta_{n}^{2}\Big) \hat{t}\Bigg)}$$

$$B_{n}$$s are extracted by using initial condition:

$$-1 = \sum_{n=0}^{\infty} B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})}$$

Or:

$$-\exp{(-\frac{Pe \cdot \hat{x}}{2})} = \sum_{n=0}^{\infty} B_{n} \sin{(\eta_{n}\hat{x})}$$

So:

$$-\exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} = \sum_{n=0}^{\infty} B_{n} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})}$$

Taking integration from $$\hat{x} = 0$$ to $$\hat{x} = 1$$:

$$-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \sum_{n=0}^{\infty} B_{n} \int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x}$$

But:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

Where:

$$\gamma_{n} = \int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d \hat{x}$$

So:

$$-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} d\hat{x} = B_{m}\int_{0}^{1} \sin^{2}{(\eta_{m}\hat{x})} d \hat{x}$$

So:

$$B_{n} = \frac{-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} d\hat{x}}{\int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d \hat{x}}$$

Finally, you can compare your results with different $$\Delta x$$ or mesh size with this analytical solution.

Update:

Just a note here. I was thinking why we should have:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

We know that:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \frac{\cos{(\eta_{n})} \cos{(\eta_{m})} (\eta_{m} \tan{(\eta_{n})} - \eta_{n} \tan{(\eta_{m})})}{(\eta_{n}^{2}- \eta_{m}^{2})}$$

When $$n \neq m$$. But we know that: $$\tan{(\eta_{n})} = -\frac{2\eta_{n}}{Pe}$$ and $$\tan{(\eta_{m})} = -\frac{2\eta_{m}}{Pe}$$:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = -\frac{2}{Pe}\frac{\cos{(\eta_{n})} \cos{(\eta_{m})} (\eta_{m} \eta_{n} - \eta_{n} \eta_{m})}{(\eta_{n}^{2}- \eta_{m}^{2})} = 0$$

So, we must have:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

Where:

$$\gamma_{n} = \int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d\hat{x} = \frac{1}{2} - \frac{\sin{(2\eta_{n})}}{4 \eta_{n}}$$

Finally:

$$B_{n} = -\frac{\eta_{n}}{\lambda_{n}^{2}(\frac{1}{2} - \frac{\sin{(2\eta_{n})}}{4 \eta_{n}})}$$

What you're looking for is an a-posteriori error estimate for your mesh study. Normally, these quantitative measures are line/area/volume- and/or time- averaged nodal or element quantities (e.g. avg. temperature of some body integrated over time). In structural analysis it is also customary to use the energy norm as a quantitative measure.

A line/area or volume-averaged quantity is simply to take the average value across a line, area or a volume. The idea is to avoid that your analysis depends on a single node or element, which might be not representative. In your case, since it's a 1D problem, you may want to integrate $$c(x)$$ from $$x_0$$ to $$x_L$$.

Energy norm is a concept used in structural analysis and thus might not be of applicability to you. You can read more about the energy norm here: https://blogs.solidworks.com/tech/2018/05/manage-solidworks-simulation-mesh-refinement-using-an-energy-norm-error-plot.html

• Thanks a lot for the response. From what I understand, let's say we consider the nodal quantities ( at positions x = 1,2,3). At each position, I've two arrays (stores the transient change in c at x positions), one for delx = 0.25 and other for delx = 1. From what has been suggested in the comment above, I can compute absolute and relative errors. And from your suggestions, I understand the time-averaged nodel quantity can also be compared. Could you please explain a bit about the area-averaged quantity and energy norm? Feb 7, 2020 at 12:14
• I edited my original contribution to answer your questions Feb 8, 2020 at 12:09
• Thank you. Please check my recent edit in the original post Feb 8, 2020 at 16:00