From: Luca Heltai Date: Tue, 9 Apr 2019 21:05:28 +0000 (+0200) Subject: ScratchData with Sacado applied to Residual X-Git-Tag: v9.1.0-rc1~210^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=da5071eb97fb2c287e666f7da846ac0f61d18421;p=dealii.git ScratchData with Sacado applied to Residual --- diff --git a/tests/meshworker/scratch_data_05.cc b/tests/meshworker/scratch_data_05.cc new file mode 100644 index 0000000000..1d72bead03 --- /dev/null +++ b/tests/meshworker/scratch_data_05.cc @@ -0,0 +1,178 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2018 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + */ + +// Compute Laplacian using mesh_loop, sacado, and ScratchData + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include "../tests.h" + +using namespace MeshWorker; + +typedef Sacado::Fad::DFad Sdouble; +typedef Sacado::Fad::DFad SSdouble; + +template +void +test() +{ + Triangulation tria; + FE_Q fe(1); + DoFHandler dh(tria); + + FunctionParser rhs_function("1"); + FunctionParser boundary_function("0"); + + AffineConstraints constraints; + + GridGenerator::hyper_cube(tria); + tria.refine_global(5); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + + DoFTools::make_hanging_node_constraints(dh, constraints); + VectorTools::interpolate_boundary_values(dh, + 0, + boundary_function, + constraints); + + constraints.close(); + + SparsityPattern sparsity; + + { + DynamicSparsityPattern dsp(dh.n_dofs(), dh.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dh, dsp); + sparsity.copy_from(dsp); + } + + SparseMatrix matrix; + matrix.reinit(sparsity); + + Vector solution(dh.n_dofs()); + Vector rhs(dh.n_dofs()); + + QGauss quad(3); + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1, 1, 1>; + + ScratchData scratch(fe, quad, cell_flags); + CopyData copy(fe.dofs_per_cell); + + auto cell = dh.begin_active(); + auto endc = dh.end(); + + typedef decltype(cell) Iterator; + + FEValuesExtractors::Scalar scalar(0); + + auto cell_worker = [&rhs_function, &solution, &scalar](const Iterator &cell, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + const auto &p = s.get_quadrature_points(); + + Sdouble dummy; + std::vector residual(fev.dofs_per_cell, Sdouble(0)); + + s.extract_local_dof_values("solution", solution, dummy); + const auto &u = s.get_values("solution", scalar, dummy); + const auto &grad_u = s.get_gradients("solution", scalar, dummy); + + c = 0; + + c.local_dof_indices[0] = s.get_local_dof_indices(); + + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + for (unsigned int q = 0; q < p.size(); ++q) + { + residual[i] += (fev.shape_grad(i, q) * grad_u[q] - + rhs_function.value(p[q]) * fev.shape_value(i, q)) * + JxW[q]; + } + + // Copy from Sacado to normal vectors + for (unsigned int i = 0; i < fev.dofs_per_cell; ++i) + { + c.vectors[0][i] = -residual[i].val(); + for (unsigned int j = 0; j < fev.dofs_per_cell; ++j) + c.matrices[0](i, j) = residual[i].dx(j); + } + }; + + auto copier = [&constraints, &matrix, &rhs](const CopyData &c) { + constraints.distribute_local_to_global( + c.matrices[0], c.vectors[0], c.local_dof_indices[0], matrix, rhs); + }; + + mesh_loop(cell, endc, cell_worker, copier, scratch, copy, assemble_own_cells); + + SparseDirectUMFPACK inv; + inv.initialize(matrix); + + inv.vmult(solution, rhs); + constraints.distribute(solution); + + deallog << "Linfty norm of solution: " << solution.linfty_norm() << std::endl; +} + + +int +main() +{ + initlog(1); + MultithreadInfo::set_thread_limit(1); // to make output deterministic + + test<2, 2>(); +} diff --git a/tests/meshworker/scratch_data_05.with_trilinos=on.with_muparser=on.with_umfpack=on.output b/tests/meshworker/scratch_data_05.with_trilinos=on.with_muparser=on.with_umfpack=on.output new file mode 100644 index 0000000000..b5ad838ee8 --- /dev/null +++ b/tests/meshworker/scratch_data_05.with_trilinos=on.with_muparser=on.with_umfpack=on.output @@ -0,0 +1,2 @@ + +DEAL::Linfty norm of solution: 0.0737281