PABLO_example_00009.cpp

Parallel 3D adaptive mesh refinement (AMR) with data using PABLO.

Parallel 3D adaptive mesh refinement (AMR) with data using PABLOThe load-balance example with data and "compact families" in example 00008 is here rerun for a three-dimensional octree mesh.

To run: ./PABLO_example_00009
To see the result visit: PABLO website

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit is distributed in the hope that it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#if BITPIT_ENABLE_MPI==1
#include <mpi.h>
#endif
 
#include "bitpit_PABLO.hpp"
 
#if BITPIT_ENABLE_MPI==1
#include "PABLO_userDataComm.hpp"
#include "PABLO_userDataLB.hpp"
#endif
 
using namespace std;
using namespace bitpit;
 
// =================================================================================== //
// =================================================================================== //

void run()
{
    int iter = 0;

    PabloUniform pablo9(0,0,0,1,3);

    int idx = 0;
    pablo9.setBalance(idx,false);

    for (iter=1; iter<6; iter++){
        pablo9.adaptGlobalRefine();
    }
 
#if BITPIT_ENABLE_MPI==1
    pablo9.loadBalance();
#endif

    double xc, yc;
    xc = yc = 0.5;
    double radius = 0.25;

    uint32_t nocts = pablo9.getNumOctants();
    uint32_t nghosts = pablo9.getNumGhosts();
    vector<double> oct_data(nocts, 0.0), ghost_data(nghosts, 0.0);

    for (unsigned int i=0; i<nocts; i++){
        vector<array<double,3> > nodes = pablo9.getNodes(i);
        array<double,3> center = pablo9.getCenter(i);
        for (int j=0; j<8; j++){
            double x = nodes[j][0];
            double y = nodes[j][1];
            if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                oct_data[i] = (pow((center[0]-xc),2.0)+pow((center[1]-yc),2.0));
            }
        }
    }

    iter = 0;
    pablo9.updateConnectivity();
    pablo9.writeTest("pablo00009_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);

    int start = 1;
    for (iter=start; iter<start+2; iter++){
        for (unsigned int i=0; i<nocts; i++){
            vector<array<double,3> > nodes = pablo9.getNodes(i);
            array<double,3> center = pablo9.getCenter(i);
            for (int j=0; j<8; j++){
                double x = nodes[j][0];
                double y = nodes[j][1];
                if ((pow((x-xc),2.0)+pow((y-yc),2.0) <= pow(radius,2.0))){
                    if (center[0]<=xc){

                        pablo9.setMarker(i,1);
                    }
                    else{

                        pablo9.setMarker(i,-1);
                    }
                }
            }
        }

        vector<double> oct_data_new;
        vector<uint32_t> mapper;
        vector<bool> isghost;
        pablo9.adapt(true);
        nocts = pablo9.getNumOctants();
        oct_data_new.resize(nocts, 0.0);

        for (uint32_t i=0; i<nocts; i++){
            pablo9.getMapping(i, mapper, isghost);
            if (pablo9.getIsNewC(i)){
                for (int j=0; j<8; j++){
                    if (isghost[j]){
                        oct_data_new[i] += ghost_data[mapper[j]]/8;
                    }
                    else{
                        oct_data_new[i] += oct_data[mapper[j]]/8;
                    }
                }
            }
            else if (pablo9.getIsNewR(i)){
                oct_data_new[i] += oct_data[mapper[0]];
            }
            else{
                oct_data_new[i] += oct_data[mapper[0]];
            }
        }
 
        oct_data = oct_data_new;
        vector<double>().swap(oct_data_new);

        pablo9.updateConnectivity();
        pablo9.writeTest("pablo00009_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);
 
    }
 
#if BITPIT_ENABLE_MPI==1
    uint8_t levels = 4;
    UserDataLB<vector<double> > data_lb(oct_data,ghost_data);
    pablo9.loadBalance(data_lb, levels);
#endif

    pablo9.updateConnectivity();
    pablo9.writeTest("pablo00009_iter"+to_string(static_cast<unsigned long long>(iter)), oct_data);
}

int main(int argc, char *argv[])
{
#if BITPIT_ENABLE_MPI==1
    MPI_Init(&argc,&argv);
#else
    BITPIT_UNUSED(argc);
    BITPIT_UNUSED(argv);
#endif
 
    int nProcs;
    int rank;
#if BITPIT_ENABLE_MPI==1
    MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#else
    nProcs = 1;
    rank   = 0;
#endif
 
    // Initialize the logger
    log::manager().initialize(log::MODE_SEPARATE, false, nProcs, rank);
    log::cout() << log::fileVerbosity(log::LEVEL_INFO);
    log::cout() << log::disableConsole();
 
    // Run the example
    try {
        run();
    } catch (const std::exception &exception) {
        log::cout() << exception.what();
        exit(1);
    }
 
#if BITPIT_ENABLE_MPI==1
    MPI_Finalize();
#endif
}