MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
diagonal_preconditioner.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
13#include "general/forall.hpp"
14#include "linalg/tensor.hpp"
15
16using mfem::internal::tensor;
17
18namespace mfem
19{
20
22{
23 gradient_operator_ = dynamic_cast<const ElasticityGradientOperator *>(&op);
24 MFEM_ASSERT(gradient_operator_ != nullptr,
25 "Operator is not ElasticityGradientOperator");
26
27 width = height = op.Height();
28
29 gradient_operator_->AssembleGradientDiagonal(Ke_diag_, K_diag_local_,
30 K_diag_);
31
32 submat_height_ = gradient_operator_->elasticity_op_.h1_fes_.GetVDim();
33 num_submats_ = gradient_operator_->elasticity_op_.h1_fes_.GetTrueVSize() /
34 gradient_operator_->elasticity_op_.h1_fes_.GetVDim();
35}
36
38{
39 const int ns = num_submats_, sh = submat_height_, nsh = ns * sh;
40
41 const auto K_diag_submats = Reshape(K_diag_.Read(), ns, sh, sh);
42 const auto X = Reshape(x.Read(), ns, dim);
43
44 auto Y = Reshape(y.Write(), ns, dim);
45
46 if (type_ == Type::Diagonal)
47 {
48 // Assuming Y and X are ordered byNODES. K_diag is ordered byVDIM.
49 mfem::forall(nsh, [=] MFEM_HOST_DEVICE (int si)
50 {
51 const int s = si / sh;
52 const int i = si % sh;
53 Y(s, i) = X(s, i) / K_diag_submats(s, i, i);
54 });
55 }
56 else if (type_ == Type::BlockDiagonal)
57 {
58 mfem::forall(ns, [=] MFEM_HOST_DEVICE (int s)
59 {
60 const auto submat = make_tensor<dim, dim>(
61 [&](int i, int j) { return K_diag_submats(s, i, j); });
62
63 const auto submat_inv = inv(submat);
64
65 const auto x_block = make_tensor<dim>([&](int i) { return X(s, i); });
66
67 tensor<real_t, dim> y_block = submat_inv * x_block;
68
69 for (int i = 0; i < dim; i++)
70 {
71 Y(s, i) = y_block(i);
72 }
73 });
74 }
75 else
76 {
77 MFEM_ABORT("Unknown ElasticityDiagonalPreconditioner::Type");
78 }
79}
80
81} // namespace mfem
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
ElasticityGradientOperator is a wrapper class to pass ElasticityOperator::AssembleGradientDiagonal an...
void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local, Vector &K_diag) const
ParFiniteElementSpace h1_fes_
H1 finite element space.
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
Vector data type.
Definition vector.hpp:80
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:482
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
void forall(int N, lambda &&body)
Definition forall.hpp:754
RefCoord s[3]
Implementation of the tensor class.