From 81896a114a6dac178312bb616b40b74e7a1b1175 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 20 May 2008 15:23:24 +0000 Subject: [PATCH] Fix a couple of signs. Document assembling cell terms (mostly). git-svn-id: https://svn.dealii.org/trunk@16136 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/doc/intro.dox | 4 +- deal.II/examples/step-33/step-33.cc | 195 ++++++++++++++++++------- 2 files changed, 147 insertions(+), 52 deletions(-) diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index dc3e4a1826..79fdb7a1ee 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -111,8 +111,8 @@ function $\mathbf z$ equals zero: 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}}), -\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) +- \int_{\Omega} \left(\mathbf{F}(\tilde{\mathbf{w}}), +\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) + \int_{\partial \Omega} \left(\mathbf{H}(\tilde{\mathbf{w}}^+), \mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right) diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 2476946222..fb24fe3235 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -1866,14 +1866,35 @@ void ConservationLaw::assemble_system () // 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). + // 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$ 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 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}$. // - // At the top, 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 + // 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 @@ -1920,7 +1941,10 @@ void ConservationLaw::assemble_system () // everything that is computed from them such // as the residual, but not the previous time // step's solution. These variables are all - // found in the first part of the function: + // found in the first part of the function, + // along with a variable that we will use to + // store the derivatives of a single + // component of the residual: template void ConservationLaw:: @@ -1942,6 +1966,8 @@ assemble_cell_term (const FEValues &fe_v, Table<3,Sacado::Fad::DFad > grad_W (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 determine // by solving a Newton step. These @@ -2034,33 +2060,68 @@ assemble_cell_term (const FEValues &fe_v, } - // Gather the flux values for all components at - // all of the quadrature points. This also - // computes the matrix of sensitivities. Perhaps - // this could be done in a better way, since this - // could be a rather large object, but for now it - // seems to work just fine. + // Next, in order to compute the cell + // contributions, we need to evaluate + // $F(\tilde{\mathbf w})$ at all quadrature + // points. To store these, we also need to + // allocate a bit of memory. Note that we + // compute the flux matrices in terms of + // autodifferentiation variables, so that + // the Jacobian contributions can later + // easily be computed from it: typedef Sacado::Fad::DFad FluxMatrix[EulerEquations::n_components][dim]; FluxMatrix *flux = new FluxMatrix[n_q_points]; - for (unsigned int q=0; q < n_q_points; ++q) - EulerEquations::flux_matrix(W_theta[q], flux[q]); + for (unsigned int q=0; q::flux_matrix (W_theta[q], flux[q]); - // We now have all of the function values/grads/fluxes, - // 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 ith residual. + // 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 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} - + // \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}, + // \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$, where integrals are + // understood to be evaluated through + // summation over quadrature points. // - // We initialy sum all contributions of the residual - // in the positive sense, so that we don't need to - // negative the Jacobian entries. Then, when we sum - // into the right_hand_side vector, + // We initialy sum all contributions of the + // residual in the positive sense, so that + // we don't need to negative the Jacobian + // entries. Then, when we sum into the + // right_hand_side vector, // we negate this residual. for (unsigned int i=0; i F_i = 0; + const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; @@ -2068,33 +2129,31 @@ assemble_cell_term (const FEValues &fe_v, // into this fad variable. At the end of the assembly // for this row, we will query for the sensitivities // to this variable and add them into the Jacobian. - Sacado::Fad::DFad F_i; for (unsigned int point=0; pointdiameter(), parameters.diffusion_power) * + + for (unsigned int d=0; ddiameter(), + parameters.diffusion_power) * grad_W[point][component_i][d] * + fe_v.shape_grad_component(i, point, component_i)[d] * fe_v.JxW(point); - // The gravity component only enters into the energy - // equation and into the vertical component of the - // velocity. + // The gravity component only + // enters into the energy equation + // and into the vertical component + // of the velocity. if (component_i == dim - 1) F_i += parameters.gravity * W_theta[point][EulerEquations::density_component] * @@ -2108,14 +2167,51 @@ assemble_cell_term (const FEValues &fe_v, fe_v.JxW(point); } - // Here we gain access to the array of sensitivities - // of the residual. We then sum these into the - // Epetra matrix. - double *values = &(F_i.fastAccessDx(0)); + // 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). The + // code to get Trilinos to add elements + // to the matrix is made a bit more + // awkward by the fact that the + // function takes plain pointers as + // arguments. The first one, taking a + // pointer to + // dofs_per_cell + // double values as its + // third argument is easy enough to + // deal with by just taking the address + // of the first element of the + // residual_derivatives + // variable. However, it also wants an + // int* for the column + // numbers to be written to; this is a + // bit more strenuous because in + // deal.II we always use unsigned + // int to represent indices + // (which are, after all, always + // non-negative), and that the + // dof_indices passed to + // this function are a + // const argument. Why + // Trilinos wants this argument + // non-const is unknown, but in any + // case to make it work we have to + // first cast away the constness, and + // then reinterpret all numbers as + // signed integers. Not pretty but + // works: + for (unsigned int k=0; kSumIntoGlobalValues(dof_indices[i], dofs_per_cell, - values, - reinterpret_cast(const_cast(&dof_indices[0]))); + &residual_derivatives[0], + reinterpret_cast( + const_cast( + &dof_indices[0]))); right_hand_side(dof_indices[i]) -= F_i.val(); } @@ -2142,7 +2238,6 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, const unsigned int boundary_id, const double face_diameter) { - Sacado::Fad::DFad F_i; const unsigned int n_q_points = fe_v.n_quadrature_points; const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell; const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell; @@ -2327,8 +2422,8 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, { if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue; - - F_i = 0; + + Sacado::Fad::DFad F_i = 0; for (unsigned int point=0; point