I'm dealing with a coupled system of three transient, non-linear convection-diffusion equations. Let's just say to simplify the problem that they take the following form: $$ -\nabla\cdot(D_{1}(u_{2},u_{3})\nabla u_{1}) = -\nabla\cdot\mathbf{\chi}(u_{2},u_{3}) $$ $$ \frac{\partial u_{2}}{\partial t} + \nabla\cdot f_{1}(u_{1},u_{2},u_{3}) - \nabla\cdot(D_{2}(u_{2},u_{3})u_{2}) = 0 $$ $$ \frac{\partial u_{3}}{\partial t} + \nabla\cdot f_{2}(u_{1},u_{2},u_{3}) - \nabla\cdot(D_{3}(u_{2},u_{3})u_{3}) = 0 $$
Assume a general set of boundary conditions for each equation: $$ u_{i} = g_{i,D}, \;\;\;\text{on}\;\Gamma_{D_{i}} $$ $$ D_{i}\nabla u_{i} = g_{i,N}, \;\;\;\text{on}\;\Gamma_{N_{i}} $$
My spatial discretization is via DGFEM and the time discretization is just Backward Euler. Anyway, I'll use the superscript $(k)$ to denote the solution at the $k$th time step. This is how I am decoupling the equations:
1) Solve for $u_{1}^{(k+1)}$ using information from the previous time step for $u_{2}$ and $u_{3}$. This effectively linearizes the equation.
$$ -\nabla\cdot(D_{1}(u_{2}^{(k)},u_{3}^{(k)})\nabla u_{1}^{(k+1)}) = -\nabla\cdot\chi(u_{2}^{(k)},u_{3}^{(k)}) $$
2) Use the new $u_{1}^{(k+1)}$ and information from the previous time step for $u_{3}$ to solve for $u_{2}^{(k+1)}$:
$$ \frac{\partial u_{2}}{\partial t} + \nabla\cdot f_{1}(u_{1}^{(k+1)},u_{2}^{(k+1)},u_{3}^{(k)}) - \nabla\cdot(D_{2}(u_{2}^{(k+1)},u_{3}^{(k)})\nabla u_{2}^{(k+1)}) = 0 $$
3) Solve for $u_{3}^{(k+1)}$:
$$ \frac{\partial u_{3}}{\partial t} + \nabla\cdot f_{2}(u_{1}^{(k+1)},u_{2}^{(k+1)},u_{3}^{(k+1)}) - \nabla\cdot(D_{3}(u_{2}^{(k+1)},u_{3}^{(k+1)})\nabla u_{3}^{(k+1)}) = 0 $$
So the first equation is linear and the last two are non-linear. I am using Newton's method to handle the non-linearity. However, there are some bugs in my code and I'm a bit stuck trying to interpret the results.
I know this is long but I'm trying to be as detailed as possible here. Just a couple of questions:
1) I test the simple case where $u_{1} = u_{2} = u_{3} = 5+t$ and I get about machine precision error using linear basis functions. Newton's method is also converging quadratically. I'm using a timestep approximately equal to $h^{p+1}$. I can also increase it to $100h^{p+1}$ which surprisingly doesn't seem to affect the convergence much. Using this time step, I have:
/* sample Newton iterations */
******** Current Time: 7.812500 ********
u2 Iteration: 1, Error: 4.155186e-01
u2 Iteration: 2, Error: 2.480219e-02
u2 Iteration: 3, Error: 9.169256e-05
u2 Iteration: 4, Error: 1.250692e-09
u2 Iteration: 5, Error: 8.042476e-16
Iterations for u2 convergence: 5
u3 Iteration: 1, Error: 4.155186e-01
u3 Iteration: 2, Error: 2.480219e-02
u3 Iteration: 3, Error: 9.169256e-05
u3 Iteration: 4, Error: 1.250692e-09
u3 Iteration: 5, Error: 7.664044e-16
Iterations for u3 convergence: 5
/* error results */
*********** L2 Error ***********
u1 error in L2-norm:
4.53162e-12
4.64396e-12
6.15316e-12
1.07398e-11
u2 error in L2-norm:
4.12483e-13
7.30247e-13
1.39962e-12
2.76142e-12
u3 error in L2-norm:
4.09102e-13
7.29944e-13
1.39935e-12
2.76121e-12
*********** H1 Error ***********
u1 error in H1-norm:
1.05844e-11
2.03320e-11
3.84458e-11
7.47941e-11
u2 error in H1-norm:
6.40033e-13
1.33164e-12
2.71254e-12
5.47814e-12
u3 error in H1-norm:
6.34670e-13
1.33089e-12
2.71211e-12
5.47862e-12
Now the bizarre part is that if I take $u_{1}=5+t$, $u_{2}=6+t$, and $u_{3}=7+t$, my results are awful. I'm thinking maybe the first case was the exception rather than the norm since it's the only example I've tested so far that seems to work. Anyway, the weird part about this one is that I'm getting about machine precision error for $u_{2}$ while the others are flatlining around 1.e0. I'm not sure how this could possible be since the solution to $u_{2}$ depends on the solutions to $u_{1}$ and $u_{3}$, which are incorrect if my code is to be believed. I observe quadratic convergence for $u_{2}$, but obviously not the others...
The time step is around $h^{p+1}$ here, but I observe similar results for smaller time steps.
******** Current Time: 0.074219 ********
u2 Iteration: 1, Error: 3.908065e-03
u2 Iteration: 2, Error: 2.219552e-06
u2 Iteration: 3, Error: 2.889154e-12
Iterations for u2 convergence: 3
u3 Iteration: 1, Error: 1.373826e-01
u3 Iteration: 2, Error: 1.564705e-02
u3 Iteration: 3, Error: 2.492157e-03
u3 Iteration: 4, Error: 7.872758e-04
u3 Iteration: 5, Error: 1.311116e-04
u3 Iteration: 6, Error: 2.997904e-05
u3 Iteration: 7, Error: 9.158730e-06
u3 Iteration: 8, Error: 1.519301e-06
u3 Iteration: 9, Error: 3.555316e-07
u3 Iteration: 10, Error: 1.067623e-07
u3 Iteration: 11, Error: 1.755777e-08
u3 Iteration: 12, Error: 4.195674e-09
u3 Iteration: 13, Error: 1.242190e-09
u3 Iteration: 14, Error: 2.024999e-10
u3 Iteration: 15, Error: 4.936768e-11
u3 Iteration: 16, Error: 1.443518e-11
u3 Iteration: 17, Error: 2.333785e-12
Iterations for u3 convergence: 17
*********** L2 Error ***********
u1 error in L2-norm:
3.02951e+00
5.43305e+00
5.96871e+00
4.25056e+00
u2 error in L2-norm:
1.45749e-13
2.71039e-13
3.69646e-13
4.95591e-13
u3 error in L2-norm:
6.29385e-01
8.30180e-01
1.11972e+00
1.13907e+00
*********** H1 Error ***********
u1 error in H1-norm:
7.86736e+00
1.81562e+01
3.67674e+01
3.97038e+01
u2 error in H1-norm:
2.33954e-13
2.00202e-12
4.95659e-12
1.03835e-11
u3 error in H1-norm:
2.00438e+00
3.91391e+00
6.24097e+00
7.06979e+00
I think maybe the most telling part is the convergence of the Newton iterations. Namely, what does it mean if I'm not observing quadratic convergence? My initial guess is just the solution at the previous time step, and since the time step is quite small, I'm pretty sure the initial guess is not the issue. I know that non-quadratic convergence can imply that the Jacobian is not being computed properly, but the first example is so similar to the second that I can't imagine that's the case...anyone have any insights here?
2) I'm a bit iffy on the spatial discretizaton of the convection term: $$ \int_{\Omega}\nabla\cdot f_{1}(u_{1},u_{2},u_{3})z = -\int_{\Omega}\nabla z\cdot f_{1}(u_{1},u_{2},u_{3}) + \int_{\Gamma}z(f_{1}\cdot\mathbf{n}) $$
(via Green's identity). The DG discretization is $$ -\sum_{E\in\mathcal{M}_{h}}\int_{E}\nabla z\cdot f_{1} + \sum_{e\in\Gamma_{\text{int}}}\int_{e}\left\{\left\{f_{1}\cdot\mathbf{n}\right\}\right\}[[z]] + \sum_{e\in\Gamma_{\text{out}}}\int_{e}(f_{1}\cdot\mathbf{n})z + \sum_{e\in\Gamma_{\text{in}}}\int_{e}(f_{1}\cdot\mathbf{n})z $$
where {{ }} is the average and [[ ]] is the jump. In all of the literature I've come across, the inflow integral is moved to the right-hand side. So if we're considering the second equation (solving for $u_{2}$), is this correct?
$$ \sum_{e\in\Gamma_{\text{in}}}\int_{e}(f_{1}(u_{1},g_{2,D},u_{3})\cdot\mathbf{n})z $$
So when computing the Jacobian for the second equation, does this inflow term not contribute to the Jacobian at all since $u_{2,D}$ is known?
I'm able to solve a single transient, non-linear convection diffusion equation using Newton's method, so I'm not entirely sure what it is that I'm missing.
