//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



#include "parallel/PetscSolverNavierStokes.h"

namespace AMDiS {

  using namespace std;

  
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

    MSG("RUN NAVIER STOKES SOLVER INIT!\n");

    KSPCreate(mpiCommGlobal, &ksp);
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPGMRES);
    KSPSetOptionsPrefix(ksp, "ns_");
    KSPSetFromOptions(ksp);
    KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
  }


  void PetscSolverNavierStokes::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");

    MSG("RUN NAVIER STOKES PRECONDITIONER INIT!\n");

    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
    vector<int> pressureComponent;
    pressureComponent.push_back(2);

    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);

    KSPSetUp(kspInterior);

    KSP *subKsp;
    int nSubKsp;
    PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);

    TEST_EXIT(nSubKsp == 2)("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");

    KSP velocityKsp = subKsp[0];
    KSP schurKsp = subKsp[1];
    PetscFree(subKsp);

    Mat A00, A01, A10, A11;
    PCFieldSplitGetSchurBlocks(pc, &A00, &A01, &A10, &A11);

    MatInfo minfo;
    PetscInt nRow, nCol;
    MatGetSize(A00, &nRow, &nCol);
    MatGetInfo(A00, MAT_GLOBAL_SUM, &minfo);
    MSG("A00: %d x %d  with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));

    MatGetSize(A01, &nRow, &nCol);
    MatGetInfo(A01, MAT_GLOBAL_SUM, &minfo);
    MSG("A01: %d x %d  with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));

    MatGetSize(A10, &nRow, &nCol);
    MatGetInfo(A10, MAT_GLOBAL_SUM, &minfo);
    MSG("A10: %d x %d  with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));

    MatGetSize(A11, &nRow, &nCol);
    MatGetInfo(A11, MAT_GLOBAL_SUM, &minfo);
    MSG("A11: %d x %d  with %d nnz\n", nRow, nCol, static_cast<int>(minfo.nz_used));

    KSPSetType(velocityKsp, KSPPREONLY);
    PC pcSub;
    KSPGetPC(velocityKsp, &pcSub);
    PCSetType(pcSub, PCLU);
    PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);

    KSPView(velocityKsp, PETSC_VIEWER_STDOUT_WORLD);
    KSPView(schurKsp, PETSC_VIEWER_STDOUT_WORLD);

  }

}