Solves:
\begin{equation} \frac{\partial u}{\partial t} - \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} = 1 \end{equation}with boundary conditions: \begin{equation} u(0) = 0 \\ u(1) = 1 \end{equation}
Changing the nl_rel_tol from 1e-6 to 1e-2 changes the number of time steps to get to 1 second from over 5000 to 781. Problem did not reach steady state by one second, however it did by two seconds. Note that this problem formulation shows undershoots during the initial time steps on the right side of the domain, most likely because we are starting from an initial condition of u(x, 0) = 0 but with a Dirichlet boundary condition of u(1, t) = 1 at the RHS boundary. Perhaps with rDG(P0P1) we would not observe these oscillations; however with DG(P1) (what we are using here), we most certainly do.
Solves:
\begin{equation} \frac{\partial u}{\partial t} - \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} = 1 \end{equation}with boundary conditions: \begin{equation} u(0) = 0 \\ \frac{\partial u}{\partial x}\Bigr|_{x=1} = 0 \end{equation}
The analytic solution to the corresponding steady-state equation is:
\begin{equation} u(x) = -(e - 1)*e^{-1} + x - e^{x - 1} + 1\\ \frac{du}{dx} = -e^{x - 1} + 1\\ \frac{d^2u}{dx^2} = -e^{x-1}\\ - \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} = e^{x-1} - e^{x-1} + 1 = 1 \end{equation}Upon simplification:
\begin{equation} u(x) = x + e^{-1} - e^{x - 1} \end{equation}as determined by sage math in the local notebook solve_1d_advection_diffusion_equation. By inspection, it's clear that the governing equation is satisfied as well as the boundary conditions. So we have a clear, correct analytic solution that we need to try and reproduce with our numerical methods in MOOSE. So far, the numerical solution is tending towards u(x) = x, which is wrong. So the question is to figure out why. Alright now we get a perfect result with DG! And we can run in Steady mode even with the FDP preconditioner. Something is clearly wrong with the user provided Jacobians. That's a debugging issue for tomorrow.
Solves:
\begin{equation} \frac{\partial u}{\partial t} - \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} = 1 \end{equation}with boundary conditions: \begin{equation} u(0) = 0 \\ \frac{\partial u}{\partial x}\Bigr|_{x=1} = 0 \end{equation}
The analytic solution to the corresponding steady-state equation is:
\begin{equation} u(x) = -(e - 1)*e^{-1} + x - e^{x - 1} + 1\\ \frac{du}{dx} = -e^{x - 1} + 1\\ \frac{d^2u}{dx^2} = -e^{x-1}\\ - \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} = e^{x-1} - e^{x-1} + 1 = 1 \end{equation}Upon simplification:
\begin{equation} u(x) = x + e^{-1} - e^{x - 1} \end{equation}as determined by sage math in the local notebook solve_1d_advection_diffusion_equation. By inspection, it's clear that the governing equation is satisfied as well as the boundary conditions. So we have a clear, correct analytic solution that we need to try and reproduce with our numerical methods in MOOSE. So far, the numerical solution is tending towards u(x) = x, which is wrong, just like for the above dg case. Given that, it seems likely that the issue in terms of accuracy is not the dg formulation. I'm tending to believe it must be the right boundary condition because recall that u(x) = x is a possible solution to the second order ODE pending boundary conditions. And the LHS BC is holding true. It's the RHS BC that is being violated. So...the issue was that I had changed my diffusion coefficient to .1, so actually the solve was perfect for CG! Now with diffusion coefficient changed to 1, I am getting the exactly correct analytical solution as predicted by Sage.
In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [3]:
data = np.loadtxt("residual_data.txt")
In [4]:
data
Out[4]:
In [7]:
iteration_num = data[:,0]
residuals = np.log(data[:,1])
plt.plot(iteration_num, residuals)
plt.show()
This is a 2D test case mean to mimic the convection-diffusion properties of the MSR model. There is a fuel region with advection and a source and then there is a moderator region with only diffusion. Boundary conditions are Dirichlet at the bottom foundary (both fuel and moderator) and then an outflow condition at the fuel top. Everywhere else is insulated. Results look good. Able to achieve steady-state for velocity ranging from 1e-3 to 1e3 while diffusion coefficient and source coefficient stay constant at 1. However, there is a negative value somewhere in the domain, which should not be the case. Ok the negative values were on the graphite_bottom portion of the domain. Negative values were observed for monomial basis functions with sigma = 6 and 60. However, with l2_lagrange, negative values were not observed for sigma = 60 or 6. Hmm, might be time to try that tomorrow with the actual MSR temperature simulation...
In [ ]: