Calculations.cpp

Go to the documentation of this file.

// ==============================================================================================
// Microflow 3D, http://www.microflow.pwr.edu.pl/
// Created by Roman Szafran on 05.06.19.
// Copyright (c) 2019 Wroclaw University of Science and Technology.
// Distributed under the Apache License, Version 2.0. You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0 or see accompanying file license.txt.
// Redistributions of source code must retain the above copyright and license notice.
// ==============================================================================================

#include "Calculations.h"

double MF::Solver_CPU::Calculations::VelocityLBVolumeSum(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr){

    double VelocitySum = 0.0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 1)  // fluid nodes only
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
            }
    }
    return (VelocitySum);
}

double MF::Solver_CPU::Calculations::VelocityLBMeanForVNodes(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 21 && Thread_it->m_NodeType <= 26) { // velocity nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
    }
    return(VelocitySum/n);
}

double MF::Solver_CPU::Calculations::VelocityFizMeanForVNodes(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 21 && Thread_it->m_NodeType <= 26) { // velocity nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
    }
    return(MF::LBPC::ParametersConversion::ULB_toPhys(VelocitySum/n));
}

double MF::Solver_CPU::Calculations::VelocityLBMeanForPNodes(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr)
{
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 31 && Thread_it->m_NodeType <= 36) { // velocity nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
    }
    return(VelocitySum/n);
}

double MF::Solver_CPU::Calculations::VelocityFizMeanForPNodes(
        const std::shared_ptr<MF::Database::ThreadArray> &ThreadArray_Ptr) {
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 31 && Thread_it->m_NodeType <= 36) { // velocity nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
    }
    return(MF::LBPC::ParametersConversion::ULB_toPhys(VelocitySum/n));
}

double MF::Solver_CPU::Calculations::VelocityLBMeanForThread(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr, std::string &ThreadName) {
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_ThreadName == ThreadName && Thread_it->m_NodeType != 61) { // Excluding BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
    }
    return(VelocitySum/n);
}

double MF::Solver_CPU::Calculations::VelocityFizMeanForInlet(const std::shared_ptr<MF::Database::ThreadArray> &ThreadArray_Ptr) {
    double VelocitySum = 0;
    unsigned int n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_IsInletBalanceSurface && Thread_it->m_NodeType != 61) { // Excluding BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
        else if (Thread_it->m_IsInletBalanceSurface) {
            n = n+(Thread_it->getNumberOfNodes()/2); // Half of number of BB nodes because we need cell mean velocity
        }
    }
    return(MF::LBPC::ParametersConversion::ULB_toPhys(VelocitySum/n));
}

double MF::Solver_CPU::Calculations::VelocityFizMeanForOutlet(
        const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    double VelocitySum = 0;
    uint64_t n = 0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_IsOutletBalanceSurface && Thread_it->m_NodeType != 61) { // Excluding BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                VelocitySum += sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                n++;
            }
        }
        else if (Thread_it->m_IsOutletBalanceSurface) {
            n = n+(Thread_it->getNumberOfNodes()/2); // Half of number of BB nodes because we need cell mean velocity
        }
    }
    return(MF::LBPC::ParametersConversion::ULB_toPhys(VelocitySum/n));
}

double MF::Solver_CPU::Calculations::DeltaMassFlowLB(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr, const std::shared_ptr<MF::Solver_CPU::CaseParameters>& CaseParameters_Ptr)
{
    double sum1 = 0.0, sum2 = 0.0, Velocity;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 21 && Thread_it->m_NodeType <= 26) { // velocity nodes type 20
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                if (CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible)
                    sum1 += (Velocity * CaseParameters_Ptr->ReferenceDensityLB_Rho0_LB);
                else //quasi compressible
                    sum1 += (Velocity * Vector4.rho);
            }
        }
        else if (Thread_it->m_NodeType >= 31 && Thread_it->m_NodeType <= 36) { // pressure nodes type 30
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                if (CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible)
                    sum2 += (Velocity * CaseParameters_Ptr->ReferenceDensityLB_Rho0_LB);
                else //quasi compressible
                    sum2 += (Velocity * Vector4.rho);
            }
        }
    }
    //auto NumberOfVelocity20Nodes = MF::Solver_CPU::Calculations::getVelocityNodeNumber_20(ThreadArray_Ptr);
    //auto NumberOfPressure30Nodes = MF::Solver_CPU::Calculations::getPressureNodeNumber_30(ThreadArray_Ptr);
    return(sum1-sum2)/sum1*100;
}

double MF::Solver_CPU::Calculations::DeltaMassFlow_fiz(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr, const std::shared_ptr<MF::Solver_CPU::CaseParameters>& CaseParameters_Ptr)
{
    double sum1 = 0.0, sum2 = 0.0, Velocity;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 21 && Thread_it->m_NodeType <= 26) { // velocity nodes type 20
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = MF::LBPC::ParametersConversion::ULB_toPhys(sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z));
                sum1 += (Velocity * MF::LBPC::ParametersConversion::RhoLB_toRhoPhys(Vector4.rho));
            }
        }
        else if (Thread_it->m_NodeType >= 31 && Thread_it->m_NodeType <= 36) { // pressure nodes type 30
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = MF::LBPC::ParametersConversion::ULB_toPhys(sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z));
                sum2 += (Velocity * MF::LBPC::ParametersConversion::RhoLB_toRhoPhys(Vector4.rho));
            }
        }
    }
    //auto NumberOfVelocity20Nodes = MF::Solver_CPU::Calculations::getVelocityNodeNumber_20(ThreadArray_Ptr);
    //auto NumberOfPressure30Nodes = MF::Solver_CPU::Calculations::getPressureNodeNumber_30(ThreadArray_Ptr);
    //double unitA = CaseParameters_Ptr->its_DXPhysical_Dx_Phys*CaseParameters_Ptr->its_DXPhysical_Dx_Phys;
    return((sum1-sum2)/sum1*100);
}

double MF::Solver_CPU::Calculations::DeltaMassFlowThread_LB(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr,
                                                               const std::shared_ptr<MF::Solver_CPU::CaseParameters>& CaseParameters_Ptr) {
    double sum1 = 0.0, sum2 = 0.0, Velocity;
    MF::Database::Vec4<double> Vector4;
    //uint64_t NrSum1 = 0 , NrSum2 = 0;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_IsInletBalanceSurface && Thread_it->m_NodeType != 61) { // inlet thread except BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                if (CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {
                    sum1 += (Velocity * CaseParameters_Ptr->ReferenceDensityLB_Rho0_LB);
                    //NrSum1++;
                }
                else {//quasi compressible
                    sum1 += (Velocity * Vector4.rho);
                    //NrSum1++;
                }
            }
        }
        else if (Thread_it->m_IsOutletBalanceSurface && Thread_it->m_NodeType != 61) { // outled thread except BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z);
                if (CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) {
                    sum2 += (Velocity * CaseParameters_Ptr->ReferenceDensityLB_Rho0_LB);
                    //NrSum2++;
                }
                else {//quasi compressible
                    sum2 += (Velocity * Vector4.rho);
                    //NrSum2++;
                }
            }
        }
    }
    return(sum1-sum2)/sum1*100;
}

double MF::Solver_CPU::Calculations::DeltaMassFlowThread_fiz(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr,
                                                               const std::shared_ptr<MF::Solver_CPU::CaseParameters>& CaseParameters_Ptr) {
    double sum1 = 0.0, sum2 = 0.0, Velocity;
    MF::Database::Vec4<double> Vector4;
    //uint64_t NrSum1 = 0 , NrSum2 = 0;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_IsInletBalanceSurface && Thread_it->m_NodeType != 61) { // inlet thread except BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = MF::LBPC::ParametersConversion::ULB_toPhys(sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z));
                sum1 += (Velocity * MF::LBPC::ParametersConversion::RhoLB_toRhoPhys(Vector4.rho));
                //NrSum1++;
            }
        }
        else if (Thread_it->m_IsOutletBalanceSurface && Thread_it->m_NodeType != 61) { // outled thread except BB 61 nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                Velocity = MF::LBPC::ParametersConversion::ULB_toPhys(sqrt(Vector4.x * Vector4.x + Vector4.y * Vector4.y + Vector4.z * Vector4.z));
                sum2 += (Velocity * MF::LBPC::ParametersConversion::RhoLB_toRhoPhys(Vector4.rho));
                //NrSum2++;
            }
        }
    }
    return((sum1-sum2)/sum1*100);
}

double MF::Solver_CPU::Calculations::getULBmaks(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr)
{
    double uLB_maks = 0.0;
    MF::Database::Vec4<double> Vector4;

    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 1) { // only fluid nodes
            for (auto & Node_it : *Thread_it->m_NodeArray_Ptr) {
                Vector4 = (*Thread_it->m_pVRLBFunction)(&Node_it);
                if (fabs(Vector4.x) > fabs(uLB_maks)) {
                    uLB_maks = Vector4.x;
                }
                if (fabs(Vector4.y) > fabs(uLB_maks)) {
                    uLB_maks = Vector4.y;
                }
                if (fabs(Vector4.z) > fabs(uLB_maks)) {
                    uLB_maks = Vector4.z;
                }
            }
        }
    }
    return(uLB_maks);
}

uint64_t MF::Solver_CPU::Calculations::getFluidNodeNumber_1(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 1) { // only fluid nodes
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getSolidNodeNumber_1(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 0) { // only fluid nodes
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getBBNodeNumber_61(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 61) { // only fluid nodes
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getVelocityNodeNumber_20(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 21 && Thread_it->m_NodeType <= 26) { // only nodes of type 20
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getPressureNodeNumber_30(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 31 && Thread_it->m_NodeType <= 36) { // only nodes of type 30
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getVelocity0NodeNumber_40(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType >= 41 && Thread_it->m_NodeType <= 46) { // only nodes of type 40
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getPeriodicNodeNumber_4(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (Thread_it->m_NodeType == 4 || Thread_it->m_NodeType == 5) { // only nodes of type 4 and 5
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    }
    return sum;
}

uint64_t MF::Solver_CPU::Calculations::getTotalNodeNumber(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    unsigned int sum = 0;
    for (auto & Thread_it : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
            sum += Thread_it->m_NodeArray_Ptr->size();
        }
    return sum;
}