<commit_before>/**
 * \file decoupledibpm.cpp
 * \brief Implementation of the class \c DecoupledIBPMSolver.
 */

// PETSc
# include <petscviewerhdf5.h>

// Decoupled solver
# include "decoupledibpm.h"


DecoupledIBPMSolver::DecoupledIBPMSolver(
        const petibm::type::Mesh &inMesh,
        const petibm::type::Boundary &inBC,
        const petibm::type::BodyPack &inBodies,
        const YAML::Node &node)
{
    initialize(inMesh, inBC, inBodies, node);
} // DecoupledIBPMSolver


PetscErrorCode DecoupledIBPMSolver::initialize(
        const petibm::type::Mesh &inMesh,
        const petibm::type::Boundary &inBC,
        const petibm::type::BodyPack &inBodies,
        const YAML::Node &node) 
{
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    
    NavierStokesSolver::initialize(inMesh, inBC, node);
    
    // continue on stageInitialize
    ierr = PetscLogStagePush(stageInitialize); CHKERRQ(ierr);
    
    // make a reference to provided BodyPack
    bodies = inBodies;

    // create force solver
    ierr = petibm::linsolver::createLinSolver("forces", settings, fSolver); CHKERRQ(ierr);

    // create extra operators for decoupled method
    ierr = createOperators(); CHKERRQ(ierr);

    // create extra vectors for decoupled method
    ierr = createVectors(); CHKERRQ(ierr);

    // set coefficient matrix to the force solver
    ierr = fSolver->setMatrix(EBNH); CHKERRQ(ierr);
    
    // register additional events
    ierr = PetscLogStageRegister("rhsForces", &stageRHSForces); CHKERRQ(ierr);
    ierr = PetscLogStageRegister("solveForces", &stageSolveForces); CHKERRQ(ierr);
    ierr = PetscLogStageRegister("integrateForces", &stageIntegrateForces); CHKERRQ(ierr);
    
    ierr = PetscLogStagePop(); CHKERRQ(ierr);

    PetscFunctionReturn(0);
} // initialize


// create extra operators for decoupled method
PetscErrorCode DecoupledIBPMSolver::createOperators()
{
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    
    Mat     R;
    Mat     MHat;
    Mat     BN;
    
    Vec     RDiag;
    Vec     MHatDiag;

    // create diagonal matrix R and get a Vec holding the diagonal
    ierr = petibm::operators::createR(mesh, R); CHKERRQ(ierr);
    ierr = MatCreateVecs(R, nullptr, &RDiag); CHKERRQ(ierr);
    ierr = MatGetDiagonal(R, RDiag); CHKERRQ(ierr);
    
    // create diagonal matrix M and get a Vec holding the diagonal
    ierr = petibm::operators::createMHead(mesh, MHat); CHKERRQ(ierr);
    ierr = MatCreateVecs(MHat, nullptr, &MHatDiag); CHKERRQ(ierr);
    ierr = MatGetDiagonal(MHat, MHatDiag); CHKERRQ(ierr);

    // create operator E
    ierr = petibm::operators::createDelta(mesh, bc, bodies, E); CHKERRQ(ierr);
    ierr = MatDiagonalScale(E, nullptr, RDiag); CHKERRQ(ierr);
    ierr = MatDiagonalScale(E, nullptr, MHatDiag); CHKERRQ(ierr);
    
    // create operator H
    ierr = MatTranspose(E, MAT_INITIAL_MATRIX, &H); CHKERRQ(ierr);
    ierr = MatDiagonalScale(H, RDiag, nullptr); CHKERRQ(ierr);
<commit_msg>fix a bug: forces with correct trend but wrong scaling<commit_after>/**
 * \file decoupledibpm.cpp
 * \brief Implementation of the class \c DecoupledIBPMSolver.
 */

// PETSc
# include <petscviewerhdf5.h>

// Decoupled solver
# include "decoupledibpm.h"


DecoupledIBPMSolver::DecoupledIBPMSolver(
        const petibm::type::Mesh &inMesh,
        const petibm::type::Boundary &inBC,
        const petibm::type::BodyPack &inBodies,
        const YAML::Node &node)
{
    initialize(inMesh, inBC, inBodies, node);
} // DecoupledIBPMSolver


PetscErrorCode DecoupledIBPMSolver::initialize(
        const petibm::type::Mesh &inMesh,
        const petibm::type::Boundary &inBC,
        const petibm::type::BodyPack &inBodies,
        const YAML::Node &node) 
{
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    
    NavierStokesSolver::initialize(inMesh, inBC, node);
    
    // continue on stageInitialize
    ierr = PetscLogStagePush(stageInitialize); CHKERRQ(ierr);
    
    // make a reference to provided BodyPack
    bodies = inBodies;

    // create force solver
    ierr = petibm::linsolver::createLinSolver("forces", settings, fSolver); CHKERRQ(ierr);

    // create extra operators for decoupled method
    ierr = createOperators(); CHKERRQ(ierr);

    // create extra vectors for decoupled method
    ierr = createVectors(); CHKERRQ(ierr);

    // set coefficient matrix to the force solver
    ierr = fSolver->setMatrix(EBNH); CHKERRQ(ierr);
    
    // register additional events
    ierr = PetscLogStageRegister("rhsForces", &stageRHSForces); CHKERRQ(ierr);
    ierr = PetscLogStageRegister("solveForces", &stageSolveForces); CHKERRQ(ierr);
    ierr = PetscLogStageRegister("integrateForces", &stageIntegrateForces); CHKERRQ(ierr);
    
    ierr = PetscLogStagePop(); CHKERRQ(ierr);

    PetscFunctionReturn(0);
} // initialize


// create extra operators for decoupled method
PetscErrorCode DecoupledIBPMSolver::createOperators()
{
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    
    Mat     R;
    Mat     MHat;
    Mat     BN;
    
    Vec     RDiag;
    Vec     MHatDiag;

    // create diagonal matrix R and get a Vec holding the diagonal
    ierr = petibm::operators::createR(mesh, R); CHKERRQ(ierr);
    ierr = MatCreateVecs(R, nullptr, &RDiag); CHKERRQ(ierr);
    ierr = MatGetDiagonal(R, RDiag); CHKERRQ(ierr);
    
    // create diagonal matrix M and get a Vec holding the diagonal
    ierr = petibm::operators::createMHead(mesh, MHat); CHKERRQ(ierr);
    ierr = MatCreateVecs(MHat, nullptr, &MHatDiag); CHKERRQ(ierr);
    ierr = MatGetDiagonal(MHat, MHatDiag); CHKERRQ(ierr);

    // create a Delta operator and its transpose
    ierr = petibm::operators::createDelta(mesh, bc, bodies, E); CHKERRQ(ierr);
    ierr = MatTranspose(E, MAT_INITIAL_MATRIX, &H); CHKERRQ(ierr);
    
    // create operator E
    ierr = MatDiagonalScale(E, nullptr, RDiag); CHKERRQ(ierr);
    ierr = MatDiagonalScale(E, nullptr, MHatDiag); CHKERRQ(ierr);
    
    // create operator H
    // Note: we do nothing here. Because if f is force (unit [m/s^2]), we need to 
    // convert it to pressure when using it in Eulerian momentum equations.
    // That is, (H*R^{-1}*f). While H = R*Delta, This means the real scattering
    // operator used in momentum equation is just a Delta operator. So we just
    // let H be the Delta operator. One thing to remember is that this is based
    // on the assumption that the size of Lagragian element is equal to the size
    // of the Eulerian grid near by.

{'fraction_non_alphanumeric': 0.1091581868640148, 'fraction_numerical': 0.0009250693802035153, 'mean_word_length': 3.4892733564013843, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3096687', 'n_tokens_mistral': 2283, 'n_tokens_neox': 2216, 'n_words': 548}