coloredlusgs.c

Go to the documentation of this file.

 
#include <omp.h>
#include "coloredlusgs.h"
#include "flux.h"
 
void parallel_pre_lusgs(int neles, int nface, double factor, \
                        double *fnorm_vol, double *dt, double *diag, double *fspr)
{
    int idx;        // Element index
    int jdx;        // Face index
    double lamf;    // Spectral radius at each face
 
    // SMP applied
    #pragma omp parallel for private(jdx, lamf)
    for (idx=0; idx<neles; idx++) {
        // Diagonals of implicit operator
        diag[idx] = 1.0/(dt[idx]*factor);
 
        for (jdx=0; jdx<nface; jdx++) {
            // Diffusive margin of wave speed of face
            lamf = fspr[neles*jdx + idx]*1.01;
 
            // Save spectral radius
            fspr[neles*jdx + idx] = lamf;
 
            // Add portion of lower and upper spectral radius
            diag[idx] += 0.5*lamf*fnorm_vol[neles*jdx + idx];
        }
    }
}
 
 
void ns_parallel_lower_sweep(int n0, int ne, int neles, int nfvars, int nface, int ndims, \
                             int *nei_ele, int *icolor, int *lcolor, double *fnorm_vol, double *vec_fnorm, \
                             double *uptsb, double *rhsb, double *dub, double *diag, double *fspr)
{   
    int _idx, idx, jdx, kdx, neib, curr_level;
    double du[nfvars], dfj[nfvars], df[nfvars], nf[ndims];
    double u[nfvars], f[nfvars];
    
    // Lower sweep via coloring
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, \
                                     du, dfj, df, nf, u, f)
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        // Initialize `df` array
        for (kdx=0; kdx<nfvars; kdx++) {
            df[kdx] = 0.0;
        }
 
        // Set of faces surrounding a cell
        for (jdx=0; jdx<nface; jdx++) {
 
            // Get face normal vector
            for (kdx=0; kdx<ndims; kdx++) {
                nf[kdx] = vec_fnorm[ndims*neles*jdx + neles*kdx + idx];
            }
            
            // Neighbor element index meet at face
            neib = nei_ele[neles*jdx + idx];
 
            // Only for lower level cell
            if (lcolor[neib] < curr_level) {
 
                for (kdx=0; kdx<nfvars; kdx++) {
                    u[kdx] = uptsb[neles*kdx + neib];
                    du[kdx] = u[kdx] + dub[neles*kdx + neib];
                }
 
                ns_flux_container(nfvars, ndims, u, nf, f);
                ns_flux_container(nfvars, ndims, du, nf, dfj);
 
                for (kdx=0; kdx<nfvars; kdx++) {
                    dfj[kdx] -= f[kdx];
                }
 
                for (kdx=0; kdx<nfvars; kdx++) {
                    df[kdx] += (dfj[kdx] - fspr[neles*jdx + idx] \
                                * dub[neles*kdx + neib])*fnorm_vol[neles*jdx + idx];
                }
            }
        }
 
        // Update dub array
        for (kdx=0; kdx<nfvars; kdx++)
            dub[neles*kdx + idx] = (rhsb[neles*kdx + idx] - 0.5*df[kdx])/diag[idx];
    }
}
 
void rans_parallel_lower_sweep(int n0, int ne, int neles, int nvars, int nfvars, int nface, int ndims, \
                               int *nei_ele, int *icolor, int *lcolor, double *fnorm_vol, double *vec_fnorm, \
                               double *uptsb, double *rhsb, double *dub, double *diag, double *fspr, double *dsrc)
{
    int _idx, idx, jdx, kdx, neib, curr_level;
    int dnv = nvars - nfvars;
    double du[nvars], dfj[dnv], df[dnv], nf[ndims];
    double u[nvars], f[dnv];
 
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, \
                                     du, dfj, df, nf, u, f)
    // Lower sweep via coloring
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        // Initialize `df` array
        for (kdx=0; kdx<dnv; kdx++) {
            df[kdx] = 0.0;
        }
        
        // Set of faces surrounding a cell
        for (jdx=0; jdx<nface; jdx++) {
 
            // Get face normal vector
            for (kdx=0; kdx<ndims; kdx++) {
                nf[kdx] = vec_fnorm[ndims*neles*jdx + neles*kdx + idx];
            }
            
            // Neighbor element index meet at face
            neib = nei_ele[neles*jdx + idx];
 
            // Only for lower level cell
            if (lcolor[neib] < curr_level) {
 
                for (kdx=0; kdx<nvars; kdx++) {
                    u[kdx] = uptsb[neles*kdx + neib];
                    du[kdx] = u[kdx];
                }
 
                for (kdx=nfvars; kdx<nvars; kdx++) {
                    du[kdx] += dub[neles*kdx + neib];
                }
                
                rans_flux_container(nfvars, ndims, dnv, u, nf, f);
                rans_flux_container(nfvars, ndims, dnv, du, nf, dfj);
 
                for (kdx=0; kdx<dnv; kdx++) {
                    dfj[kdx] -= f[kdx];
                }
 
                for (kdx=0; kdx<dnv; kdx++) {
                    df[kdx] += (dfj[kdx] - fspr[neles*jdx + idx] \
                                * dub[neles*(kdx+nfvars) + neib]) * fnorm_vol[neles*jdx + idx];
                }
            }
        }
 
        // Update dub array
        for (kdx=0; kdx<dnv; kdx++) {
            dub[neles*(kdx+nfvars) + idx] = (rhsb[neles*(kdx+nfvars)+idx] - \
                                            0.5*df[kdx])/(diag[idx]+dsrc[neles*(kdx+nfvars)+idx]);
        }
    }
}
 
 
void ns_parallel_upper_sweep(int n0, int ne, int neles, int nfvars, int nface, int ndims, \
                             int *nei_ele, int *icolor, int *lcolor, double *fnorm_vol, double *vec_fnorm, \
                             double *uptsb, double *rhsb, double *dub, double *diag, double *fspr)
{   
    int _idx, idx, jdx, kdx, neib, curr_level;
    double du[nfvars], dfj[nfvars], df[nfvars], nf[ndims];
    double u[nfvars], f[nfvars];
 
    // Upper sweep via coloring
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, \
                                     du, dfj, df, nf, u, f)
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        // Initialize `df` array
        for (kdx=0; kdx<nfvars; kdx++) {
            df[kdx] = 0.0;
        }
        
        // Set of faces surrounding a cell
        for (jdx=0; jdx<nface; jdx++) {
            
            // Get face normal vector
            for (kdx=0; kdx<ndims; kdx++) {
                nf[kdx] = vec_fnorm[ndims*neles*jdx + neles*kdx + idx];
            }
            
            // Neighbor element index meet at face
            neib = nei_ele[neles*jdx + idx];
 
            // Only for upper level cell
            if (lcolor[neib] > curr_level) {
 
                for (kdx=0; kdx<nfvars; kdx++) {
                    u[kdx] = uptsb[neles*kdx + neib];
                    du[kdx] = u[kdx] + rhsb[neles*kdx + neib];
                }
                
                ns_flux_container(nfvars, ndims, u, nf, f);
                ns_flux_container(nfvars, ndims, du, nf, dfj);
 
                for (kdx=0; kdx<nfvars; kdx++) {
                    dfj[kdx] -= f[kdx];
                }
                
                for (kdx=0; kdx<nfvars; kdx++) {
                    df[kdx] += (dfj[kdx] - fspr[neles*jdx + idx] \
                                * rhsb[neles*kdx + neib])*fnorm_vol[neles*jdx + idx];
                }
            }
        }
 
        // Update rhsb array
        for (kdx=0; kdx<nfvars; kdx++) {
            rhsb[neles*kdx + idx] = dub[neles*kdx + idx] - 0.5*df[kdx]/diag[idx];
        }
    }
}
 
 
void rans_parallel_upper_sweep(int n0, int ne, int neles, int nvars, int nfvars, int nface, int ndims, \
                               int *nei_ele, int *icolor, int *lcolor, double *fnorm_vol, double *vec_fnorm, \
                               double *uptsb, double *rhsb, double *dub, double *diag, double *fspr, double *dsrc)
{
    int _idx, idx, jdx, kdx, neib, curr_level;
    int dnv = nvars - nfvars;
    double du[nvars], dfj[dnv], df[dnv], nf[ndims];
    double u[nvars], f[dnv];
 
    #pragma omp parallel for private(idx, jdx, kdx, neib, curr_level, \
                                     du, dfj, df, nf, u, f)
    // Upper sweep via coloring
    for (_idx=n0; _idx<ne; _idx++) {
        idx = icolor[_idx];
        curr_level = lcolor[idx];
 
        // Initialize `df` array
        for (kdx=0; kdx<dnv; kdx++) {
            df[kdx] = 0.0;
        }
        
        // Set of faces surrounding a cell
        for (jdx=0; jdx<nface; jdx++) {
 
            // Get face normal vector
            for (kdx=0; kdx<ndims; kdx++) {
                nf[kdx] = vec_fnorm[ndims*neles*jdx + neles*kdx + idx];
            }
 
            // Neighbor element index meet at face
            neib = nei_ele[neles*jdx + idx];
 
            // Only for upper level cell
            if (lcolor[neib] > curr_level) {
 
                for (kdx=0; kdx<nvars; kdx++) {
                    u[kdx] = uptsb[neles*kdx + neib];
                    du[kdx] = u[kdx];
                }
 
                for (kdx=nfvars; kdx<nvars; kdx++) {
                    du[kdx] += rhsb[neles*kdx + neib];
                }
 
                rans_flux_container(nfvars, ndims, dnv, u, nf, f);
                rans_flux_container(nfvars, ndims, dnv, du, nf, dfj);
 
                for (kdx=0; kdx<dnv; kdx++) {
                    dfj[kdx] -= f[kdx];
                }
 
                for (kdx=0; kdx<dnv; kdx++) {
                    df[kdx] += (dfj[kdx] - fspr[neles*jdx+idx] \
                                * rhsb[neles*(kdx+nfvars)+neib])*fnorm_vol[neles*jdx+idx];
                }
            }
        }
 
        // Update rhsb array
        for (kdx=0; kdx<dnv; kdx++) {
            rhsb[neles*(kdx+nfvars)+idx] = dub[neles*(kdx+nfvars)+idx] - \
                                        0.5*df[kdx]/(diag[idx] + dsrc[neles*(kdx+nfvars)+idx]);
        }
    }
}
 
 
void parallel_update(int neles, int nvars, double *uptsb, double *rhsb)
{
    int idx, kdx;   
    #pragma omp parallel for private(kdx)
    // Iterate for all cell
    for (idx=0; idx<neles; idx++) {
        // Update conservative variables
        for (kdx=0; kdx<nvars; kdx++) {
            // Indexing 2D array as 1D
            uptsb[neles*kdx + idx] += rhsb[neles*kdx + idx];
        }
    }
}