coloredblusgs.c

Go to the documentation of this file.

 
#include "coloredblusgs.h"
#include "flux.h"
#include "inverse.h"
#include <stdio.h>
#include <omp.h>
 
#define nsmatdim NFVARS*NFVARS
#define ransmatdim NTURBVARS*NTURBVARS
 
void ns_parallel_pre_blusgs(UCFD_INT neles, UCFD_INT nface, UCFD_FLOAT factor,
                          UCFD_FLOAT *fnorm_vol, UCFD_FLOAT *dt, UCFD_FLOAT *diag, UCFD_FLOAT *fjmat)
{
    UCFD_INT idx, jdx, kdx, row, col;
    UCFD_FLOAT fv, dti;
    UCFD_FLOAT dmat[NFVARS][NFVARS];
 
    #pragma omp parallel for private(jdx, kdx, row, col, fv, dmat, dti)
    for (idx=0; idx<neles; idx++) {
        // Initialize diagonal matrix
        for (row=0; row<NFVARS; row++) {
            for (col=0; col<NFVARS; col++)
                dmat[row][col] = 0.0;
        }
        
        // Computes diagonal matrix based on neighbor cells
        for (jdx=0; jdx<nface; jdx++) {
            fv = fnorm_vol[neles*jdx + idx];
            for (row=0; row<NFVARS; row++) {
                for (col=0; col<NFVARS; col++) {
                    dmat[row][col] \
                        += fjmat[idx+neles*jdx+nface*neles*col+NFVARS*nface*neles*row]*fv;
                }
            }
        }
 
        // Complete implicit operator
        dti = 1.0/(dt[idx]*factor);
        for (kdx=0; kdx<NFVARS; kdx++) {
            dmat[kdx][kdx] += dti;
        }
        
        // LU decomposition for inverse process
        ludcmp(NFVARS, dmat);
 
        // Allocate temporal matrix to diag array
        for (row=0; row<NFVARS; row++) {
            for (col=0; col<NFVARS; col++) {
                diag[idx+neles*col+neles*NFVARS*row] = dmat[row][col];
            }
        }
    }
}
 
 
void rans_parallel_pre_blusgs(UCFD_INT neles, UCFD_INT nface, UCFD_FLOAT factor,
                            UCFD_FLOAT *fnorm_vol, UCFD_FLOAT *uptsb, UCFD_FLOAT *dt,
                            UCFD_FLOAT *tdiag, UCFD_FLOAT *tjmat, UCFD_FLOAT *dsrc)
{
    UCFD_INT idx, jdx, kdx, row, col;
    UCFD_FLOAT fv;
    UCFD_FLOAT tmat[NTURBVARS][NTURBVARS];
    UCFD_FLOAT uf[NVARS], dsrcf[NVARS];
 
    #pragma omp parallel for private(jdx, kdx, row, col, fv, tmat, uf, dsrcf)
    for (idx=0; idx<neles; idx++) {
        // Initialize diagonal matrix
        for (row=0; row<NTURBVARS; row++) {
            for (col=0; col<NTURBVARS; col++)
                tmat[row][col] = 0.0;
        }
        
        for (kdx=0; kdx<NVARS; kdx++) {
            uf[kdx] = uptsb[idx+neles*kdx];
            dsrcf[kdx] = dsrc[idx+neles*kdx];
        }
        
        // Computes diagonal matrix based on neighbor cells
        for (jdx=0; jdx<nface; jdx++) {
            fv = fnorm_vol[neles*jdx + idx];
            for (row=0; row<NTURBVARS; row++) {
                for (col=0; col<NTURBVARS; col++) {
                    tmat[row][col] \
                        += tjmat[idx+neles*jdx+nface*neles*col+NTURBVARS*nface*neles*row]*fv;
                }
            }
        }
 
        // Computes Source term Jacobian
        if (rans_source_jacobian(uf, tmat, dsrcf) == UCFD_STATUS_NOT_SUPPORTED)
            printf("Error::Invalid `NTURBVARS` value\n");
 
        // Complete implicit operator
        for (kdx=0; kdx<NTURBVARS; kdx++) {
            tmat[kdx][kdx] += 1.0/(dt[idx]*factor);
        }
        
        // LU decomposition for inverse process
        ludcmp(NTURBVARS, tmat);
 
        // Allocate temporal matrix to diag array
        for (row=0; row<NTURBVARS; row++) {
            for (col=0; col<NTURBVARS; col++) {
                tdiag[idx+neles*col+neles*NTURBVARS*row] = tmat[row][col];
            }
        }
    }
}
 
 
void ns_parallel_block_sweep(UCFD_INT n0, UCFD_INT ne, UCFD_INT neles, UCFD_INT nface,
                             UCFD_INT *nei_ele, UCFD_INT *icolor, UCFD_INT *lcolor, UCFD_FLOAT *fnorm_vol,
                             UCFD_FLOAT *rhsb, UCFD_FLOAT *dub, UCFD_FLOAT *diag, UCFD_FLOAT *fjmat)
{
    UCFD_INT _idx, idx, jdx, kdx, neib, curr_level;
    UCFD_INT row, col;
    UCFD_FLOAT val, fv;
    UCFD_FLOAT rhs[NFVARS], dmat[NFVARS][NFVARS];
 
    // Lower/Upper sweep via coloring
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, row, col, rhs, dmat, val, fv)
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        // Initialize
        for (kdx=0; kdx<NFVARS; kdx++) {
            rhs[kdx] = rhsb[idx+kdx*neles];
        }
 
        for (row=0; row<NFVARS; row++) {
            for (col=0; col<NFVARS; col++) {
                dmat[row][col] = diag[idx+neles*col+NFVARS*neles*row];
            }
        }
 
        for (jdx=0; jdx<nface; jdx++) {
            neib = nei_ele[idx+neles*jdx];
 
            if (lcolor[neib] != curr_level) {
                fv = fnorm_vol[idx+neles*jdx];
 
                for (row=0; row<NFVARS; row++) {
                    val = 0.0;
                    for (col=0; col<NFVARS; col++) {
                        val += fjmat[idx+neles*jdx+nface*neles*col+NFVARS*nface*neles*row] \
                                * dub[neib+neles*col];
                    }
                    rhs[row] -= val*fv;
                }
            }
        }
 
        lusub(NFVARS, dmat, rhs);
 
        for (kdx=0; kdx<NFVARS; kdx++) {
            dub[idx+neles*kdx] = rhs[kdx];
        }
    }
}
 
 
void rans_parallel_block_sweep(UCFD_INT n0, UCFD_INT ne, UCFD_INT neles, UCFD_INT nface,
                               UCFD_INT *nei_ele, UCFD_INT *icolor, UCFD_INT *lcolor, UCFD_FLOAT *fnorm_vol,
                               UCFD_FLOAT *rhsb, UCFD_FLOAT *dub, UCFD_FLOAT *tdiag, UCFD_FLOAT *tjmat)
{
    UCFD_INT _idx, idx, jdx, kdx, neib, curr_level;
    UCFD_INT row, col;
    UCFD_FLOAT val, fv;
    UCFD_FLOAT rhs[NTURBVARS], tmat[NTURBVARS][NTURBVARS];
 
    // Lower/Upper sweep via coloring
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, row, col, rhs, tmat, val, fv)
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        //Initialize
        for (kdx=0; kdx<NTURBVARS; kdx++) {
            rhs[kdx] = rhsb[idx+(kdx+NFVARS)*neles];
        }
 
        for (row=0; row<NTURBVARS; row++) {
            for (col=0; col<NTURBVARS; col++) {
                tmat[row][col] = tdiag[idx+neles*col+neles*NTURBVARS*row];
            }
        }
 
        for (jdx=0; jdx<nface; jdx++) {
            neib = nei_ele[idx+neles*jdx];
 
            if (lcolor[neib] != curr_level) {
                fv = fnorm_vol[idx+neles*jdx];
                
                for (row=0; row<NTURBVARS; row++) {
                    val = 0.0;
                    for (col=0; col<NTURBVARS; col++) {
                        val += tjmat[idx+neles*jdx+nface*neles*col+NTURBVARS*nface*neles*row] \
                                * dub[neib+neles*(col+NFVARS)];
                    }
                    rhs[row] -= val*fv;
                }
            }
        }
 
        // Compute inverse of diagonal matrix multiplication
        lusub(NTURBVARS, tmat, rhs);
 
        // Update dub array
        for (kdx=0; kdx<NTURBVARS; kdx++) {
            dub[idx+neles*(kdx+NFVARS)] = rhs[kdx];
        }
    }
}
 
 
void blusgs_parallel_ns_update(UCFD_INT neles, UCFD_FLOAT *uptsb, UCFD_FLOAT *dub, UCFD_FLOAT *subres)
{
    UCFD_INT idx, kdx;
 
    #pragma omp parallel for private(kdx)
    for (idx=0; idx<neles; idx++) {
        for (kdx=0; kdx<NFVARS; kdx++) {
            uptsb[idx+neles*kdx] += dub[idx+neles*kdx];
 
            // Initialize dub array
            dub[idx+neles*kdx] = 0.0;
        }
        // Initialize sub-residual array
        subres[idx] = 0.0;
    }
}
 
 
void blusgs_parallel_update(UCFD_INT neles, UCFD_FLOAT *uptsb, UCFD_FLOAT *dub, UCFD_FLOAT *subres)
{
    UCFD_INT idx, kdx;
 
    #pragma omp parallel for private(kdx)
    for (idx=0; idx<neles; idx++) {
        for (kdx=0; kdx<NVARS; kdx++) {
            uptsb[idx+neles*kdx] += dub[idx+neles*kdx];
 
            // Initialize dub array
            dub[idx+neles*kdx] = 0.0;
        }
        // Initialize sub-residual array
        subres[idx] = 0.0;
    }
}