From e59363c7a4194be461dddec4be1330b4994fb243 Mon Sep 17 00:00:00 2001 From: Lei Qiao Date: Wed, 18 Feb 2015 12:11:33 -0600 Subject: [PATCH] Crank-Nicolson scheme fixed: The linear interpolation between timestep should be upon flux rather than independent vector. When working at implicit Euler time stepping, the current version generats an output file that identical to output of the base version as expected. The algorithem documentation is also updated while the output part is not. This is because the correction of energy source term make the code diverge much sooner and the result doesn't looks beautiful. --- examples/step-33/doc/intro.dox | 34 +++--- examples/step-33/doc/results.dox | 17 +++ examples/step-33/step-33.cc | 200 +++++++++++++++++++------------ 3 files changed, 162 insertions(+), 89 deletions(-) diff --git a/examples/step-33/doc/intro.dox b/examples/step-33/doc/intro.dox index 86c57e266b..f124cb4000 100644 --- a/examples/step-33/doc/intro.dox +++ b/examples/step-33/doc/intro.dox @@ -128,27 +128,33 @@ Hartmann's PhD thesis ("Adaptive Finite Element Methods for the Compressible Euler Equations", PhD thesis, University of Heidelberg, 2002). We use a time stepping scheme to substitute the time derivative in the -above equations. At each time step, our full discretization is thus +above equations. For simplicity, we define $ \mathbf{B}({\mathbf{w}_{n})(\mathbf z) $ as the spacial residual at time step $n$ : + +@f{eqnarray*} + \mathbf{B}({\mathbf{w}_{n})(\mathbf z) &=& +- \int_{\Omega} \left(\mathbf{F}(\mathbf{w_n}), +\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w_n} , \nabla \mathbf{z}) +\\ +&& ++ +\int_{\partial \Omega} \left(\mathbf{H}(\mathbf{w_n}^+, +\mathbf{w}^-(\mathbf{w_n}^+), \mathbf{n}), \mathbf{z}\right) +- +\int_{\partial \Omega} \left(\mathbf{G}(\mathbf{w_n}), +\mathbf{z}\right) . +@f} + +At each time step, our full discretization is thus that the residual applied to any test function $\mathbf z$ equals zero: @f{eqnarray*} R(\mathbf{W}_{n+1})(\mathbf z) &=& \int_{\Omega} \left(\frac{\mathbf{w}_{n+1} - \mathbf{w}_n}{\delta t}, -\mathbf{z}\right) -- \int_{\Omega} \left(\mathbf{F}(\tilde{\mathbf{w}}), -\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) -\\ -&& \qquad -+ -\int_{\partial \Omega} \left(\mathbf{H}(\tilde{\mathbf{w}}^+), -\mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right) -- -\int_{\partial \Omega} \left(\mathbf{G}(\mathbf{w}), -\mathbf{z}\right) -\\ +\mathbf{z}\right)+ +\theta \mathbf{B}({\mathbf{w}_{n+1}) + (1-\theta) \mathbf{B}({\mathbf{w}_{n}) \\ & = & 0 @f} -where $\tilde{\mathbf{w}} = \theta \mathbf{w}_{n+1} + (1-\theta) \mathbf{w}_n$ for $0 \leq \theta \leq 1$ and +where $ \theta \in [0,1] $ and $\mathbf{w}_i = \sum_k \mathbf{W}_i^k \mathbf{\phi}_k$. Choosing $\theta=0$ results in the explicit (forward) Euler scheme, $\theta=1$ in the stable implicit (backward) Euler scheme, and $\theta=\frac 12$ diff --git a/examples/step-33/doc/results.dox b/examples/step-33/doc/results.dox index 414e8e063e..d1136ec5e6 100644 --- a/examples/step-33/doc/results.dox +++ b/examples/step-33/doc/results.dox @@ -263,6 +263,23 @@ build the Newton matrix each time, the resulting scheme may well be faster. +

Cache the explicit part of residual

+ +The residual calulated in ConservationLaw::assemble_cell_term function +read + $R_i = \left(\frac{\mathbf{w}^{k}_{n+1} - \mathbf{w}_n}{\delta t} + , \mathbf{z}_i \right)_K + + \theta \mathbf{B}({\mathbf{w}^{k}_{n+1})(\mathbf{z}_i)_K + + (1-\theta) \mathbf{B}({\mathbf{w}_{n}) (\mathbf{z}_i)_K $ +means that we calculate the spacial residual twice at one Newton +iteration step: once respect to the current solution ${\mathbf{w}^{k}_{n+1}$ +and another respect to the last time step solution $\mathbf{w}_{n}$ which +remains the same during all Newton interations through one timestep. +Cache up the explicit part of residual + $ \mathbf{B}({\mathbf{w}_{n}) (\mathbf{z}_i)_K}$ +during Newton iteration will save lots of labor. + +

Other conservation laws

Finally, as a direction beyond the immediate solution of the Euler diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 6ae98203de..4987288a96 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -266,16 +266,16 @@ namespace Step33 // numerical flux function to enforce boundary conditions. This routine // is the basic Lax-Friedrich's flux with a stabilization parameter // $\alpha$. It's form has also been given already in the introduction: - template + template static void numerical_normal_flux (const Point &normal, const InputVector &Wplus, const InputVector &Wminus, const double alpha, - Sacado::Fad::DFad (&normal_flux)[n_components]) + number (&normal_flux)[n_components]) { - Sacado::Fad::DFad iflux[n_components][dim]; - Sacado::Fad::DFad oflux[n_components][dim]; + number iflux[n_components][dim]; + number oflux[n_components][dim]; compute_flux_matrix (Wplus, iflux); compute_flux_matrix (Wminus, oflux); @@ -434,7 +434,7 @@ namespace Step33 // component here so that the average of the velocities is // orthogonal to the surface normal. This creates sensitivities of // across the velocity components. - Sacado::Fad::DFad vdotn = 0; + typename DataVector::value_type vdotn = 0; for (unsigned int d = 0; d < dim; d++) { vdotn += Wplus[d]*normal_vector[d]; @@ -1618,25 +1618,31 @@ namespace Step33 // residual, adding its negative to the right hand side vector, and adding // its derivative with respect to the local variables to the Jacobian // (i.e. the Newton matrix). Recall that the cell contributions to the - // residual read $F_i = \left(\frac{\mathbf{w}_{n+1} - \mathbf{w}_n}{\delta - // t},\mathbf{z}_i\right)_K - \left(\mathbf{F}(\tilde{\mathbf{w}}), - // \nabla\mathbf{z}_i\right)_K + h^{\eta}(\nabla \mathbf{w} , \nabla - // \mathbf{z}_i)_K - (\mathbf{G}(\tilde{\mathbf w}), \mathbf{z}_i)_K$ where - // $\tilde{\mathbf w}$ is represented by the variable W_theta, - // $\mathbf{z}_i$ is the $i$th test function, and the scalar product - // $\left(\mathbf{F}(\tilde{\mathbf{w}}), \nabla\mathbf{z}\right)_K$ is + // residual read + // $R_i = \left(\frac{\mathbf{w}^{k}_{n+1} - \mathbf{w}_n}{\delta t} , + // \mathbf{z}_i \right)_K $ $ + + // \theta \mathbf{B}({\mathbf{w}^{k}_{n+1})(\mathbf{z}_i)_K $ $ + + // (1-\theta) \mathbf{B}({\mathbf{w}_{n}) (\mathbf{z}_i)_K $ where + // $\mathbf{B}({\mathbf{w})(\mathbf{z}_i)_K = + // - \left(\mathbf{F}(\mathbf{w}),\nabla\mathbf{z}_i\right)_K $ $ + // + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}_i)_K $ $ + // - (\mathbf{G}(\mathbf {w}), \mathbf{z}_i)_K $ for both + // ${\mathbf{w} = \mathbf{w}^k_{n+1}$ and ${\mathbf{w} = \mathbf{w}_{n}}$ , + // $\mathbf{z}_i$ is the $i$th vector valued test function. + // Furthermore, the scalar product + // $\left(\mathbf{F}(\mathbf{w}), \nabla\mathbf{z}_i\right)_K$ is // understood as $\int_K \sum_{c=1}^{\text{n\_components}} - // \sum_{d=1}^{\text{dim}} \mathbf{F}(\tilde{\mathbf{w}})_{cd} - // \frac{\partial z_c}{x_d}$. + // \sum_{d=1}^{\text{dim}} \mathbf{F}(\mathbf{w})_{cd} + // \frac{\partial z^c_i}{x_d}$ where $z^c_i$ is the $c$th component of + // the $i$th test function. + // // // At the top of this function, we do the usual housekeeping in terms of // allocating some local variables that we will need later. In particular, // we will allocate variables that will hold the values of the current // solution $W_{n+1}^k$ after the $k$th Newton iteration (variable - // W), the previous time step's solution $W_{n}$ (variable - // W_old), as well as the linear combination $\theta W_{n+1}^k - // + (1-\theta)W_n$ that results from choosing different time stepping - // schemes (variable W_theta). + // W) and the previous time step's solution $W_{n}$ (variable + // W_old). // // In addition to these, we need the gradients of the current variables. It // is a bit of a shame that we have to compute these; we almost don't. The @@ -1655,7 +1661,7 @@ namespace Step33 // Table class also supports. // // Secondly, we want to use automatic differentiation. To this end, we use - // the Sacado::Fad::DFad template for everything that is a computed from the + // the Sacado::Fad::DFad template for everything that is computed from the // variables with respect to which we would like to compute // derivatives. This includes the current solution and gradient at the // quadrature points (which are linear combinations of the degrees of @@ -1679,12 +1685,12 @@ namespace Step33 Table<2,double> W_old (n_q_points, EulerEquations::n_components); - Table<2,Sacado::Fad::DFad > - W_theta (n_q_points, EulerEquations::n_components); - Table<3,Sacado::Fad::DFad > grad_W (n_q_points, EulerEquations::n_components, dim); + Table<3,double> + grad_W_old(n_q_points, EulerEquations::n_components, dim); + std::vector residual_derivatives (dofs_per_cell); // Next, we have to define the independent variables that we will try to @@ -1709,8 +1715,8 @@ namespace Step33 independent_local_dof_values[i].diff (i, dofs_per_cell); // After all these declarations, let us actually compute something. First, - // the values of W, W_old, W_theta, - // and grad_W, which we can compute from the local DoF values + // the values of W, W_old, grad_W + // and grad_W_old, which we can compute from the local DoF values // by using the formula $W(x_q)=\sum_i \mathbf W_i \Phi_i(x_q)$, where // $\mathbf W_i$ is the $i$th entry of the (local part of the) solution // vector, and $\Phi_i(x_q)$ the value of the $i$th vector-valued shape @@ -1729,9 +1735,11 @@ namespace Step33 { W[q][c] = 0; W_old[q][c] = 0; - W_theta[q][c] = 0; for (unsigned int d=0; d FluxMatrix[EulerEquations::n_components][dim]; FluxMatrix *flux = new FluxMatrix[n_q_points]; + typedef double FluxMatrixOld[EulerEquations::n_components][dim]; + FluxMatrixOld *flux_old = new FluxMatrixOld[n_q_points]; + typedef Sacado::Fad::DFad ForcingVector[EulerEquations::n_components]; ForcingVector *forcing = new ForcingVector[n_q_points]; + typedef double ForcingVectorOld[EulerEquations::n_components]; + ForcingVectorOld *forcing_old = new ForcingVectorOld[n_q_points]; + for (unsigned int q=0; q::compute_flux_matrix (W_theta[q], flux[q]); - EulerEquations::compute_forcing_vector (W_theta[q], forcing[q]); + EulerEquations::compute_flux_matrix (W_old[q], flux_old[q]); + EulerEquations::compute_forcing_vector (W_old[q], forcing_old[q]); + EulerEquations::compute_flux_matrix (W[q], flux[q]); + EulerEquations::compute_forcing_vector (W[q], forcing[q]); } // We now have all of the pieces in place, so perform the assembly. We // have an outer loop through the components of the system, and an inner // loop over the quadrature points, where we accumulate contributions to - // the $i$th residual $F_i$. The general formula for this residual is + // the $i$th residual $R_i$. The general formula for this residual is // given in the introduction and at the top of this function. We can, // however, simplify it a bit taking into account that the $i$th // (vector-valued) test function $\mathbf{z}_i$ has in reality only a // single nonzero component (more on this topic can be found in the @ref // vector_valued module). It will be represented by the variable // component_i below. With this, the residual term can be - // re-written as $F_i = \left(\frac{(\mathbf{w}_{n+1} - + // re-written as + // @f{eqnarray*} + // R_i &=& + // \left(\frac{(\mathbf{w}_{n+1} - // \mathbf{w}_n)_{\text{component\_i}}}{\delta - // t},(\mathbf{z}_i)_{\text{component\_i}}\right)_K$ $- - // \sum_{d=1}^{\text{dim}} \left(\mathbf{F} - // (\tilde{\mathbf{w}})_{\text{component\_i},d}, + // t},(\mathbf{z}_i)_{\text{component\_i}}\right)_K \\ + // &-& \sum_{d=1}^{\text{dim}} \left( \theta \mathbf{F} + // ({\mathbf{w^k_{n+1}}})_{\text{component\_i},d} + (1-\theta) + // \mathbf{F} ({\mathbf{w_{n}}})_{\text{component\_i},d} , // \frac{\partial(\mathbf{z}_i)_{\text{component\_i}}} {\partial - // x_d}\right)_K$ $+ \sum_{d=1}^{\text{dim}} h^{\eta} \left(\frac{\partial - // \mathbf{w}_{\text{component\_i}}}{\partial x_d} , \frac{\partial - // (\mathbf{z}_i)_{\text{component\_i}}}{\partial x_d} \right)_K$ - // $-(\mathbf{G}(\tilde{\mathbf{w}} )_{\text{component\_i}}, - // (\mathbf{z}_i)_{\text{component\_i}})_K$, where integrals are + // x_d}\right)_K \\ + // &+& \sum_{d=1}^{\text{dim}} h^{\eta} \left( \theta \frac{\partial + // \mathbf{w^k_{n+1}}_{\text{component\_i}}}{\partial x_d} + (1-\theta) + // \frac{\partial \mathbf{w_n}_{\text{component\_i}}}{\partial x_d} , + // \frac{\partial (\mathbf{z}_i)_{\text{component\_i}}}{\partial x_d} \right)_K\\ + // &-& \left( \theta\mathbf{G}({\mathbf{w}^k_n+1} )_{\text{component\_i}} + + // (1-\theta)\mathbf{G}({\mathbf{w}_n} )_{\text{component\_i}} , + // (\mathbf{z}_i)_{\text{component\_i}} \right)_K , + // @f} + // where integrals are // understood to be evaluated through summation over quadrature points. // // We initially sum all contributions of the residual in the positive @@ -1804,7 +1828,7 @@ namespace Step33 // this residual. for (unsigned int i=0; i F_i = 0; + Sacado::Fad::DFad R_i = 0; const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; @@ -1817,24 +1841,27 @@ namespace Step33 for (unsigned int point=0; pointdiameter(), + R_i += 1.0*std::pow(fe_v.get_cell()->diameter(), parameters.diffusion_power) * - grad_W[point][component_i][d] * + ( parameters.theta * grad_W[point][component_i][d] + + (1.0-parameters.theta) * grad_W_old[point][component_i][d] ) * fe_v.shape_grad_component(i, point, component_i)[d] * fe_v.JxW(point); - F_i -= forcing[point][component_i] * + R_i -= ( parameters.theta * forcing[point][component_i] + + (1.0 - parameters.theta) * forcing_old[point][component_i] ) * fe_v.shape_value_component(i, point, component_i) * fe_v.JxW(point); } @@ -1842,24 +1869,27 @@ namespace Step33 // At the end of the loop, we have to add the sensitivities to the // matrix and subtract the residual from the right hand side. Trilinos // FAD data type gives us access to the derivatives using - // F_i.fastAccessDx(k), so we store the data in a + // R_i.fastAccessDx(k), so we store the data in a // temporary array. This information about the whole row of local dofs // is then added to the Trilinos matrix at once (which supports the // data types we have chosen). for (unsigned int k=0; kfe_v variable now is of type FEFaceValues or // FESubfaceValues: Table<2,Sacado::Fad::DFad > - Wplus (n_q_points, EulerEquations::n_components), - Wminus (n_q_points, EulerEquations::n_components); + Wplus (n_q_points, EulerEquations::n_components), + Wminus (n_q_points, EulerEquations::n_components); + Table<2,double> + Wplus_old(n_q_points, EulerEquations::n_components), + Wminus_old(n_q_points, EulerEquations::n_components); for (unsigned int q=0; q::compute_Wminus (parameters.boundary_conditions[boundary_id].kind, fe_v.normal_vector(q), Wplus[q], boundary_values[q], Wminus[q]); + // Here we assume that boundary type, boundary normal vector and boundary data values + // maintain the same during time advancing. + EulerEquations::compute_Wminus (parameters.boundary_conditions[boundary_id].kind, + fe_v.normal_vector(q), + Wplus_old[q], + boundary_values[q], + Wminus_old[q]); + } } @@ -1988,6 +2028,9 @@ namespace Step33 typedef Sacado::Fad::DFad NormalFlux[EulerEquations::n_components]; NormalFlux *normal_fluxes = new NormalFlux[n_q_points]; + typedef double NormalFluxOld[EulerEquations::n_components]; + NormalFluxOld *normal_fluxes_old = new NormalFluxOld[n_q_points]; + double alpha; switch (parameters.stabilization_kind) @@ -2004,9 +2047,14 @@ namespace Step33 } for (unsigned int q=0; q::numerical_normal_flux(fe_v.normal_vector(q), Wplus[q], Wminus[q], alpha, normal_fluxes[q]); + EulerEquations::numerical_normal_flux(fe_v.normal_vector(q), + Wplus_old[q], Wminus_old[q], alpha, + normal_fluxes_old[q]); + } // Now assemble the face term in exactly the same way as for the cell // contributions in the previous function. The only difference is that if @@ -2017,34 +2065,36 @@ namespace Step33 for (unsigned int i=0; i F_i = 0; + Sacado::Fad::DFad R_i = 0; for (unsigned int point=0; point