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 "errorindicator.hpp"
       5              : 
       6              : #include <mfem/general/forall.hpp>
       7              : 
       8              : namespace palace
       9              : {
      10              : 
      11            0 : void ErrorIndicator::AddIndicator(const Vector &indicator)
      12              : {
      13            0 :   if (n == 0)
      14              :   {
      15            0 :     local = indicator;
      16            0 :     n = 1;
      17            0 :     return;
      18              :   }
      19              : 
      20              :   // The average local indicator is used rather than the indicator for the maximum
      21              :   // error to drive the adaptation, to account for a local error that might be marginally
      22              :   // important to many solves, rather than only large in one solve.
      23              :   MFEM_ASSERT(local.Size() == indicator.Size(),
      24              :               "Unexpected size mismatch for ErrorIndicator::AddIndicator!");
      25              : 
      26              :   // The local indicators must be squared before combining, so that the global error
      27              :   // calculation is valid:
      28              :   //                            E = √(1/N ∑ₙ ∑ₖ ηₖₙ²)
      29              :   // from which it follows that:
      30              :   //                            E² = 1/N ∑ₙ ∑ₖ ηₖₙ²
      31              :   //                               = 1/N ∑ₙ Eₙ²
      32              :   // Namely the average of the global error indicators included in the reduction.
      33              :   // Squaring both sides means the summation can be rearranged, and then the local error
      34              :   // indicators become:
      35              :   //                            eₖ = √(1/N ∑ₙ ηₖₙ²)
      36            0 :   const bool use_dev = local.UseDevice() || indicator.UseDevice();
      37              :   const int N = local.Size();
      38            0 :   const int Dn = n;
      39            0 :   const auto *DI = indicator.Read();
      40              :   auto *DL = local.ReadWrite();
      41              :   mfem::forall_switch(
      42            0 :       use_dev, N, [=] MFEM_HOST_DEVICE(int i)
      43            0 :       { DL[i] = std::sqrt((DL[i] * DL[i] * Dn + DI[i] * DI[i]) / (Dn + 1)); });
      44              : 
      45              :   // More samples have been added, update for the running average.
      46            0 :   n += 1;
      47              : }
      48              : 
      49              : }  // namespace palace