Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #include "floquetcorrection.hpp"
       5              : 
       6              : #include <limits>
       7              : #include <mfem.hpp>
       8              : #include "fem/bilinearform.hpp"
       9              : #include "fem/fespace.hpp"
      10              : #include "fem/integrator.hpp"
      11              : #include "linalg/iterative.hpp"
      12              : #include "linalg/jacobi.hpp"
      13              : #include "linalg/rap.hpp"
      14              : #include "models/materialoperator.hpp"
      15              : 
      16              : namespace palace
      17              : {
      18              : 
      19              : template <typename VecType>
      20            0 : FloquetCorrSolver<VecType>::FloquetCorrSolver(const MaterialOperator &mat_op,
      21              :                                               FiniteElementSpace &nd_fespace,
      22              :                                               FiniteElementSpace &rt_fespace, double tol,
      23            0 :                                               int max_it, int print)
      24              : {
      25              :   // Create the mass and cross product operators for Floquet correction.
      26              :   {
      27              :     constexpr bool skip_zeros = false;
      28              :     BilinearForm a(rt_fespace);
      29            0 :     a.AddDomainIntegrator<VectorFEMassIntegrator>();
      30            0 :     std::unique_ptr<Operator> m = a.Assemble(skip_zeros);
      31              :     if constexpr (std::is_same<OperType, ComplexOperator>::value)
      32              :     {
      33            0 :       M = std::make_unique<ComplexParOperator>(std::move(m), nullptr, rt_fespace);
      34              :     }
      35              :     else
      36              :     {
      37              :       M = std::make_unique<ParOperator>(std::move(m), rt_fespace);
      38              :     }
      39              :   }
      40              : 
      41              :   {
      42            0 :     MaterialPropertyCoefficient f(mat_op.MaxCeedAttribute());
      43            0 :     f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetFloquetCross(), 1.0);
      44              :     constexpr bool skip_zeros = false;
      45              :     BilinearForm a(nd_fespace, rt_fespace);
      46            0 :     a.AddDomainIntegrator<VectorFEMassIntegrator>(f);
      47            0 :     std::unique_ptr<Operator> m = a.Assemble(skip_zeros);
      48              :     if constexpr (std::is_same<OperType, ComplexOperator>::value)
      49              :     {
      50            0 :       Cross = std::make_unique<ComplexParOperator>(std::move(m), nullptr, nd_fespace,
      51            0 :                                                    rt_fespace, false);
      52              :     }
      53              :     else
      54              :     {
      55              :       Cross = std::make_unique<ParOperator>(std::move(m), nd_fespace, rt_fespace, false);
      56              :     }
      57            0 :   }
      58              : 
      59              :   // Setup the linear solver.
      60            0 :   auto pcg = std::make_unique<CgSolver<OperType>>(rt_fespace.GetComm(), print);
      61            0 :   pcg->SetInitialGuess(0);
      62              :   pcg->SetRelTol(tol);
      63              :   pcg->SetAbsTol(std::numeric_limits<double>::epsilon());
      64              :   pcg->SetMaxIter(max_it);
      65            0 :   auto jac = std::make_unique<JacobiSmoother<OperType>>(rt_fespace.GetComm());
      66            0 :   ksp = std::make_unique<BaseKspSolver<OperType>>(std::move(pcg), std::move(jac));
      67            0 :   ksp->SetOperators(*M, *M);
      68              : 
      69            0 :   rhs.SetSize(rt_fespace.GetTrueVSize());
      70            0 :   rhs.UseDevice(true);
      71            0 : }
      72              : 
      73              : template <typename VecType>
      74            0 : void FloquetCorrSolver<VecType>::Mult(const VecType &x, VecType &y) const
      75              : {
      76            0 :   Cross->Mult(x, rhs);
      77            0 :   ksp->Mult(rhs, y);
      78            0 : }
      79              : 
      80              : template <typename VecType>
      81            0 : void FloquetCorrSolver<VecType>::AddMult(const VecType &x, VecType &y, ScalarType a) const
      82              : {
      83            0 :   this->Mult(x, rhs);
      84            0 :   rhs *= a;
      85              :   y += rhs;
      86            0 : }
      87              : 
      88              : template class FloquetCorrSolver<ComplexVector>;
      89              : 
      90              : }  // namespace palace