ERF_InitForEnsemble.cpp File Reference

#include <ERF.H>
#include <ERF_TileNoZ.H>
#include <AMReX_PlotFileUtil.H>
#include <filesystem>

Include dependency graph for ERF_InitForEnsemble.cpp:

Functions
void	ApplyNeumannBCs (const Geometry &geom, MultiFab &mf_cc)
void	NormalizeMultiFabRMS_PerComponent (MultiFab &mf_cc_pert)
void	ReadCustomDataFile (const std::string &filename_custom, int &nx, int &ny, int &nz, int &ng, int &ncomp, std::array< Real, 3 > &problo_ext, std::array< Real, 3 > &probhi_ext, Vector< Real > &data_rho, Vector< Real > &data_theta, Vector< Real > &data_xvel, Vector< Real > &data_yvel, Vector< Real > &data_zvel, Vector< Real > &data_qv, Vector< Real > &data_qc, Vector< Real > &data_qrain)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int	idx (int i, int j, int k, int nx, int ny)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	interp_trilinear (const Real *f, int i, int j, int k, Real tx, Real ty, Real tz, int nx, int ny, int nz)
void	InterpolateToFineMF (const Vector< Real > &data_rho, const Vector< Real > &data_theta, const Vector< Real > &data_xvel, const Vector< Real > &data_yvel, const Vector< Real > &data_zvel, const Vector< Real > &data_qv, const Vector< Real > &data_qc, const Vector< Real > &data_qrain, int nx, int ny, int nz, const std::array< Real, 3 > &problo, const std::array< Real, 3 > &probhi, MultiFab &mf_fine, const Geometry &geom_fine)
void	MakeFinalMultiFabs (const MultiFab &mf_cc_fine, MultiFab &cons_pert, MultiFab &xvel_pert, MultiFab &yvel_pert, MultiFab &zvel_pert, const int n_qstate_moist)
void	AddPertToBckgnd (MultiFab &mf_cc_fine, const MultiFab &mf_cc_pert, const Real &ens_pert_amplitude)

Variables
Vector< Real >	perturb_scale

Function Documentation

AddPertToBckgnd()

void AddPertToBckgnd	(	MultiFab &	mf_cc_fine,
		const MultiFab &	mf_cc_pert,
		const Real &	ens_pert_amplitude
	)

Add normalized perturbations to the interpolated background state.

Parameters

mf_cc_fine	Background state modified in place
mf_cc_pert	Perturbation field to apply
ens_pert_amplitude	Relative perturbation amplitude

{
    const int ncomp = mf_cc_fine.nComp();
 
    // Optional safety check (recommended)
    AMREX_ALWAYS_ASSERT(mf_cc_pert.nComp() == ncomp);
 
    // Create a GPU-accessible array and copy the host vector into it
    GpuArray<Real, 8> scale_gpu; // Use the maximum possible ncomp, e.g., 8
    for (int n = 0; n < ncomp; ++n) {
        scale_gpu[n] = perturb_scale[n];
    }
 
    for (MFIter mfi(mf_cc_fine, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
 
        auto const& bg   = mf_cc_fine.array(mfi);
        auto const& pert = mf_cc_pert.const_array(mfi);
 
        amrex::ParallelFor(bx, ncomp,
        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
        {
            Real ens_amp = ens_pert_amplitude*scale_gpu[n]*std::abs(bg(i,j,k,n));
            bg(i,j,k,n) += ens_amp*pert(i,j,k,n);
        });
    }
}

Referenced by ERF::create_background_state_for_ensemble().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyNeumannBCs()

void ApplyNeumannBCs	(	const Geometry &	geom,
		MultiFab &	mf_cc
	)

Fill periodic ghost cells and apply first-order extrapolation at boundaries.

Parameters

geom	Geometry defining domain bounds and periodicity
mf_cc	Cell-centered MultiFab whose ghost cells are filled in place

{
 
     // -------------------------------------------------
    // 2. Fill interior + periodic ghost cells
    // -------------------------------------------------
    mf_cc.FillBoundary(geom.periodicity());
    // -------------------------------------------------
    // 3. Apply FOExtrap (Neumann) at domain boundaries
    // -------------------------------------------------
    const Box& domain = geom.Domain();
 
    const bool periodic_x = geom.isPeriodic(0);
    const bool periodic_y = geom.isPeriodic(1);
    const bool periodic_z = geom.isPeriodic(2);
 
    for (MFIter mfi(mf_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& gbx = mfi.growntilebox();   // includes ghost cells
        const Box& vbx = mfi.validbox();
 
        auto const& arr = mf_cc.array(mfi);
        int ncomp = mf_cc.nComp();
 
        ParallelFor(gbx, ncomp,
        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
        {
            if (vbx.contains(i,j,k)) return;
 
            int ii = i;
            int jj = j;
            int kk = k;
 
             if (!periodic_x) {
                ii = amrex::max(domain.smallEnd(0),
                                amrex::min(i, domain.bigEnd(0)));
            }
 
            if (!periodic_y) {
                jj = amrex::max(domain.smallEnd(1),
                                amrex::min(j, domain.bigEnd(1)));
            }
 
            if (!periodic_z) {
                kk = amrex::max(domain.smallEnd(2),
                                amrex::min(k, domain.bigEnd(2)));
            }
 
            arr(i,j,k,n) = arr(ii,jj,kk,n);
        });
    }
}

Referenced by ERF::apply_gaussian_smoothing_to_perturbations(), and ERF::create_background_state_for_ensemble().

Here is the call graph for this function:

Here is the caller graph for this function:

idx()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int idx	(	int	i,
		int	j,
		int	k,
		int	nx,
		int	ny
	)

Convert a three-dimensional cell index into a flattened array index.

Parameters

i	x-index
j	y-index
k	z-index
nx	Number of cells in the x-direction
ny	Number of cells in the y-direction

Returns

Flattened index into a row-major 3D array

{
    return i + nx * (j + ny * k);
}

Referenced by NOAHMP::Advance_With_State(), br_shift(), IBFaceSet::build(), ChopGrids2D(), clamp_interp_index(), closest_index(), Compute_InvDt_Compressible(), SurfaceLayer::compute_sfc_params_from_lsm_fluxes(), ERF::estTimeStep(), ERF::FillPlot3DVars(), ERF::FillSurfaceStateMultiFabs(), NOAHMP::Init(), ERF::init_stuff(), ERF::InitData_post(), MRISplitIntegrator< T >::initialize_data(), interp_trilinear(), NOAHMP::Lsm_DataIndex(), SLM::Lsm_DataIndex(), NOAHMP::Lsm_FluxIndex(), SLM::Lsm_FluxIndex(), obs_nudging::obs_nudging_tendency(), OpenBCCornerBox(), parse_fixed_width_int(), ERF::read_box_for_refinement(), ReadBndryPlanes::read_file(), ReadBndryPlanes::read_input_files(), ReadCustomBinaryIC(), ReadCustomDataFile(), save_qmoist_restart_vars(), ObsNudging::update_targets(), SatMethods::wv_sat_qsat_ice(), SatMethods::wv_sat_qsat_trans(), SatMethods::wv_sat_qsat_water(), SatMethods::wv_sat_svp_ice(), SatMethods::wv_sat_svp_trans(), and SatMethods::wv_sat_svp_water().

Here is the caller graph for this function:

interp_trilinear()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real interp_trilinear	(	const Real *	f,
		int	i,
		int	j,
		int	k,
		Real	tx,
		Real	ty,
		Real	tz,
		int	nx,
		int	ny,
		int	nz
	)

Trilinearly interpolate one scalar component from a flattened array.

Parameters

f	Flattened source data
i	Lower x-index of the interpolation cell
j	Lower y-index of the interpolation cell
k	Lower z-index of the interpolation cell
tx	Fractional x-coordinate within the cell
ty	Fractional y-coordinate within the cell
tz	Fractional z-coordinate within the cell
nx	Number of cells in the x-direction
ny	Number of cells in the y-direction
nz	Number of cells in the z-direction

Returns

Interpolated scalar value

{
    int i1 = amrex::min(i+1, nx-1);
    int j1 = amrex::min(j+1, ny-1);
    int k1 = amrex::min(k+1, nz-1);
 
    Real c000 = f[idx(i ,j ,k ,nx,ny)];
    Real c100 = f[idx(i1,j ,k ,nx,ny)];
    Real c010 = f[idx(i ,j1,k ,nx,ny)];
    Real c110 = f[idx(i1,j1,k ,nx,ny)];
    Real c001 = f[idx(i ,j ,k1,nx,ny)];
    Real c101 = f[idx(i1,j ,k1,nx,ny)];
    Real c011 = f[idx(i ,j1,k1,nx,ny)];
    Real c111 = f[idx(i1,j1,k1,nx,ny)];
 
    Real c00 = c000*(1-tx) + c100*tx;
    Real c10 = c010*(1-tx) + c110*tx;
    Real c01 = c001*(1-tx) + c101*tx;
    Real c11 = c011*(1-tx) + c111*tx;
 
    Real c0 = c00*(1-ty) + c10*ty;
    Real c1 = c01*(1-ty) + c11*ty;
 
    return c0*(1-tz) + c1*tz;
}

Referenced by InterpolateToFineMF().

Here is the call graph for this function:

Here is the caller graph for this function:

InterpolateToFineMF()

void InterpolateToFineMF	(	const Vector< Real > &	data_rho,
		const Vector< Real > &	data_theta,
		const Vector< Real > &	data_xvel,
		const Vector< Real > &	data_yvel,
		const Vector< Real > &	data_zvel,
		const Vector< Real > &	data_qv,
		const Vector< Real > &	data_qc,
		const Vector< Real > &	data_qrain,
		int	nx,
		int	ny,
		int	nz,
		const std::array< Real, 3 > &	problo,
		const std::array< Real, 3 > &	probhi,
		MultiFab &	mf_fine,
		const Geometry &	geom_fine
	)

Interpolate coarse custom background data onto the fine cell-centered grid.

Parameters

data_rho	Coarse density values
data_theta	Coarse potential temperature values
data_xvel	Coarse x-velocity values
data_yvel	Coarse y-velocity values
data_zvel	Coarse z-velocity values
data_qv	Coarse water vapor values
data_qc	Coarse cloud water values
data_qrain	Coarse rain water values
nx	Number of coarse cells in the x-direction
ny	Number of coarse cells in the y-direction
nz	Number of coarse cells in the z-direction
problo	Physical lower bounds of the coarse data
probhi	Physical upper bounds of the coarse data
mf_fine	Fine-grid cell-centered MultiFab to fill
geom_fine	Geometry of the fine grid

{
    // coarse spacing
    Real dx_c[3];
    dx_c[0] = (probhi[0] - problo[0]) / nx;
    dx_c[1] = (probhi[1] - problo[1]) / ny;
    dx_c[2] = (probhi[2] - problo[2]) / nz;
 
    const auto problo_f = geom_fine.ProbLoArray();
    const auto dx_f     = geom_fine.CellSizeArray();
 
    // Step 1: declare device vectors with correct size
    amrex::Gpu::DeviceVector<Real> d_rho(data_rho.size());
    amrex::Gpu::DeviceVector<Real> d_theta(data_theta.size());
    amrex::Gpu::DeviceVector<Real> d_xvel(data_xvel.size());
    amrex::Gpu::DeviceVector<Real> d_yvel(data_yvel.size());
    amrex::Gpu::DeviceVector<Real> d_zvel(data_zvel.size());
    amrex::Gpu::DeviceVector<Real> d_qv(data_qv.size());
    amrex::Gpu::DeviceVector<Real> d_qc(data_qc.size());
    amrex::Gpu::DeviceVector<Real> d_qrain(data_qrain.size());
 
    // Step 2: copy data from host to device
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_rho.begin(), data_rho.end(),
                      d_rho.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_theta.begin(), data_theta.end(),
                      d_theta.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_xvel.begin(), data_xvel.end(),
                      d_xvel.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_yvel.begin(), data_yvel.end(),
                      d_yvel.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_zvel.begin(), data_zvel.end(),
                      d_zvel.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_qv.begin(), data_qv.end(),
                      d_qv.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_qc.begin(), data_qc.end(),
                      d_qc.begin());
 
    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice,
                      data_qrain.begin(), data_qrain.end(),
                      d_qrain.begin());
 
    const Real* rho_ptr   = d_rho.data();
    const Real* theta_ptr = d_theta.data();
    const Real* xvel_ptr  = d_xvel.data();
    const Real* yvel_ptr  = d_yvel.data();
    const Real* zvel_ptr  = d_zvel.data();
    const Real* qv_ptr    = d_qv.data();
    const Real* qc_ptr    = d_qc.data();
    const Real* qrain_ptr = d_qrain.data();
 
    // -------------------------------
    // GPU kernel over MultiFab
    // -------------------------------
    for (MFIter mfi(mf_fine); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.validbox();
        auto arr = mf_fine.array(mfi);
 
        amrex::ParallelFor(bx,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // physical location (fine cell center)
            Real x = problo_f[0] + (i + myhalf) * dx_f[0];
            Real y = problo_f[1] + (j + myhalf) * dx_f[1];
            Real z = problo_f[2] + (k + myhalf) * dx_f[2];
 
            // map to coarse index space (cell centers)
            Real rx = (x - problo[0]) / dx_c[0] - myhalf;
            Real ry = (y - problo[1]) / dx_c[1] - myhalf;
            Real rz = (z - problo[2]) / dx_c[2] - myhalf;
 
            // clamp coordinates into valid coarse cell-center range
            rx = amrex::max(Real(0), amrex::min(rx, Real(nx-1)));
            ry = amrex::max(Real(0), amrex::min(ry, Real(ny-1)));
            rz = amrex::max(Real(0), amrex::min(rz, Real(nz-1)));
 
            int ic = static_cast<int>(amrex::Math::floor(rx));
            int jc = static_cast<int>(amrex::Math::floor(ry));
            int kc = static_cast<int>(amrex::Math::floor(rz));
 
            Real tx = rx - ic;
            Real ty = ry - jc;
            Real tz = rz - kc;
 
            // optional: avoid degenerate interpolation at upper edge
            if (ic == nx-1) tx = Real(0);
            if (jc == ny-1) ty = Real(0);
            if (kc == nz-1) tz = Real(0);
 
            //printf("The values are x, y, z, rx, ry, rz = %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", x, y, z, rx, ry, rz);
 
            // interpolate each component using device trilinear
            arr(i,j,k,0) = interp_trilinear(rho_ptr,   ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,1) = interp_trilinear(theta_ptr, ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,2) = interp_trilinear(xvel_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,3) = interp_trilinear(yvel_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,4) = interp_trilinear(zvel_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,5) = interp_trilinear(qv_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,6) = interp_trilinear(qc_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
            arr(i,j,k,7) = interp_trilinear(qrain_ptr,  ic,jc,kc, tx,ty,tz, nx,ny,nz);
 
            /*printf("Values are %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n",
                    arr(i,j,k,0), arr(i,j,k,1), arr(i,j,k,2), arr(i,j,k,3),
                    arr(i,j,k,4), arr(i,j,k,5), arr(i,j,k,6), arr(i,j,k,7));*/
        });
    }
}

Referenced by ERF::create_background_state_for_ensemble().

Here is the call graph for this function:

Here is the caller graph for this function:

MakeFinalMultiFabs()

void MakeFinalMultiFabs	(	const MultiFab &	mf_cc_fine,
		MultiFab &	cons_pert,
		MultiFab &	xvel_pert,
		MultiFab &	yvel_pert,
		MultiFab &	zvel_pert,
		const int	n_qstate_moist
	)

Split cell-centered interpolated background data into ERF state and face velocities.

Parameters

mf_cc_fine	Fine-grid cell-centered source data
cons_pert	Conserved-state perturbation MultiFab to fill
xvel_pert	x-face velocity perturbation MultiFab to fill
yvel_pert	y-face velocity perturbation MultiFab to fill
zvel_pert	z-face velocity perturbation MultiFab to fill

{
 
    for (MFIter mfi(cons_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
 
        auto const& mf_cc_fine_arr  = mf_cc_fine.const_array(mfi);
        auto const& cons_pert_arr = cons_pert.array(mfi);
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            Real tmp_rho   = mf_cc_fine_arr(i,j,k,0);
            Real tmp_theta = mf_cc_fine_arr(i,j,k,1);
            Real tmp_qv    = mf_cc_fine_arr(i,j,k,5);
            Real tmp_qc    = mf_cc_fine_arr(i,j,k,6);
            Real tmp_qrain = mf_cc_fine_arr(i,j,k,7);
            cons_pert_arr(i,j,k,Rho_comp)      = tmp_rho;
            cons_pert_arr(i,j,k,RhoTheta_comp) = tmp_rho*tmp_theta;
            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)  = std::max(tmp_rho*tmp_qv,0.0_rt);
            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)  = std::max(tmp_rho*tmp_qc,0.0_rt);
            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)  = std::max(tmp_rho*tmp_qrain,0.0_rt);
        });
    }
 
    // --- X-faces (component 2) ---
    for (MFIter mfi(xvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        auto const& uface = xvel_pert.array(mfi);
        auto const& cc    = mf_cc_fine.const_array(mfi);
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            uface(i,j,k) = myhalf * (cc(i-1,j,k,2) + cc(i,j,k,2));
        });
    }
 
    // --- Y-faces (component 3) ---
    for (MFIter mfi(yvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        auto const& vface = yvel_pert.array(mfi);
        auto const& cc    = mf_cc_fine.const_array(mfi);
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            vface(i,j,k) = myhalf * (cc(i,j-1,k,3) + cc(i,j,k,3));
        });
    }
 
    // --- Z-faces (component 4) ---
    for (MFIter mfi(zvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        auto const& wface = zvel_pert.array(mfi);
        auto const& cc    = mf_cc_fine.const_array(mfi);
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            wface(i,j,k) = myhalf * (cc(i,j,k-1,4) + cc(i,j,k,4));
        });
    }
}

Referenced by ERF::create_background_state_for_ensemble().

Here is the call graph for this function:

Here is the caller graph for this function:

NormalizeMultiFabRMS_PerComponent()

void NormalizeMultiFabRMS_PerComponent

(

MultiFab &

mf_cc_pert

)

Normalize each component of a perturbation MultiFab by its global RMS value.

Parameters

mf_cc_pert

MultiFab whose components are normalized in place

{
    const int ncomp = mf_cc_pert.nComp();
 
    for (int n = 0; n < ncomp; ++n)
    {
        // 1. Set up AMReX reduction (sum of squares + count)
        ReduceOps<ReduceOpSum, ReduceOpSum> reduce_op;
        ReduceData<Real, Long> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
 
        // 2. Loop over tiles and accumulate
        for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            const Box& bx = mfi.tilebox();
            auto const& arr = mf_cc_pert.const_array(mfi);
 
            reduce_op.eval(bx, reduce_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> ReduceTuple
                {
                    Real v = arr(i, j, k, n);
                    return { v * v, 1L };
                });
        }
 
        // 3. Retrieve results (includes implicit GPU sync + host copy)
        auto rv      = reduce_data.value(reduce_op);
        Real h_sumsq = amrex::get<0>(rv);
        Long h_count = amrex::get<1>(rv);
 
        // 4. Sum across MPI ranks
        ParallelDescriptor::ReduceRealSum(h_sumsq);
        ParallelDescriptor::ReduceLongSum(h_count);
 
        // 5. Compute RMS and normalize
        if (h_count > 0)
        {
            Real rms = std::sqrt(h_sumsq / static_cast<Real>(h_count));
            if (rms > zero) {
                mf_cc_pert.mult(one / rms, n, 1);
            }
        }
    }
}

Referenced by ERF::apply_gaussian_smoothing_to_perturbations().

Here is the caller graph for this function:

ReadCustomDataFile()

void ReadCustomDataFile	(	const std::string &	filename_custom,
		int &	nx,
		int &	ny,
		int &	nz,
		int &	ng,
		int &	ncomp,
		std::array< Real, 3 > &	problo_ext,
		std::array< Real, 3 > &	probhi_ext,
		Vector< Real > &	data_rho,
		Vector< Real > &	data_theta,
		Vector< Real > &	data_xvel,
		Vector< Real > &	data_yvel,
		Vector< Real > &	data_zvel,
		Vector< Real > &	data_qv,
		Vector< Real > &	data_qc,
		Vector< Real > &	data_qrain
	)

Read binary custom ensemble background data into host vectors.

Parameters

filename_custom	Path to the custom binary data file
nx	Number of cells in the x-direction read from the file
ny	Number of cells in the y-direction read from the file
nz	Number of cells in the z-direction read from the file
ng	Number of ghost cells read from the file
ncomp	Number of stored components read from the file
problo_ext	Physical lower bounds read from the file
probhi_ext	Physical upper bounds read from the file
data_rho	Density values read from the file
data_theta	Potential temperature values read from the file
data_xvel	x-velocity values read from the file
data_yvel	y-velocity values read from the file
data_zvel	z-velocity values read from the file
data_qv	Water vapor values read from the file
data_qc	Cloud water values read from the file
data_qrain	Rain water values read from the file

{
    std::ifstream ifs(filename_custom, std::ios::binary);
    if (!ifs.is_open()) {
        Abort("Failed to open file " + filename_custom + " for reading");
    }
 
    // ----------------------------
    // Read header
    // ----------------------------
    ifs.read(reinterpret_cast<char*>(&nx), sizeof(int));
    ifs.read(reinterpret_cast<char*>(&ny), sizeof(int));
    ifs.read(reinterpret_cast<char*>(&nz), sizeof(int));
 
    ifs.read(reinterpret_cast<char*>(&ng), sizeof(int));
    ifs.read(reinterpret_cast<char*>(&ncomp), sizeof(int));
 
    ifs.read(reinterpret_cast<char*>(&problo_ext[0]), sizeof(Real));
    ifs.read(reinterpret_cast<char*>(&problo_ext[1]), sizeof(Real));
    ifs.read(reinterpret_cast<char*>(&problo_ext[2]), sizeof(Real));
 
    ifs.read(reinterpret_cast<char*>(&probhi_ext[0]), sizeof(Real));
    ifs.read(reinterpret_cast<char*>(&probhi_ext[1]), sizeof(Real));
    ifs.read(reinterpret_cast<char*>(&probhi_ext[2]), sizeof(Real));
 
    const std::size_t ncell = static_cast<std::size_t>(nx) * ny * nz;
 
    // ----------------------------
    // Allocate storage
    // ----------------------------
    data_rho.resize(ncell);
    data_theta.resize(ncell);
    data_xvel.resize(ncell);
    data_yvel.resize(ncell);
    data_zvel.resize(ncell);
    data_qv.resize(ncell);
    data_qc.resize(ncell);
    data_qrain.resize(ncell);
 
    // ----------------------------
    // Read data
    // ----------------------------
    std::size_t idx = 0;
 
    for (int k = 0; k < nz; ++k)
    {
        for (int j = 0; j < ny; ++j)
        {
            for (int i = 0; i < nx; ++i)
            {
                // Skip coordinates
                Real x, y, z;
                ifs.read(reinterpret_cast<char*>(&x), sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&y), sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&z), sizeof(Real));
 
                // Read components (fixed order)
                ifs.read(reinterpret_cast<char*>(&data_rho[idx]),   sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_theta[idx]), sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_xvel[idx]),  sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_yvel[idx]),  sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_zvel[idx]),  sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_qv[idx]),  sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_qc[idx]),  sizeof(Real));
                ifs.read(reinterpret_cast<char*>(&data_qrain[idx]),  sizeof(Real));
 
                /*if(ParallelDescriptor::IOProcessor()) {
                    std::cout << "Values are " << data_rho[idx] << " " << data_theta[idx] << " "
                                               << data_xvel[idx] << " " << data_yvel[idx] << " "
                                               << data_zvel[idx] << " " << data_qv[idx] << " "
                                               << data_qc[idx] << " " << data_qrain[idx] << std::endl;
                }*/
                ++idx;
            }
        }
    }
 
    ifs.close();
}

Referenced by ERF::create_background_state_for_ensemble().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

perturb_scale

Vector<Real> perturb_scale

Initial value:

=
{
    0.01,       
    1.0,        
    1.0,        
    1.0,        
    0.5,        
    1.0e-6,     
    1.0e-6,      
    1.0e-6      
}

Referenced by AddPertToBckgnd().