Detailed Description

Definition at line 20 of file ParallelVectorMatrixAssembler.h.

#include <ParallelVectorMatrixAssembler.h>

Collaboration diagram for ProcessLib::Assembly::ParallelVectorMatrixAssembler:

Public Member Functions
	ParallelVectorMatrixAssembler (AbstractJacobianAssembler &jacobian_assembler)

void	assembleWithJacobian (BaseLib::PolymorphicRandomAccessContainerView< LocalAssemblerInterface > const &local_assemblers, std::vector< std::size_t > const &active_elements, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &xs, std::vector< GlobalVector * > const &x_prevs, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)

Private Attributes
AbstractJacobianAssembler &	jacobian_assembler_

LocalMatrixOutput	local_matrix_output_

GlobalMatrixOutput	global_matrix_output_

int const	num_threads_

Constructor & Destructor Documentation

◆ ParallelVectorMatrixAssembler()

ProcessLib::Assembly::ParallelVectorMatrixAssembler::ParallelVectorMatrixAssembler ( AbstractJacobianAssembler & jacobian_assembler )

explicit

Definition at line 110 of file ParallelVectorMatrixAssembler.cpp.

    : jacobian_assembler_{jacobian_assembler},
      num_threads_(getNumberOfThreads())
{
}

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::Assembly::ParallelVectorMatrixAssembler::assembleWithJacobian	(	BaseLib::PolymorphicRandomAccessContainerView< LocalAssemblerInterface > const &	local_assemblers,
		std::vector< std::size_t > const &	active_elements,
		std::vector< NumLib::LocalToGlobalIndexMap const * > const &	dof_tables,
		const double	t,
		double const	dt,
		std::vector< GlobalVector * > const &	xs,
		std::vector< GlobalVector * > const &	x_prevs,
		int const	process_id,
		GlobalMatrix &	M,
		GlobalMatrix &	K,
		GlobalVector &	b,
		GlobalMatrix &	Jac )

Definition at line 117 of file ParallelVectorMatrixAssembler.cpp.

{
    // checks //////////////////////////////////////////////////////////////////
    if (process_id != 0)
    {
        OGS_FATAL("Process id is not 0 but {}", process_id);
    }
 
    if (dof_tables.size() != 1)
    {
        OGS_FATAL("More than 1 dof table");
    }
    auto const& dof_table = *(dof_tables.front());
 
    if (xs.size() != 1)
    {
        OGS_FATAL("More than 1 solution vector");
    }
    auto const& x = *xs.front();
 
    if (x_prevs.size() != 1)
    {
        OGS_FATAL("More than 1 x_prev vector");
    }
    auto const& x_prev = *x_prevs.front();
 
    // algorithm ///////////////////////////////////////////////////////////////
 
    auto stats = CumulativeStats<MultiStats>::create();
    ConcurrentMatrixView M_view(M);
    ConcurrentMatrixView K_view(K);
    ConcurrentMatrixView b_view(b);
    ConcurrentMatrixView Jac_view(Jac);
 
    ThreadException exception;
    bool assembly_error = false;
#pragma omp parallel num_threads(num_threads_)
    {
#ifdef _OPENMP
#pragma omp single nowait
        {
            INFO("Number of threads: {}", omp_get_num_threads());
        }
#endif
 
        // temporary data only stored here in order to avoid frequent memory
        // reallocations.
        std::vector<double> local_M_data;
        std::vector<double> local_K_data;
        std::vector<double> local_b_data;
        std::vector<double> local_Jac_data;
        std::vector<GlobalIndexType> indices;
 
        // copy to avoid concurrent access
        auto const jac_asm = jacobian_assembler_.copy();
        auto stats_this_thread = stats->clone();
 
        MultiMatrixElementCache cache{M_view, K_view, b_view, Jac_view,
                                      stats_this_thread->data};
 
        // TODO corner case: what if all elements on a submesh are deactivated?
        if (active_elements.empty())
        {
            // due to MSVC++ error:
            // error C3016: 'element_id': index variable in OpenMP 'for'
            // statement must have signed integral type
            std::ptrdiff_t const n_loc_asm =
                static_cast<std::ptrdiff_t>(local_assemblers.size());
 
#pragma omp for nowait
            for (std::ptrdiff_t element_id = 0; element_id < n_loc_asm;
                 ++element_id)
            {
                if (assembly_error)
                {
                    continue;
                }
                auto& loc_asm = local_assemblers[element_id];
 
                try
                {
                    assembleWithJacobianOneElement(
                        element_id, loc_asm, dof_table, t, dt, x, x_prev,
                        local_M_data, local_K_data, local_b_data,
                        local_Jac_data, indices, *jac_asm, cache);
                }
                catch (...)
                {
                    exception.capture();
                    assembly_error = true;
                    continue;
                }
 
                local_matrix_output_(t, process_id, element_id, local_M_data,
                                     local_K_data, local_b_data,
                                     &local_Jac_data);
            }
        }
        else
        {
            // due to MSVC++ error:
            // error C3016: 'i': index variable in OpenMP 'for' statement must
            // have signed integral type
            std::ptrdiff_t const n_act_elem =
                static_cast<std::ptrdiff_t>(active_elements.size());
 
#pragma omp for nowait
            for (std::ptrdiff_t i = 0; i < n_act_elem; ++i)
            {
                if (assembly_error)
                {
                    continue;
                }
 
                auto const element_id = active_elements[i];
                auto& loc_asm = local_assemblers[element_id];
 
                try
                {
                    assembleWithJacobianOneElement(
                        element_id, loc_asm, dof_table, t, dt, x, x_prev,
                        local_M_data, local_K_data, local_b_data,
                        local_Jac_data, indices, *jac_asm, cache);
                }
                catch (...)
                {
                    exception.capture();
                    assembly_error = true;
                    continue;
                }
 
                local_matrix_output_(t, process_id, element_id, local_M_data,
                                     local_K_data, local_b_data,
                                     &local_Jac_data);
            }
        }
    }  // OpenMP parallel section
 
    stats->print();
 
    global_matrix_output_(t, process_id, M, K, b, &Jac);
    exception.rethrow();
}