"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/Kernel/Operation_HPDDM.cpp" (29 Dec 2018, 3111 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Operation_HPDDM.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 3.0.4_vs_3.1.0.

    1 // GetDP - Copyright (C) 1997-2019 P. Dular and C. Geuzaine, University of Liege
    2 //
    3 // See the LICENSE.txt file for license information. Please report all
    4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
    5 //
    6 // Contributed by Bertrand Thierry
    7 
    8 #include <map>
    9 #include <string>
   10 #include <vector>
   11 #include "ProData.h"
   12 #include "DofData.h"
   13 #include "Message.h"
   14 #include "GetDPConfig.h"
   15 
   16 #if defined(HAVE_HPDDM) && defined(HAVE_PETSC)
   17 
   18 #define MUMPSSUB // direct mumps solver for subdomain solves
   19 #define DMUMPS // direct mumps solver for the coarse solves
   20 #define HPDDM_NUMBERING 'C' // 0-based numering
   21 #define HPDDM_NO_REGEX
   22 #undef NONE
   23 #include <HPDDM.hpp>
   24 
   25 typedef std::complex<double> K;
   26 //typedef double K;
   27 
   28 void Operation_HPDDMSolve(struct Operation *Operation_P,
   29                           struct DofData *DofData_P)
   30 {
   31   // in: idom (will always be == MPI_COMM_RANK currently, but it would be nice
   32   //     to generalize)
   33   // in: D~{idom}() list of neighbors -> o
   34   // in: A~{idom} system for idom (i.e. access to 1 formualation & dofdata -> local dof keys)
   35   //     as well as RHS b~{idom}: thsea are both in "unassembled form", i.e. the local
   36   //     subdomain contributions, with homogeneous Neumann BCs
   37   // in: Optional: B~{idom} same as A~{idom} but with optimized BCs
   38   // in: mapping of indices from common Dofs between A~{idom} and neighbors A~{D~{idom}()}
   39 
   40   // TODO BT: get info from getdp list D~{idom}()
   41   std::vector<int> o(1);
   42   if(Message::GetCommRank() == 1)
   43     o.push_back(2);
   44   else
   45     o.push_back(1);
   46 
   47   // TODO CG: contruct mapping and store somewhere
   48   std::vector<std::vector<int> > mapping(o.size());
   49 
   50   // For PETSc ; more interfaces will follow :-)
   51   Mat PetscA = DofData_P->A.M;
   52   Vec PetscB = DofData_P->b.V;
   53   PetscInt n;
   54   const PetscInt* ia;
   55   const PetscInt* ja;
   56   PetscScalar* array;
   57   PetscBool done;
   58   MatGetRowIJ(PetscA, 0, PETSC_FALSE, PETSC_FALSE, &n, &ia, &ja, &done);
   59   K *b = new K[n];
   60   K *sol = new K[n];
   61   int *ic = new int[n + 1];
   62   MatSeqAIJGetArray(PetscA, &array);
   63   std::copy_n(ia, n + 1, ic);
   64   K *c = new K[ia[n]];
   65   int *jc = new int[ia[n]];
   66   std::copy_n(array, ia[n], c);
   67   std::copy_n(ja, ia[n], jc);
   68   MatSeqAIJRestoreArray(PetscA, &array);
   69   MatRestoreRowIJ(PetscA, 0, PETSC_FALSE, PETSC_FALSE, &n, &ia, &ja, &done);
   70   VecGetArray(PetscB, &array);
   71   std::copy_n(array, n, b);
   72   VecRestoreArray(PetscB, &array);
   73   // last arg is for HPDDM to own the pointers -> should be false, which would
   74   // allow to not copy PETSc data
   75   HPDDM::MatrixCSR<K> *A = new HPDDM::MatrixCSR<K>(n, n, ic[n], c, ic, jc, false, true);
   76 
   77   // Schwarz
   78   HpSchwarz<K> *S = new HpSchwarz<K>();
   79   S->Subdomain::initialize(A, o, mapping); // add communication (last arg)
   80   S->callNumfact(); // if B~{idom} provided, add it here
   81   int mu = 1; // #rhs
   82   int it = HPDDM::IterativeMethod::solve(*S, b, sol, mu, S->getCommunicator());
   83 
   84   // TODO: FETI et BDDM
   85 }
   86 
   87 #else
   88 
   89 void Operation_HPDDMSolve(struct Operation *Operation_P,
   90                           struct DofData *DofData_P)
   91 {
   92   Message::Error("This version of GetDP is not compiled with HPDDM support");
   93 }
   94 
   95 #endif