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 "ExternalPETScProblem.h"
      11             : #include "SystemBase.h"
      12             : #include "libmesh/petsc_vector.h"
      13             : #include "libmesh/petsc_solver_exception.h"
      14             : 
      15             : registerMooseObject("ExternalPetscSolverApp", ExternalPETScProblem);
      16             : 
      17             : InputParameters
      18         740 : ExternalPETScProblem::validParams()
      19             : {
      20         740 :   InputParameters params = ExternalProblem::validParams();
      21        1480 :   params.addRequiredParam<VariableName>("sync_variable",
      22             :                                         "The variable PETSc external solution will be synced to");
      23         740 :   return params;
      24           0 : }
      25             : 
      26         370 : ExternalPETScProblem::ExternalPETScProblem(const InputParameters & params)
      27             :   : ExternalProblem(params),
      28         370 :     _sync_to_var_name(getParam<VariableName>("sync_variable")),
      29             :     // Require ExternalPetscSolverApp
      30         370 :     _external_petsc_app(static_cast<ExternalPetscSolverApp &>(_app)),
      31         370 :     _ts(_external_petsc_app.getPetscTS()),
      32             :     // RestartableData is required for recovering when PETSc solver runs as a master app
      33         740 :     _petsc_sol(declareRestartableData<Vec>("petsc_sol")),
      34         740 :     _petsc_sol_old(declareRestartableData<Vec>("petsc_sol_old")),
      35        1110 :     _petsc_udot(declareRestartableData<Vec>("petsc_udot"))
      36             : {
      37             :   DM da;
      38         370 :   LibmeshPetscCall(TSGetDM(_ts, &da));
      39             :   // Create a global solution vector
      40         370 :   LibmeshPetscCall(DMCreateGlobalVector(da, &_petsc_sol));
      41             :   // This is required because libMesh incorrectly treats the PETSc parallel vector as a ghost vector
      42             :   // We should be able to remove this line of code once libMesh is updated
      43         370 :   LibmeshPetscCall(VecMPISetGhost(_petsc_sol, 0, nullptr));
      44             :   // The solution at the previous time step
      45         370 :   LibmeshPetscCall(VecDuplicate(_petsc_sol, &_petsc_sol_old));
      46             :   // Udot
      47         370 :   LibmeshPetscCall(VecDuplicate(_petsc_sol, &_petsc_udot));
      48             :   // RHS
      49         370 :   LibmeshPetscCall(VecDuplicate(_petsc_sol, &_petsc_rhs));
      50             :   // Form an initial condition
      51         370 :   LibmeshPetscCall(FormInitialSolution(_ts, _petsc_sol, NULL));
      52         370 :   LibmeshPetscCall(VecCopy(_petsc_sol, _petsc_sol_old));
      53         370 :   LibmeshPetscCall(VecSet(_petsc_udot, 0));
      54         370 : }
      55             : 
      56         740 : ExternalPETScProblem::~ExternalPETScProblem()
      57             : {
      58             :   // Destroy all handles of external Petsc solver
      59         370 :   auto ierr = VecDestroy(&_petsc_sol);
      60         370 :   ierr = VecDestroy(&_petsc_sol_old);
      61         370 :   ierr = VecDestroy(&_petsc_udot);
      62         370 :   ierr = VecDestroy(&_petsc_rhs);
      63             :   // Don't throw during destruction, just abort
      64         370 :   CHKERRABORT(this->comm().get(), ierr);
      65         740 : }
      66             : 
      67             : void
      68        5925 : ExternalPETScProblem::externalSolve()
      69             : {
      70        5925 :   _console << "PETSc External Solve!" << std::endl;
      71             :   // "_petsc_sol_old" is the solution of the current time step, and "_petsc_sol" will be updated
      72             :   // to store the solution of the next time step after this call.
      73             :   // This call advances a time step so that there is an opportunity to
      74             :   // exchange information with MOOSE simulations.
      75        5925 :   LibmeshPetscCall(externalPETScDiffusionFDMSolve(
      76             :       _ts, _petsc_sol_old, _petsc_sol, dt(), time(), &_petsc_converged));
      77        5925 : }
      78             : 
      79             : // This function is called when MOOSE time stepper actually moves to the next time step
      80             : // "PostTimeStep" may not be called for certain cases (e.g., auto_advance=false)
      81             : void
      82        4216 : ExternalPETScProblem::advanceState()
      83             : {
      84        4216 :   FEProblemBase::advanceState();
      85             :   // Compute udot using a backward Euler method
      86             :   // If the external code uses a different method,
      87             :   // udot should be retrieved from the external solver
      88        4216 :   LibmeshPetscCall(VecCopy(_petsc_sol, _petsc_udot));
      89        4216 :   LibmeshPetscCall(VecAXPY(_petsc_udot, -1., _petsc_sol_old));
      90        4216 :   LibmeshPetscCall(VecScale(_petsc_udot, 1. / dt()));
      91             :   // Save current solution because we are moving to the next time step
      92        4216 :   LibmeshPetscCall(VecCopy(_petsc_sol, _petsc_sol_old));
      93        4216 : }
      94             : 
      95             : Real
      96           0 : ExternalPETScProblem::computeResidualL2Norm()
      97             : {
      98           0 :   LibmeshPetscCall(
      99             :       TSComputeIFunction(_ts, time(), _petsc_sol, _petsc_udot, _petsc_rhs, PETSC_FALSE));
     100             :   PetscReal norm;
     101           0 :   LibmeshPetscCall(VecNorm(_petsc_rhs, NORM_2, &norm));
     102             : 
     103           0 :   return norm;
     104             : }
     105             : 
     106             : void
     107       11850 : ExternalPETScProblem::syncSolutions(Direction direction)
     108             : {
     109       11850 :   if (direction == Direction::FROM_EXTERNAL_APP)
     110             :   {
     111        5925 :     _console << "syncSolutions from external petsc App" << std::endl;
     112             :     DM da;
     113             :     // xs: start grid point in x direction on local
     114             :     // ys: start grid point in y direciton on local
     115             :     // xm: number of grid points in x direciton on local
     116             :     // ym: number of grid points in y direction on local
     117             :     // Mx: number of grid points in x direction on all processors
     118             :     PetscInt i, j, xs, ys, xm, ym, Mx;
     119             :     PetscScalar ** _petsc_sol_array;
     120        5925 :     LibmeshPetscCall(TSGetDM(_ts, &da));
     121        5925 :     LibmeshPetscCall(DMDAGetInfo(da,
     122             :                                  PETSC_IGNORE,
     123             :                                  &Mx,
     124             :                                  PETSC_IGNORE,
     125             :                                  PETSC_IGNORE,
     126             :                                  PETSC_IGNORE,
     127             :                                  PETSC_IGNORE,
     128             :                                  PETSC_IGNORE,
     129             :                                  PETSC_IGNORE,
     130             :                                  PETSC_IGNORE,
     131             :                                  PETSC_IGNORE,
     132             :                                  PETSC_IGNORE,
     133             :                                  PETSC_IGNORE,
     134             :                                  PETSC_IGNORE));
     135        5925 :     LibmeshPetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
     136        5925 :     LibmeshPetscCall(DMDAVecGetArray(da, _petsc_sol, &_petsc_sol_array));
     137             : 
     138             :     // Take the solution from PETSc, and sync it to one MOOSE variable
     139             :     // We currently support one variable only but it is straightforward
     140             :     // to have multiple moose variables
     141        5925 :     MeshBase & to_mesh = mesh().getMesh();
     142        5925 :     auto & sync_to_var = getVariable(
     143        5925 :         0, _sync_to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
     144             : 
     145       41532 :     for (j = ys; j < ys + ym; j++)
     146      383240 :       for (i = xs; i < xs + xm; i++)
     147             :       {
     148      347633 :         Node * to_node = to_mesh.node_ptr(i + j * Mx);
     149             :         // For the current example, we need to update only one variable.
     150             :         // This line of code is used to make sure users won't make a mistake in the demo input file.
     151             :         // If multiple variables need to be transfered for some use cases, users should
     152             :         // loop over variables and copy necessary data.
     153      347633 :         if (to_node->n_comp(sync_to_var.sys().number(), sync_to_var.number()) > 1)
     154           0 :           mooseError("Does not support multiple components");
     155             : 
     156      347633 :         dof_id_type dof = to_node->dof_number(sync_to_var.sys().number(), sync_to_var.number(), 0);
     157             :         // Copy the solution to the right location
     158      347633 :         sync_to_var.sys().solution().set(dof, _petsc_sol_array[j][i]);
     159             :       }
     160             : 
     161        5925 :     sync_to_var.sys().solution().close();
     162             : 
     163        5925 :     LibmeshPetscCall(DMDAVecRestoreArray(da, _petsc_sol, &_petsc_sol_array));
     164             : 
     165             :     // Make the solution and the current local solution consistent
     166        5925 :     sync_to_var.sys().update();
     167             :   }
     168        5925 :   else if (direction == Direction::TO_EXTERNAL_APP)
     169             :   {
     170        5925 :     _console << "syncSolutions to external petsc App " << std::endl;
     171             :     // We could the similar thing to sync the solution back to PETSc.
     172             :   }
     173       11850 : }