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       35512 : ComputeDiracThread::ComputeDiracThread(FEProblemBase & feproblem,
      24             :                                        const std::set<TagID> & tags,
      25       35512 :                                        bool is_jacobian)
      26             :   : ThreadedElementLoop<DistElemRange>(feproblem),
      27       35512 :     _is_jacobian(is_jacobian),
      28       71024 :     _nl(feproblem.currentNonlinearSystem()),
      29       35512 :     _tags(tags),
      30       35512 :     _dirac_kernels(_nl.getDiracKernelWarehouse())
      31             : {
      32       35512 : }
      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       35512 : ComputeDiracThread::~ComputeDiracThread() {}
      45             : 
      46             : void
      47       35512 : 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       35512 :   _tid = 0;
      52       35512 : }
      53             : 
      54             : void
      55       29282 : ComputeDiracThread::subdomainChanged()
      56             : {
      57       29282 :   _fe_problem.subdomainSetup(_subdomain, _tid);
      58             : 
      59       29282 :   std::set<MooseVariableFEBase *> needed_moose_vars;
      60       29282 :   _dirac_kernels.updateVariableDependency(needed_moose_vars, _tid);
      61             : 
      62             :   // Update material dependencies
      63       29282 :   std::unordered_set<unsigned int> needed_mat_props;
      64       29282 :   _dirac_kernels.updateMatPropDependency(needed_mat_props, _tid);
      65             : 
      66       29282 :   _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
      67       29282 :   _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       29282 :   if (!_tags.size() || _tags.size() == _fe_problem.numMatrixTags())
      72       19104 :     _dirac_warehouse = &_dirac_kernels;
      73             :   // If we have one tag only,  We call tag based storage
      74       10178 :   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       20274 :     _dirac_warehouse = _is_jacobian ? &(_dirac_kernels.getMatrixTagsObjectWarehouse(_tags, _tid))
      81       10096 :                                     : &(_dirac_kernels.getVectorTagsObjectWarehouse(_tags, _tid));
      82       29282 : }
      83             : 
      84             : void
      85      304022 : ComputeDiracThread::onElement(const Elem * elem)
      86             : {
      87      304022 :   const bool has_dirac_kernels_on_elem = _fe_problem.reinitDirac(elem, _tid);
      88      304022 :   if (!has_dirac_kernels_on_elem)
      89           0 :     return;
      90             : 
      91      304022 :   std::set<MooseVariableFEBase *> needed_moose_vars;
      92      304022 :   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      616416 :   for (const auto & dirac_kernel : dkernels)
      98             :   {
      99             :     // If any of the DiracKernels have had getMaterialProperty()
     100             :     // called, we need to reinit Materials.
     101      318936 :     if (dirac_kernel->getMaterialPropertyCalled())
     102             :     {
     103        6542 :       _fe_problem.reinitMaterials(_subdomain, _tid, /*swap_stateful=*/false);
     104        6542 :       break;
     105             :     }
     106             :   }
     107             : 
     108      623890 :   for (const auto & dirac_kernel : dkernels)
     109             :   {
     110      319868 :     if (!dirac_kernel->hasPointsOnElem(elem))
     111       12856 :       continue;
     112      307012 :     else if (!_is_jacobian)
     113             :     {
     114      269435 :       dirac_kernel->computeResidual();
     115      269435 :       continue;
     116             :     }
     117             : 
     118             :     // Get a list of coupled variables from the SubProblem
     119             :     const auto & coupling_entries =
     120       37577 :         dirac_kernel->subProblem().assembly(_tid, _nl.number()).couplingEntries();
     121             : 
     122             :     // Loop over the list of coupled variable pairs
     123      111049 :     for (const auto & it : coupling_entries)
     124             :     {
     125       73472 :       const MooseVariableFEBase * const ivariable = it.first;
     126       73472 :       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       73472 :       if (dirac_kernel->variable().number() == ivariable->number() &&
     134      111256 :           ivariable->activeOnSubdomain(_subdomain) && jvariable->activeOnSubdomain(_subdomain) &&
     135       37784 :           (jvariable->numberOfDofs() > 0))
     136             :       {
     137       37784 :         dirac_kernel->prepareShapes(jvariable->number());
     138       37784 :         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      304022 : }
     147             : 
     148             : void
     149      304022 : ComputeDiracThread::postElement(const Elem * /*elem*/)
     150             : {
     151      304022 :   if (!_is_jacobian)
     152      266607 :     _fe_problem.addResidual(_tid);
     153             :   else
     154       37415 :     _fe_problem.addJacobian(_tid);
     155      304022 : }
     156             : 
     157             : void
     158       35512 : ComputeDiracThread::post()
     159             : {
     160       35512 :   _fe_problem.clearActiveElementalMooseVariables(_tid);
     161       35512 :   _fe_problem.clearActiveMaterialProperties(_tid);
     162       35512 : }
     163             : 
     164             : void
     165           0 : ComputeDiracThread::join(const ComputeDiracThread & /*y*/)
     166             : {
     167           0 : }
     168             : 
     169             : void
     170       35512 : ComputeDiracThread::printGeneralExecutionInformation() const
     171             : {
     172       35512 :   if (!_fe_problem.shouldPrintExecution(_tid))
     173       32071 :     return;
     174        3441 :   const auto & console = _fe_problem.console();
     175        3441 :   console << "[DBG] Executing Dirac Kernels on " << _fe_problem.getCurrentExecuteOnFlag().name()
     176        3441 :           << std::endl;
     177             : }
     178             : 
     179             : void
     180       29282 : ComputeDiracThread::printBlockExecutionInformation() const
     181             : {
     182       34186 :   if (!_fe_problem.shouldPrintExecution(_tid) || _blocks_exec_printed.count(_subdomain) ||
     183        4904 :       !_dirac_warehouse->hasActiveBlockObjects(_subdomain, _tid))
     184       24378 :     return;
     185             : 
     186        4904 :   const auto & dkernels = _dirac_warehouse->getActiveBlockObjects(_subdomain, _tid);
     187        4904 :   const auto & console = _fe_problem.console();
     188        4904 :   console << "[DBG] Ordering of DiracKernels on subdomain " << _subdomain << std::endl;
     189        9808 :   printExecutionOrdering<DiracKernelBase>(dkernels, false);
     190        4904 :   _blocks_exec_printed.insert(_subdomain);
     191             : }