Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ComputeDiracThread.h"
      11             : 
      12             : // Moose Includes
      13             : #include "ParallelUniqueId.h"
      14             : #include "DiracKernel.h"
      15             : #include "Problem.h"
      16             : #include "NonlinearSystem.h"
      17             : #include "MooseVariableFE.h"
      18             : #include "Assembly.h"
      19             : #include "ThreadedElementLoop.h"
      20             : 
      21             : #include "libmesh/threads.h"
      22             : 
      23       38401 : ComputeDiracThread::ComputeDiracThread(FEProblemBase & feproblem,
      24             :                                        const std::set<TagID> & tags,
      25       38401 :                                        bool is_jacobian)
      26             :   : ThreadedElementLoop<DistElemRange>(feproblem),
      27       38401 :     _is_jacobian(is_jacobian),
      28       76802 :     _nl(feproblem.currentNonlinearSystem()),
      29       38401 :     _tags(tags),
      30       38401 :     _dirac_kernels(_nl.getDiracKernelWarehouse())
      31             : {
      32       38401 : }
      33             : 
      34             : // Splitting Constructor
      35           0 : ComputeDiracThread::ComputeDiracThread(ComputeDiracThread & x, Threads::split split)
      36             :   : ThreadedElementLoop<DistElemRange>(x, split),
      37           0 :     _is_jacobian(x._is_jacobian),
      38           0 :     _nl(x._nl),
      39           0 :     _tags(x._tags),
      40           0 :     _dirac_kernels(x._dirac_kernels)
      41             : {
      42           0 : }
      43             : 
      44       38401 : ComputeDiracThread::~ComputeDiracThread() {}
      45             : 
      46             : void
      47       38401 : ComputeDiracThread::pre()
      48             : {
      49             :   // Force TID=0 because we run this object _NON THREADED_
      50             :   // Take this out if we ever get Dirac's working with threads!
      51       38401 :   _tid = 0;
      52       38401 : }
      53             : 
      54             : void
      55       33071 : ComputeDiracThread::subdomainChanged()
      56             : {
      57       33071 :   _fe_problem.subdomainSetup(_subdomain, _tid);
      58             : 
      59       33071 :   std::set<MooseVariableFEBase *> needed_moose_vars;
      60       33071 :   _dirac_kernels.updateVariableDependency(needed_moose_vars, _tid);
      61             : 
      62             :   // Update material dependencies
      63       33071 :   std::unordered_set<unsigned int> needed_mat_props;
      64       33071 :   _dirac_kernels.updateMatPropDependency(needed_mat_props, _tid);
      65             : 
      66       33071 :   _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
      67       33071 :   _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
      68             : 
      69             :   // If users pass a empty vector or a full size of vector,
      70             :   // we take all kernels
      71       33071 :   if (!_tags.size() || _tags.size() == _fe_problem.numMatrixTags())
      72       22033 :     _dirac_warehouse = &_dirac_kernels;
      73             :   // If we have one tag only,  We call tag based storage
      74       11038 :   else if (_tags.size() == 1)
      75           0 :     _dirac_warehouse = _is_jacobian
      76           0 :                            ? &(_dirac_kernels.getMatrixTagObjectWarehouse(*(_tags.begin()), _tid))
      77           0 :                            : &(_dirac_kernels.getVectorTagObjectWarehouse(*(_tags.begin()), _tid));
      78             :   // This one may be expensive, and hopefully we do not use it so often
      79             :   else
      80       21984 :     _dirac_warehouse = _is_jacobian ? &(_dirac_kernels.getMatrixTagsObjectWarehouse(_tags, _tid))
      81       10946 :                                     : &(_dirac_kernels.getVectorTagsObjectWarehouse(_tags, _tid));
      82       33071 : }
      83             : 
      84             : void
      85      341949 : ComputeDiracThread::onElement(const Elem * elem)
      86             : {
      87      341949 :   const bool has_dirac_kernels_on_elem = _fe_problem.reinitDirac(elem, _tid);
      88      341949 :   if (!has_dirac_kernels_on_elem)
      89           0 :     return;
      90             : 
      91      341949 :   std::set<MooseVariableFEBase *> needed_moose_vars;
      92      341949 :   const auto & dkernels = _dirac_warehouse->getActiveObjects(_tid);
      93             : 
      94             :   // Only call reinitMaterials() if one or more DiracKernels has
      95             :   // actually called getMaterialProperty().  Loop over all the
      96             :   // DiracKernels and check whether this is the case.
      97      695023 :   for (const auto & dirac_kernel : dkernels)
      98             :   {
      99             :     // If any of the DiracKernels have had getMaterialProperty()
     100             :     // called, we need to reinit Materials.
     101      360413 :     if (dirac_kernel->getMaterialPropertyCalled())
     102             :     {
     103        7339 :       _fe_problem.reinitMaterials(_subdomain, _tid, /*swap_stateful=*/false);
     104        7339 :       break;
     105             :     }
     106             :   }
     107             : 
     108      703390 :   for (const auto & dirac_kernel : dkernels)
     109             :   {
     110      361441 :     if (!dirac_kernel->hasPointsOnElem(elem))
     111       15762 :       continue;
     112      345679 :     else if (!_is_jacobian)
     113             :     {
     114      303444 :       dirac_kernel->computeResidual();
     115      303444 :       continue;
     116             :     }
     117             : 
     118             :     // Get a list of coupled variables from the SubProblem
     119             :     const auto & coupling_entries =
     120       42235 :         dirac_kernel->subProblem().assembly(_tid, _nl.number()).couplingEntries();
     121             : 
     122             :     // Loop over the list of coupled variable pairs
     123      124743 :     for (const auto & it : coupling_entries)
     124             :     {
     125       82508 :       const MooseVariableFEBase * const ivariable = it.first;
     126       82508 :       const MooseVariableFEBase * const jvariable = it.second;
     127             : 
     128             :       // A variant of the check that is in
     129             :       // ComputeFullJacobianThread::computeJacobian().  We
     130             :       // only want to call computeOffDiagJacobian() if both
     131             :       // variables are active on this subdomain, and the
     132             :       // off-diagonal variable actually has dofs.
     133       82508 :       if (dirac_kernel->variable().number() == ivariable->number() &&
     134      124927 :           ivariable->activeOnSubdomain(_subdomain) && jvariable->activeOnSubdomain(_subdomain) &&
     135       42419 :           (jvariable->numberOfDofs() > 0))
     136             :       {
     137       42419 :         dirac_kernel->prepareShapes(jvariable->number());
     138       42419 :         dirac_kernel->computeOffDiagJacobian(jvariable->number());
     139             :       }
     140             :     }
     141             :   }
     142             : 
     143             :   // Note that we do not call swapBackMaterials() here as they were
     144             :   // never swapped in the first place.  This avoids messing up
     145             :   // stored values of stateful material properties.
     146      341949 : }
     147             : 
     148             : void
     149      341949 : ComputeDiracThread::postElement(const Elem * /*elem*/)
     150             : {
     151      341949 :   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     152      341949 :   if (!_is_jacobian)
     153      299901 :     _fe_problem.addResidual(_tid);
     154             :   else
     155       42048 :     _fe_problem.addJacobian(_tid);
     156      341949 : }
     157             : 
     158             : void
     159       38401 : ComputeDiracThread::post()
     160             : {
     161       38401 :   _fe_problem.clearActiveElementalMooseVariables(_tid);
     162       38401 :   _fe_problem.clearActiveMaterialProperties(_tid);
     163       38401 : }
     164             : 
     165             : void
     166           0 : ComputeDiracThread::join(const ComputeDiracThread & /*y*/)
     167             : {
     168           0 : }
     169             : 
     170             : void
     171       38401 : ComputeDiracThread::printGeneralExecutionInformation() const
     172             : {
     173       38401 :   if (!_fe_problem.shouldPrintExecution(_tid))
     174       34193 :     return;
     175        4208 :   const auto & console = _fe_problem.console();
     176        4208 :   console << "[DBG] Executing Dirac Kernels on " << _fe_problem.getCurrentExecuteOnFlag().name()
     177        4208 :           << std::endl;
     178             : }
     179             : 
     180             : void
     181       33071 : ComputeDiracThread::printBlockExecutionInformation() const
     182             : {
     183       39503 :   if (!_fe_problem.shouldPrintExecution(_tid) || _blocks_exec_printed.count(_subdomain) ||
     184        6432 :       !_dirac_warehouse->hasActiveBlockObjects(_subdomain, _tid))
     185       26639 :     return;
     186             : 
     187        6432 :   const auto & dkernels = _dirac_warehouse->getActiveBlockObjects(_subdomain, _tid);
     188        6432 :   const auto & console = _fe_problem.console();
     189        6432 :   console << "[DBG] Ordering of DiracKernels on subdomain " << _subdomain << std::endl;
     190        6432 :   printExecutionOrdering<DiracKernelBase>(dkernels, false);
     191        6432 :   _blocks_exec_printed.insert(_subdomain);
     192             : }