io.h File Reference

Public interface for data input/output routines. More...

#include "variables.h"
#include "logging.h"
#include "Boundaries.h"

Include dependency graph for io.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
PetscErrorCode	ReadGridGenerationInputs (UserCtx *user)
	Parses command-line options for a programmatically generated grid for a SINGLE block.
PetscErrorCode	PopulateFinestUserGridResolutionFromOptions (UserCtx *finest_users, PetscInt nblk)
	Parses grid resolution arrays (-im, -jm, -km) once and applies them to all finest-grid blocks.
PetscErrorCode	ReadGridFile (UserCtx *user)
	Sets grid dimensions from a file for a SINGLE block using a one-time read cache.
PetscErrorCode	VerifyPathExistence (const char *path, PetscBool is_dir, PetscBool is_optional, const char *description, PetscBool *exists)
	A parallel-safe helper to verify the existence of a generic file or directory path.
PetscBool	ShouldWriteDataOutput (const SimCtx *simCtx, PetscInt completed_step)
	Returns whether full field/restart output should be written for the.
PetscErrorCode	ReadSimulationFields (UserCtx *user, PetscInt ti)
	Reads binary field data for velocity, pressure, and other required vectors.
PetscErrorCode	ReadFieldData (UserCtx *user, const char *field_name, Vec field_vec, PetscInt ti, const char *ext)
	Reads data for a specific field from a file into the provided vector.
PetscErrorCode	ReadStatisticalFields (UserCtx *user, PetscInt ti)
	Reads statistical fields used for time-averaged simulations.
PetscErrorCode	ReadLESFields (UserCtx *user, PetscInt ti)
	Reads LES-related fields used in turbulence modeling.
PetscErrorCode	ReadRANSFields (UserCtx *user, PetscInt ti)
	Reads RANS-related fields for turbulence modeling.
PetscErrorCode	WriteFieldData (UserCtx *user, const char *field_name, Vec field_vec, PetscInt ti, const char *ext)
	Writes data from a specific PETSc vector to a file.
PetscErrorCode	WriteSimulationFields (UserCtx *user)
	Writes simulation fields to files.
PetscErrorCode	WriteStatisticalFields (UserCtx *user)
	Writes statistical fields for averaging purposes.
PetscErrorCode	WriteLESFields (UserCtx *user)
	Writes LES-related fields.
PetscErrorCode	WriteRANSFields (UserCtx *user)
	Writes RANS-related fields.
PetscErrorCode	WriteSwarmField (UserCtx *user, const char *field_name, PetscInt ti, const char *ext)
	Writes data from a specific field in a PETSc Swarm to a file.
PetscErrorCode	WriteSwarmIntField (UserCtx *user, const char *field_name, PetscInt ti, const char *ext)
	Writes integer data from a specific PETSc Swarm field to a file.
PetscErrorCode	WriteAllSwarmFields (UserCtx *user)
	Writes a predefined set of PETSc Swarm fields to files.
PetscInt	ReadDataFileToArray (const char *filename, double **data_out, PetscInt *Nout, MPI_Comm comm)
	Reads a simple ASCII data file containing one numeric value per line.
PetscInt	CreateVTKFileFromMetadata (const char *filename, const VTKMetaData *meta, MPI_Comm comm)
	Creates a VTK file from prepared metadata and field payloads.
PetscErrorCode	VecToArrayOnRank0 (Vec inVec, PetscInt *N, double **arrayOut)
	Gathers the contents of a distributed PETSc Vec into a single array on rank 0.
PetscErrorCode	SwarmFieldToArrayOnRank0 (DM swarm, const char *field_name, PetscInt *n_total_particles, PetscInt *n_components, void **gathered_array)
	Gathers any DMSwarm field from all ranks to a single, contiguous array on rank 0.
PetscErrorCode	ReadSwarmField (UserCtx *user, const char *field_name, PetscInt ti, const char *ext)
	Reads data from a file into a specified field of a PETSc DMSwarm.
PetscErrorCode	ReadSwarmIntField (UserCtx *user, const char *field_name, PetscInt ti, const char *ext)
	Reads integer swarm data by using ReadFieldData and casting the result.
PetscErrorCode	ReadAllSwarmFields (UserCtx *user, PetscInt ti)
	Reads multiple fields (positions, velocity, CellID, and weight) into a DMSwarm.
PetscErrorCode	ReadPositionsFromFile (PetscInt timeIndex, UserCtx *user, double **coordsArray, PetscInt *Ncoords)
	Reads coordinate data (for particles) from file into a PETSc Vec, then gathers it to rank 0.
PetscErrorCode	ReadFieldDataToRank0 (PetscInt timeIndex, const char *fieldName, UserCtx *user, double **scalarArray, PetscInt *Nscalars)
	Reads a named field from file into a PETSc Vec, then gathers it to rank 0.
PetscErrorCode	DisplayBanner (SimCtx *simCtx)
	Displays a structured banner summarizing the simulation configuration.
PetscErrorCode	StringToBCFace (const char *str, BCFace *face_out)
	Converts a face-token string (e.g., "-Xi", "+Eta") to the internal BCFace enum.
PetscErrorCode	StringToBCType (const char *str, BCType *type_out)
	Converts a mathematical BC type string (e.g., "PERIODIC", "WALL") to BCType.
PetscErrorCode	StringToBCHandlerType (const char *str, BCHandlerType *handler_out)
	Converts a BC handler token (implementation strategy) to BCHandlerType.
PetscErrorCode	ValidateBCHandlerForBCType (BCType type, BCHandlerType handler)
	Validates that a selected handler is compatible with a mathematical BC type.
void	FreeBC_ParamList (BC_Param *head)
	Frees an entire linked list of boundary-condition parameters.
PetscErrorCode	GetBCParamReal (BC_Param *params, const char *key, PetscReal *value_out, PetscBool *found)
	Searches a BC_Param linked list for a key and returns its value as a double.
PetscErrorCode	GetBCParamBool (BC_Param *params, const char *key, PetscBool *value_out, PetscBool *found)
	Searches a BC_Param linked list for a key and returns its value as a bool.
PetscErrorCode	ParseAllBoundaryConditions (UserCtx *user, const char *bcs_input_filename)
	Parses the boundary conditions file to configure the type, handler, and any associated parameters for all 6 global faces of the domain.
PetscErrorCode	DeterminePeriodicity (SimCtx *simCtx)
	Scans all block-specific boundary condition files to determine a globally consistent periodicity for each dimension, reusing the core type parser.
void	TrimWhitespace (char *str)
	Helper function to trim leading/trailing whitespace from a string.
PetscErrorCode	ParsePostProcessingSettings (SimCtx *simCtx)
	Initializes post-processing settings from a config file and command-line overrides.
PetscErrorCode	ParseScalingInformation (SimCtx *simCtx)
	Parses physical scaling parameters from command-line options.

Detailed Description

Public interface for data input/output routines.

This header declares functions responsible for parsing grid geometry information, either from command-line options for programmatically generated grids or by reading the header of a grid definition file.

Definition in file io.h.

Function Documentation

ReadGridGenerationInputs()

PetscErrorCode ReadGridGenerationInputs

(

UserCtx *

user

)

Parses command-line options for a programmatically generated grid for a SINGLE block.

This function reads all per-block array options related to grid geometry, such as dimensions (-im), domain bounds (-xMins), and stretching ratios (-rxs). It then populates the fields of the provided UserCtx struct using its internal block index user->_this.

Parameters

user	Pointer to the UserCtx for a specific block. The function will populate the geometric fields (IM, Min_X, rx, etc.) within this struct.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Parses command-line options for a programmatically generated grid for a SINGLE block.

Local to this translation unit.

Definition at line 85 of file io.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
    PetscInt       nblk = simCtx->block_number;
    PetscInt       block_index = user->_this;
    PetscBool      found;
 
    // Temporary arrays to hold the parsed values for ALL blocks
    PetscInt  *IMs = NULL, *JMs = NULL, *KMs = NULL, *cgrids = NULL;
    PetscReal *xMins = NULL, *xMaxs = NULL, *rxs = NULL;
    PetscReal *yMins = NULL, *yMaxs = NULL, *rys = NULL;
    PetscReal *zMins = NULL, *zMaxs = NULL, *rzs = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Reading generated grid inputs for block %d.\n", simCtx->rank, block_index);
 
    if (block_index >= nblk) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block index %d is out of range for nblk=%d", block_index, nblk);
    }
 
    // --- Allocate temporary storage for all array options ---
    ierr = PetscMalloc4(nblk, &IMs, nblk, &JMs, nblk, &KMs, nblk, &cgrids); CHKERRQ(ierr);
    ierr = PetscMalloc6(nblk, &xMins, nblk, &xMaxs, nblk, &rxs, nblk, &yMins, nblk, &yMaxs, nblk, &rys); CHKERRQ(ierr);
    ierr = PetscMalloc3(nblk, &zMins, nblk, &zMaxs, nblk, &rzs); CHKERRQ(ierr);
 
    // --- Set default values for the temporary arrays ---
    for (PetscInt i = 0; i < nblk; ++i) {
        IMs[i] = 10; JMs[i] = 10; KMs[i] = 10; cgrids[i] = 0;
        xMins[i] = 0.0; xMaxs[i] = 1.0; rxs[i] = 1.0;
        yMins[i] = 0.0; yMaxs[i] = 1.0; rys[i] = 1.0;
        zMins[i] = 0.0; zMaxs[i] = 1.0; rzs[i] = 1.0;
    }
 
    // --- Parse the array options from the command line / control file ---
    PetscInt count;
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-im", IMs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-jm", JMs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-km", KMs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-xMins", xMins, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-xMaxs", xMaxs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-rxs", rxs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-yMins", yMins, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-yMaxs", yMaxs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-rys", rys, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-zMins", zMins, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-zMaxs", zMaxs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetRealArray(NULL, NULL, "-rzs", rzs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-cgrids", cgrids, &count, &found); CHKERRQ(ierr);
 
    // --- Assign the parsed values to the specific UserCtx struct passed in ---
    user->IM = IMs[block_index];
    user->JM = JMs[block_index];
    user->KM = KMs[block_index];
    user->Min_X = xMins[block_index];
    user->Max_X = xMaxs[block_index];
    user->rx = rxs[block_index];
    user->Min_Y = yMins[block_index];
    user->Max_Y = yMaxs[block_index];
    user->ry = rys[block_index];
    user->Min_Z = zMins[block_index];
    user->Max_Z = zMaxs[block_index];
    user->rz = rzs[block_index];
    user->cgrid = cgrids[block_index];
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Block %d grid generation inputs set: IM=%d, JM=%d, KM=%d\n",
              simCtx->rank, block_index, user->IM, user->JM, user->KM);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Block %d bounds: X=[%.2f, %.2f], Y=[%.2f, %.2f], Z=[%.2f, %.2f]\n",
              simCtx->rank, block_index, user->Min_X, user->Max_X, user->Min_Y, user->Max_Y, user->Min_Z, user->Max_Z);
 
    // --- Clean up temporary storage ---
    ierr = PetscFree4(IMs, JMs, KMs, cgrids); CHKERRQ(ierr);
    ierr = PetscFree6(xMins, xMaxs, rxs, yMins, yMaxs, rys); CHKERRQ(ierr);
    ierr = PetscFree3(zMins, zMaxs, rzs); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PopulateFinestUserGridResolutionFromOptions()

PetscErrorCode PopulateFinestUserGridResolutionFromOptions	(	UserCtx *	finest_users,
		PetscInt	nblk
	)

Parses grid resolution arrays (-im, -jm, -km) once and applies them to all finest-grid blocks.

This helper centralizes one-time resolution ingestion for analytical grid setup. It fills IM/JM/KM in each element of the finest-level UserCtx array.

Parameters

finest_users	Pointer to the finest-level UserCtx array (length nblk).
nblk	Number of blocks in the finest-level array.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Parses grid resolution arrays (-im, -jm, -km) once and applies them to all finest-grid blocks.

Local to this translation unit.

Definition at line 170 of file io.c.

{
    PetscErrorCode ierr;
    PetscBool      found;
    PetscInt       *IMs = NULL, *JMs = NULL, *KMs = NULL;
    SimCtx         *simCtx = NULL;
 
    PetscFunctionBeginUser;
 
    if (!finest_users) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "finest_users cannot be NULL.");
    }
    if (nblk <= 0) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nblk must be positive. Got %d.", nblk);
    }
    simCtx = finest_users[0].simCtx;
 
    ierr = PetscMalloc3(nblk, &IMs, nblk, &JMs, nblk, &KMs); CHKERRQ(ierr);
    for (PetscInt i = 0; i < nblk; ++i) {
        IMs[i] = 10; JMs[i] = 10; KMs[i] = 10;
    }
 
    PetscInt count;
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-im", IMs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-jm", JMs, &count, &found); CHKERRQ(ierr);
    count = nblk; ierr = PetscOptionsGetIntArray(NULL, NULL, "-km", KMs, &count, &found); CHKERRQ(ierr);
 
    for (PetscInt bi = 0; bi < nblk; ++bi) {
        finest_users[bi].IM = IMs[bi];
        finest_users[bi].JM = JMs[bi];
        finest_users[bi].KM = KMs[bi];
        if (simCtx) {
            LOG_ALLOW(LOCAL, LOG_DEBUG,
                      "Rank %d: Preloaded analytical grid resolution for block %d: IM=%d, JM=%d, KM=%d\n",
                      simCtx->rank, bi, finest_users[bi].IM, finest_users[bi].JM, finest_users[bi].KM);
        }
    }
 
    ierr = PetscFree3(IMs, JMs, KMs); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ReadGridFile()

PetscErrorCode ReadGridFile

(

UserCtx *

user

)

Sets grid dimensions from a file for a SINGLE block using a one-time read cache.

This function uses a static-variable pattern to ensure the grid file header is read only once, collectively, by all processes on the first call. Subsequent calls simply retrieve the pre-loaded and broadcasted data for the specified block.

Parameters

user	Pointer to the UserCtx for a specific block. This function will populate the IM, JM, and KM fields.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Sets grid dimensions from a file for a SINGLE block using a one-time read cache.

Local to this translation unit.

Definition at line 219 of file io.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
    PetscInt block_index = user->_this;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // --- One-Time Read and Broadcast Logic ---
    if (!g_file_has_been_read) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_INFO, "First call to ReadGridFile. Reading and broadcasting grid file header from '%s'...\n", simCtx->grid_file);
        PetscMPIInt rank = simCtx->rank;
        PetscInt    nblk = simCtx->block_number;
 
        if (rank == 0) {
            FILE *fd = fopen(simCtx->grid_file, "r");
            if (!fd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file: %s", simCtx->grid_file);
 
            // Read and validate the canonical PICGRID header.
            char firstTok[32] = {0};
            if (fscanf(fd, "%31s", firstTok) != 1)
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Empty grid file: %s", simCtx->grid_file);
            if (strcmp(firstTok, "PICGRID") != 0)
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ,
                        "Grid file %s must begin with the canonical PICGRID header.", simCtx->grid_file);
            if (fscanf(fd, "%d", &g_nblk_from_file) != 1)
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Expected number of blocks after \"PICGRID\" in %s", simCtx->grid_file);
            if (g_nblk_from_file != nblk) {
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Mismatch: -nblk is %d but grid file specifies %d blocks.", nblk, g_nblk_from_file);
            }
 
            ierr = PetscMalloc3(nblk, &g_IMs_from_file, nblk, &g_JMs_from_file, nblk, &g_KMs_from_file); CHKERRQ(ierr);
            for (PetscInt i = 0; i < nblk; ++i) {
                if (fscanf(fd, "%d %d %d\n", &g_IMs_from_file[i], &g_JMs_from_file[i], &g_KMs_from_file[i]) != 3) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Expected 3 integers for block %d in %s", i, simCtx->grid_file);
                }
            }
            fclose(fd);
        }
 
        // Broadcast nblk to verify (optional, good practice)
        ierr = MPI_Bcast(&g_nblk_from_file, 1, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        // Allocate on other ranks before receiving the broadcast
        if (rank != 0) {
            ierr = PetscMalloc3(nblk, &g_IMs_from_file, nblk, &g_JMs_from_file, nblk, &g_KMs_from_file); CHKERRQ(ierr);
        }
        
        // Broadcast the data arrays
        ierr = MPI_Bcast(g_IMs_from_file, nblk, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Bcast(g_JMs_from_file, nblk, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Bcast(g_KMs_from_file, nblk, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        g_file_has_been_read = PETSC_TRUE;
        LOG_ALLOW(GLOBAL, LOG_INFO, "Grid file header read and broadcast complete.\n");
    }
 
    // --- Per-Block Assignment Logic (runs on every call) ---
    user->IM = g_IMs_from_file[block_index];
    user->JM = g_JMs_from_file[block_index];
    user->KM = g_KMs_from_file[block_index];
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Set file inputs for Block %d: IM=%d, JM=%d, KM=%d\n",
              simCtx->rank, block_index, user->IM, user->JM, user->KM);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

VerifyPathExistence()

PetscErrorCode VerifyPathExistence	(	const char *	path,
		PetscBool	is_dir,
		PetscBool	is_optional,
		const char *	description,
		PetscBool *	exists
	)

A parallel-safe helper to verify the existence of a generic file or directory path.

This function centralizes the logic for checking arbitrary paths. Only Rank 0 performs the filesystem check, and the result is broadcast to all other processes. This ensures collective and synchronized decision-making across all ranks. It is intended for configuration files, source directories, etc., where the path is known completely.

Parameters

[in]	path	The full path to the file or directory to check.
[in]	is_dir	PETSC_TRUE if checking for a directory, PETSC_FALSE for a file.
[in]	is_optional	PETSC_TRUE if the path is optional (results in a warning), PETSC_FALSE if mandatory (results in an error).
[in]	description	A user-friendly description of the path for logging (e.g., "Grid file").
[out]	exists	The result of the check (identical on all ranks).

Returns

PetscErrorCode

A parallel-safe helper to verify the existence of a generic file or directory path.

Local to this translation unit.

Definition at line 741 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    MPI_Comm       comm = PETSC_COMM_WORLD;
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr);
 
    if (rank == 0) {
        if (is_dir) {
            ierr = PetscTestDirectory(path, 'r', exists); CHKERRQ(ierr);
        } else {
            ierr = PetscTestFile(path, 'r', exists); CHKERRQ(ierr);
        }
 
        if (!(*exists)) {
            if (is_optional) {
                LOG_ALLOW(GLOBAL, LOG_WARNING, "Optional %s not found at: %s (using defaults/ignoring).\n", description, path);
            } else {
                LOG_ALLOW(GLOBAL, LOG_ERROR, "Mandatory %s not found at: %s\n", description, path);
            }
        } else {
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "Found %s: %s\n", description, path);
        }
    }
 
    // Broadcast the result from Rank 0
    PetscMPIInt exists_int = (rank == 0) ? (PetscMPIInt)(*exists) : 0;
    ierr = MPI_Bcast(&exists_int, 1, MPI_INT, 0, comm); CHKERRMPI(ierr);
    *exists = (PetscBool)exists_int;
 
    // Collective error for mandatory files
    if (!(*exists) && !is_optional) {
        SETERRQ(comm, PETSC_ERR_FILE_OPEN, "Mandatory %s not found. Rank 0 expected it at '%s'. Check path and permissions.", description, path);
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ShouldWriteDataOutput()

PetscBool ShouldWriteDataOutput	(	const SimCtx *	simCtx,
		PetscInt	completed_step
	)

Returns whether full field/restart output should be written for the.

completed timestep.

Parameters

simCtx	Simulation context controlling the operation.
completed_step	Completed step index used by the decision helper.

Returns

PetscBool indicating the result of ShouldWriteDataOutput().

Returns whether full field/restart output should be written for the.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

ShouldWriteDataOutput()

Definition at line 70 of file io.c.

{
    if (!simCtx) {
        return PETSC_FALSE;
    }
    return (PetscBool)(simCtx->tiout > 0 && completed_step > 0 && completed_step % simCtx->tiout == 0);
}

Here is the caller graph for this function:

ReadSimulationFields()

PetscErrorCode ReadSimulationFields	(	UserCtx *	user,
		PetscInt	ti
	)

Reads binary field data for velocity, pressure, and other required vectors.

Reads contravariant velocity (Ucont) from vfield, Cartesian velocity (Ucat) from ufield, pressure (P), node state (Nvert_o), and optionally statistical quantities, LES, and RANS fields from binary files. Logs missing files but continues execution.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the simulation context.
[in]	ti	Timestep index used to construct restart file names.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads binary field data for velocity, pressure, and other required vectors.

Local to this translation unit.

Definition at line 1129 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
    const char *source_path = NULL;
    const char *eulerian_ext = "dat";
 
    if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
        source_path = simCtx->pps->source_dir;
        if (simCtx->pps->eulerianExt[0] != '\0') {
            eulerian_ext = simCtx->pps->eulerianExt;
        }
    } else if(simCtx->exec_mode == EXEC_MODE_SOLVER){
        source_path = simCtx->restart_dir;
    } else{
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Invalid execution mode for reading simulation fields.");
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read simulation fields.\n");
 
    // Set the current I/O directory context
    ierr  = PetscSNPrintf(simCtx->_io_context_buffer, sizeof(simCtx->_io_context_buffer), "%s/%s",source_path,simCtx->euler_subdir); CHKERRQ(ierr);
    simCtx->current_io_directory = simCtx->_io_context_buffer;
 
    // Read Cartesian velocity field
    ierr = ReadFieldData(user, "ufield", user->Ucat, ti, eulerian_ext); CHKERRQ(ierr);
 
    // Read contravariant velocity field
    ierr = ReadFieldData(user, "vfield", user->Ucont, ti, eulerian_ext); CHKERRQ(ierr);
 
    // Read pressure field
    ierr = ReadFieldData(user, "pfield", user->P, ti, eulerian_ext); CHKERRQ(ierr);
 
    // Read node state field (nvert)
    ierr = ReadFieldData(user, "nvfield", user->Nvert, ti, eulerian_ext); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL,LOG_INFO,"Successfully read all mandatory fields. \n");
 
    if(simCtx->np>0){    
    // Read Particle Count field
    if(!user->ParticleCount){
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "ParticleCount Vec is NULL but np>0");
    }
    ierr = ReadOptionalField(user, "ParticleCount", "Particle Count", user->ParticleCount, ti, eulerian_ext); CHKERRQ(ierr);
    
    ierr = ReadOptionalField(user, "psifield", "Scalar Psi Field", user->Psi, ti, eulerian_ext); CHKERRQ(ierr);
    }
    else{
        LOG_ALLOW(GLOBAL, LOG_INFO, "No particles in simulation, skipping Particle fields reading.\n");
    }
    // Process LES fields if enabled
    if (simCtx->les) {
      ierr = ReadLESFields(user,ti); CHKERRQ(ierr);
    }
 
    // Process RANS fields if enabled
    if (simCtx->rans) {
      ierr = ReadRANSFields(user,ti); CHKERRQ(ierr);
    }
 
    // Process statistical fields if averaging is enabled
    if (simCtx->averaging) {
      ierr = ReadStatisticalFields(user,ti); CHKERRQ(ierr);
    }
 
    simCtx->current_io_directory = NULL; // Clear the I/O context after reading
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished reading simulation fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadFieldData()

PetscErrorCode ReadFieldData	(	UserCtx *	user,
		const char *	field_name,
		Vec	field_vec,
		PetscInt	ti,
		const char *	ext
	)

Reads data for a specific field from a file into the provided vector.

This function uses the field name to construct the file path and reads the data from the corresponding file into the provided PETSc vector.

Parameters

[in]	user	Pointer to the UserCtx structure containing simulation context.
[in]	field_name	Name of the field (e.g., "ufield", "vfield", "pfield").
[out]	field_vec	PETSc vector to store the field data.
[in]	ti	Time index for constructing the file name.
[in]	ext	File extension (e.g., "dat").

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads data for a specific field from a file into the provided vector.

Local to this translation unit.

Definition at line 901 of file io.c.

{
   PetscErrorCode ierr;
   char           filename[PETSC_MAX_PATH_LEN];
   MPI_Comm       comm;
   PetscMPIInt    rank,size;
   SimCtx         *simCtx = user->simCtx;
 
 
   PetscFunctionBeginUser;
   PROFILE_FUNCTION_BEGIN;
 
    if(!simCtx->current_io_directory){
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "I/O context directory was not set before calling ReadFieldData().");
    }
 
 
   ierr = PetscObjectGetComm((PetscObject)field_vec,&comm);CHKERRQ(ierr);
   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
 
   const char *source_path = NULL;
   source_path = simCtx->current_io_directory;
 
   if(!source_path){
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "source_path was not set for the current execution mode.");
   }
   /* ---------------------------------------------------------------------
    * Compose <path-to-file>/<field_name><step with 5 digits>_0.<ext>
    * (all restart files are written by rank-0 with that naming scheme).
    * ------------------------------------------------------------------ */
   ierr = PetscSNPrintf(filename,sizeof(filename),
                        "%s/%s%05" PetscInt_FMT "_0.%s",
                        source_path,field_name,ti,ext);CHKERRQ(ierr);
 
   LOG_ALLOW(GLOBAL,LOG_DEBUG,
             "Attempting to read <%s> on rank %d/%d\n",
             filename,(int)rank,(int)size);
 
   /* ======================================================================
    * 1.  SERIAL JOB  – just hand the Vec to VecLoad()
    * ==================================================================== */
   if(size==1)
   {
      PetscViewer viewer;
      PetscBool   found;
      Vec temp_vec;
      PetscInt expectedSize,loadedSize;
 
      ierr = PetscTestFile(filename,'r',&found);CHKERRQ(ierr);
      if(!found) SETERRQ(comm,PETSC_ERR_FILE_OPEN,
                         "Restart/Source file not found: %s",filename);
 
      ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
// ---- START MODIFICATION ----
      // DO NOT load directly into field_vec, as this can resize it, which is
      // illegal for DMSwarm "view" vectors. Instead, load into a temporary vector.
      ierr = VecCreate(PETSC_COMM_SELF, &temp_vec); CHKERRQ(ierr);
      ierr = VecLoad(temp_vec,viewer);CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
 
      // Sanity check: ensure the file size matches the expected vector size.
      ierr = VecGetSize(field_vec, &expectedSize);CHKERRQ(ierr);
      ierr = VecGetSize(temp_vec, &loadedSize);CHKERRQ(ierr);
      if (loadedSize != expectedSize) {
         SETERRQ(comm,PETSC_ERR_FILE_UNEXPECTED,
                 "File %s holds %d entries – expected %d for field '%s'",
                 filename, loadedSize, expectedSize, field_name);
      }
 
      // Now, safely copy the data from the temporary vector to the final destination.
      ierr = VecCopy(temp_vec, field_vec);CHKERRQ(ierr);
 
      // Clean up the temporary vector.
      ierr = VecDestroy(&temp_vec);CHKERRQ(ierr);
 
      // ---- END MODIFICATION ----
      
      /* create EMPTY sequential Vec – VecLoad() will size it correctly   */
      /*
      ierr = VecCreate(PETSC_COMM_SELF,&seq_vec);CHKERRQ(ierr);
      ierr = VecSetType(seq_vec,VECSEQ);CHKERRQ(ierr);
 
      ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,filename,
                                   FILE_MODE_READ,&viewer);CHKERRQ(ierr);
 
      ierr = VecLoad(field_vec,viewer);CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
      */
      LOG_ALLOW(GLOBAL,LOG_INFO,
                "Loaded <%s> (serial path)\n",filename);
 
      PROFILE_FUNCTION_END;          
      PetscFunctionReturn(0);
   }
 
   /* ======================================================================
    * 2.  PARALLEL JOB
    * ==================================================================== */
   PetscInt globalSize;
   ierr = VecGetSize(field_vec,&globalSize);CHKERRQ(ierr);
 
   DM         dm = NULL;
   const char *dmtype = NULL;
   Vec        nat = NULL;                 /* Natural-ordered vector for DMDA */
 
   /* -------------------- rank-0 : read the sequential file -------------- */
   Vec            seq_vec = NULL;      /* only valid on rank-0            */
   const PetscScalar *seqArray = NULL; /* borrowed pointer on rank-0 only */
 
   if(rank==0)
   {
      PetscViewer viewer;
      PetscBool   found;
 
      ierr = PetscTestFile(filename,'r',&found);CHKERRQ(ierr);
      if(!found) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,
                         "Restart file not found: %s",filename);
 
      /* create EMPTY sequential Vec – VecLoad() will size it correctly   */
      ierr = VecCreate(PETSC_COMM_SELF,&seq_vec);CHKERRQ(ierr);
      ierr = VecSetType(seq_vec,VECSEQ);CHKERRQ(ierr);
 
      ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,filename,
                                   FILE_MODE_READ,&viewer);CHKERRQ(ierr);
      ierr = VecLoad(seq_vec,viewer);CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
 
      /* size sanity-check */
      PetscInt loaded;
      ierr = VecGetSize(seq_vec,&loaded);CHKERRQ(ierr);
      if(loaded != globalSize)
         SETERRQ(comm,PETSC_ERR_FILE_UNEXPECTED,
                 "File %s holds %d entries – expected %d",
                 filename,loaded,globalSize);
 
      /* borrow array for later Bcast */
      ierr = VecGetArrayRead(seq_vec,&seqArray);CHKERRQ(ierr);
 
      LOG_ALLOW(GLOBAL,LOG_TRACE,
                "Rank 0 successfully loaded <%s>\n",filename);
   }
 
   /* -------------------- Check if this is a DMDA vector ----------------- */
   ierr = VecGetDM(field_vec, &dm); CHKERRQ(ierr);
   if (dm) { ierr = DMGetType(dm, &dmtype); CHKERRQ(ierr); }
 
   if (dmtype && !strcmp(dmtype, DMDA)) {
      /* ==================================================================
       * DMDA PATH: File is in natural ordering, need to convert to global
       * ================================================================== */
 
      /* Create natural vector */
      ierr = DMDACreateNaturalVector(dm, &nat); CHKERRQ(ierr);
 
      /* Scatter from rank 0's seq_vec to all ranks' natural vector */
      VecScatter scatter;
      Vec nat_seq = NULL;  /* Sequential natural vector on rank 0 */
 
      ierr = VecScatterCreateToZero(nat, &scatter, &nat_seq); CHKERRQ(ierr);
 
      /* Reverse scatter: from rank 0 to all ranks */
      ierr = VecScatterBegin(scatter, (rank == 0 ? seq_vec : nat_seq), nat,
                             INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr);
      ierr = VecScatterEnd(scatter, (rank == 0 ? seq_vec : nat_seq), nat,
                           INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr);
 
      /* Convert natural → global ordering */
      ierr = DMDANaturalToGlobalBegin(dm, nat, INSERT_VALUES, field_vec); CHKERRQ(ierr);
      ierr = DMDANaturalToGlobalEnd(dm, nat, INSERT_VALUES, field_vec); CHKERRQ(ierr);
 
      /* Cleanup */
      ierr = VecScatterDestroy(&scatter); CHKERRQ(ierr);
      ierr = VecDestroy(&nat_seq); CHKERRQ(ierr);
      ierr = VecDestroy(&nat); CHKERRQ(ierr);
 
   } else {
      /* ==================================================================
       * NON-DMDA PATH: Use broadcast and direct copy (assumes global ordering)
       * ================================================================== */
 
      PetscScalar *buffer = NULL;
      if (rank == 0) {
         buffer = (PetscScalar *)seqArray;
      } else {
         ierr = PetscMalloc1(globalSize, &buffer); CHKERRQ(ierr);
      }
 
      ierr = MPI_Bcast(buffer, (int)globalSize, MPIU_SCALAR, 0, comm); CHKERRQ(ierr);
 
      /* Copy slice based on ownership range */
      PetscInt rstart, rend, loc;
      PetscScalar *locArray;
 
      ierr = VecGetOwnershipRange(field_vec, &rstart, &rend); CHKERRQ(ierr);
      loc = rend - rstart;
 
      ierr = VecGetArray(field_vec, &locArray); CHKERRQ(ierr);
      ierr = PetscMemcpy(locArray, buffer + rstart, loc * sizeof(PetscScalar)); CHKERRQ(ierr);
      ierr = VecRestoreArray(field_vec, &locArray); CHKERRQ(ierr);
 
      if (rank != 0) {
         ierr = PetscFree(buffer); CHKERRQ(ierr);
      }
   }
 
   /* -------------------- tidy up ---------------------------------------- */
   if (rank == 0) {
      ierr = VecRestoreArrayRead(seq_vec, &seqArray); CHKERRQ(ierr);
      ierr = VecDestroy(&seq_vec); CHKERRQ(ierr);
   }
 
   LOG_ALLOW(GLOBAL,LOG_INFO,
             "Loaded <%s> (parallel path)\n",filename);
 
   PROFILE_FUNCTION_END;          
   PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ReadStatisticalFields()

PetscErrorCode ReadStatisticalFields	(	UserCtx *	user,
		PetscInt	ti
	)

Reads statistical fields used for time-averaged simulations.

Reads statistical quantities such as velocity sums and pressure sums. Logs missing files and initializes fields to zero if files are unavailable.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the simulation context.
[in]	ti	Timestep index used to construct statistics file names.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads statistical fields used for time-averaged simulations.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

ReadStatisticalFields()

Definition at line 1208 of file io.c.

{
    PetscErrorCode ierr;
    const char *eulerian_ext = "dat";
    if (user->simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR && user->simCtx->pps->eulerianExt[0] != '\0') {
        eulerian_ext = user->simCtx->pps->eulerianExt;
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read statistical fields.\n");
 
    ierr = ReadOptionalField(user, "su0", "Velocity Sum",       user->Ucat_sum,        ti, eulerian_ext); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "su1", "Velocity Cross Sum", user->Ucat_cross_sum,  ti, eulerian_ext); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "su2", "Velocity Square Sum",user->Ucat_square_sum, ti, eulerian_ext); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "sp",  "Pressure Sum",       user->P_sum,           ti, eulerian_ext); CHKERRQ(ierr);    
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished reading statistical fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadLESFields()

PetscErrorCode ReadLESFields	(	UserCtx *	user,
		PetscInt	ti
	)

Reads LES-related fields used in turbulence modeling.

Reads the Smagorinsky constant (Cs) and transfers data to local vectors. Logs missing files and initializes fields to zero if files are unavailable.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the simulation context.
[in]	ti	Time index for constructing the file name.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads LES-related fields used in turbulence modeling.

Local to this translation unit.

Definition at line 1232 of file io.c.

{
    PetscErrorCode ierr;
    const char *eulerian_ext = "dat";
    if (user->simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR && user->simCtx->pps->eulerianExt[0] != '\0') {
        eulerian_ext = user->simCtx->pps->eulerianExt;
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read LES fields.\n");
 
    ierr = ReadOptionalField(user, "Nu_t", "Turbulent Viscosity", user->Nu_t, ti, eulerian_ext); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "cs", "Smagorinsky Constant (Cs)", user->CS, ti, eulerian_ext); CHKERRQ(ierr);
 
    ierr = UpdateLocalGhosts(user, "CS"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Nu_t"); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished reading LES fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadRANSFields()

PetscErrorCode ReadRANSFields	(	UserCtx *	user,
		PetscInt	ti
	)

Reads RANS-related fields for turbulence modeling.

Reads K_Omega fields (used in RANS modeling) and initializes them if files are unavailable.

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the simulation context.
[in]	ti	Time index for constructing the file name.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads RANS-related fields for turbulence modeling.

Local to this translation unit.

Definition at line 1257 of file io.c.

{
    PetscErrorCode ierr;
    const char *eulerian_ext = "dat";
    if (user->simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR && user->simCtx->pps->eulerianExt[0] != '\0') {
        eulerian_ext = user->simCtx->pps->eulerianExt;
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read RANS fields.\n");
 
    ierr = ReadOptionalField(user, "kfield", "K-Omega RANS", user->K_Omega, ti, eulerian_ext); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "Nu_t", "Turbulent Viscosity", user->Nu_t, ti, eulerian_ext); CHKERRQ(ierr);
 
    VecCopy(user->K_Omega, user->K_Omega_o);
 
    ierr = UpdateLocalGhosts(user, "K_Omega"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "Nu_t"); CHKERRQ(ierr);
    ierr = UpdateLocalGhosts(user, "K_Omega_o"); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished reading RANS fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteFieldData()

PetscErrorCode WriteFieldData	(	UserCtx *	user,
		const char *	field_name,
		Vec	field_vec,
		PetscInt	ti,
		const char *	ext
	)

Writes data from a specific PETSc vector to a file.

This function uses the field name to construct the file path and writes the data from the provided PETSc vector to the corresponding file.

Parameters

[in]	user	Pointer to the UserCtx structure containing simulation context.
[in]	field_name	Name of the field (e.g., "ufield", "vfield", "pfield").
[in]	field_vec	PETSc vector containing the field data to write.
[in]	ti	Time index for constructing the file name.
[in]	ext	File extension (e.g., "dat").

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes data from a specific PETSc vector to a file.

Local to this translation unit.

Definition at line 1443 of file io.c.

{
    PetscErrorCode ierr;
    MPI_Comm       comm;
    PetscMPIInt    rank, size;
 
    const PetscInt placeholder_int = 0;                    /* keep legacy name */
    char           filename[PETSC_MAX_PATH_LEN];
    SimCtx  *simCtx=user->simCtx;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if(!simCtx->current_io_directory){
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "I/O context directory was not set before calling WriteFieldData().");
    }
 
    /* ------------------------------------------------------------ */
    /*                  Basic communicator information              */
    /* ------------------------------------------------------------ */
    ierr = PetscObjectGetComm((PetscObject)field_vec,&comm);CHKERRQ(ierr);
    ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
    ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
 
    ierr = PetscSNPrintf(filename,sizeof(filename),
                         "%s/%s%05" PetscInt_FMT "_%d.%s",
                         simCtx->current_io_directory,field_name,ti,placeholder_int,ext);CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL,LOG_DEBUG,
              " Preparing to write <%s> on rank %d/%d\n",
              filename,rank,size);
 
    /* ------------------------------------------------------------ */
    /*                       1.  Serial path                        */
    /* ------------------------------------------------------------ */
    if (size == 1) {
        PetscViewer viewer;
 
        ierr = PetscViewerBinaryOpen(comm,filename,
                                     FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
        ierr = VecView(field_vec,viewer);CHKERRQ(ierr);
        ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
 
        LOG_ALLOW(GLOBAL,LOG_INFO,
                  " Wrote <%s> (serial path)\n",filename);
        PetscFunctionReturn(0);
    }
 
    /* ------------------------------------------------------------ */
    /*                      2.  Parallel path                       */
    /* ------------------------------------------------------------ */
    VecScatter scatter;
    Vec        seq_vec=NULL;               /* created by PETSc, lives only on rank 0 */
    DM         dm = NULL;
    const char *dmtype = NULL;
    Vec        nat = NULL;                 /* Natural-ordered vector for DMDA */
 
    /* 2.1  Check if vector has DMDA and convert to natural ordering if needed */
    ierr = VecGetDM(field_vec, &dm); CHKERRQ(ierr);
    if (dm) { ierr = DMGetType(dm, &dmtype); CHKERRQ(ierr); }
 
    if (dmtype && !strcmp(dmtype, DMDA)) {
        /* DMDA path: convert to natural ordering first */
        ierr = DMDACreateNaturalVector(dm, &nat); CHKERRQ(ierr);
        ierr = DMDAGlobalToNaturalBegin(dm, field_vec, INSERT_VALUES, nat); CHKERRQ(ierr);
        ierr = DMDAGlobalToNaturalEnd(dm, field_vec, INSERT_VALUES, nat); CHKERRQ(ierr);
 
        /* Gather the natural-ordered vector */
        ierr = VecScatterCreateToZero(nat, &scatter, &seq_vec); CHKERRQ(ierr);
    } else {
        /* Non-DMDA path: direct gather in global ordering */
        ierr = VecScatterCreateToZero(field_vec, &scatter, &seq_vec); CHKERRQ(ierr);
    }
 
    /* 2.2  Gather distributed → sequential (on rank 0)             */
    ierr = VecScatterBegin(scatter, (nat ? nat : field_vec), seq_vec,
                           INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecScatterEnd(scatter, (nat ? nat : field_vec), seq_vec,
                         INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
 
    /* 2.3  Rank 0 writes the file                                  */
    if (rank == 0) {
        PetscViewer viewer;
 
        /* (optional) value diagnostics */
        PetscReal vmin,vmax;
        ierr = VecMin(seq_vec,NULL,&vmin);CHKERRQ(ierr);
        ierr = VecMax(seq_vec,NULL,&vmax);CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL,LOG_DEBUG,
                  " <%s> range = [%.4e … %.4e]\n",
                  field_name,(double)vmin,(double)vmax);
 
        ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,filename,
                                     FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
        ierr = VecView(seq_vec,viewer);CHKERRQ(ierr);
        ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
 
        LOG_ALLOW(GLOBAL,LOG_INFO,
                  " Wrote <%s> (parallel path)\n",filename);
    }
 
    /* 2.4  Cleanup                                                 */
    ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr);
    ierr = VecDestroy(&seq_vec);CHKERRQ(ierr);
    if (nat) { ierr = VecDestroy(&nat); CHKERRQ(ierr); }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

WriteSimulationFields()

PetscErrorCode WriteSimulationFields

(

UserCtx *

user

)

Writes simulation fields to files.

This function writes contravariant velocity, Cartesian velocity, pressure, and node state fields to their respective binary files. It also conditionally writes LES, RANS, and statistical fields if they are enabled.

Parameters

[in]

user

Pointer to the UserCtx structure containing simulation context.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes simulation fields to files.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

WriteSimulationFields()

Definition at line 1563 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to write simulation fields.\n");
 
    // Set the current IO directory
    ierr = PetscSNPrintf(simCtx->_io_context_buffer, sizeof(simCtx->_io_context_buffer),
                         "%s/%s", simCtx->output_dir, simCtx->euler_subdir);CHKERRQ(ierr);
    simCtx->current_io_directory = simCtx->_io_context_buffer;
 
    // Write contravariant velocity field
    ierr = WriteFieldData(user, "vfield", user->Ucont, simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write Cartesian velocity field
    ierr = WriteFieldData(user, "ufield", user->Ucat, simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write pressure field
    ierr = WriteFieldData(user, "pfield", user->P, simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write node state field (nvert)
    ierr = WriteFieldData(user, "nvfield", user->Nvert, simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write ParticleCountPerCell if enabled.
    if(simCtx->np>0){
    ierr = WriteFieldData(user, "ParticleCount",user->ParticleCount,simCtx->step,"dat"); CHKERRQ(ierr);
    ierr = WriteFieldData(user, "psifield", user->Psi, simCtx->step, "dat"); CHKERRQ(ierr);
    }
    
    // Write LES fields if enabled
    if (simCtx->les) {
        ierr = WriteLESFields(user); CHKERRQ(ierr);
    }
 
    // Write RANS fields if enabled
    if (simCtx->rans) {
        ierr = WriteRANSFields(user); CHKERRQ(ierr);
    }
 
    // Write statistical fields if averaging is enabled
    if (simCtx->averaging) {
        ierr = WriteStatisticalFields(user); CHKERRQ(ierr);
    }
 
    simCtx->current_io_directory = NULL;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished writing simulation fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteStatisticalFields()

PetscErrorCode WriteStatisticalFields

(

UserCtx *

user

)

Writes statistical fields for averaging purposes.

This function writes data for fields such as Ucat_sum, Ucat_cross_sum, Ucat_square_sum, and P_sum to their respective binary files.

Parameters

[in]

user

Pointer to the UserCtx structure containing simulation context.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes statistical fields for averaging purposes.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

WriteStatisticalFields()

Definition at line 1622 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to write statistical fields.\n");
 
    ierr = WriteFieldData(user, "su0", user->Ucat_sum, simCtx->step, "dat"); CHKERRQ(ierr);
    ierr = WriteFieldData(user, "su1", user->Ucat_cross_sum, simCtx->step, "dat"); CHKERRQ(ierr);
    ierr = WriteFieldData(user, "su2", user->Ucat_square_sum, simCtx->step, "dat"); CHKERRQ(ierr);
    ierr = WriteFieldData(user, "sp", user->P_sum, simCtx->step, "dat"); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished writing statistical fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteLESFields()

PetscErrorCode WriteLESFields

(

UserCtx *

user

)

Writes LES-related fields.

This function writes LES-related fields such as Cs (Smagorinsky constant) to their respective binary files.

Parameters

[in]

user

Pointer to the UserCtx structure containing simulation context.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes LES-related fields.

Local to this translation unit.

Definition at line 1644 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to write LES fields.\n");
 
 
    DMLocalToGlobalBegin(user->da, user->lCs, INSERT_VALUES, user->CS);
    DMLocalToGlobalEnd(user->da, user->lCs, INSERT_VALUES, user->CS);
 
    DMLocalToGlobalBegin(user->da, user->lNu_t, INSERT_VALUES, user->Nu_t);
    DMLocalToGlobalEnd(user->da, user->lNu_t, INSERT_VALUES, user->Nu_t);
 
    ierr = WriteFieldData(user, "Nu_t", user->Nu_t, simCtx->step, "dat"); CHKERRQ(ierr);
    ierr = WriteFieldData(user, "cs", user->CS, simCtx->step, "dat"); CHKERRQ(ierr);
   
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished writing LES fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteRANSFields()

PetscErrorCode WriteRANSFields

(

UserCtx *

user

)

Writes RANS-related fields.

This function writes RANS-related fields such as K_Omega to their respective binary files.

Parameters

[in]

user

Pointer to the UserCtx structure containing simulation context.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes RANS-related fields.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

WriteRANSFields()

Definition at line 1674 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to write RANS fields.\n");
 
    ierr = WriteFieldData(user, "kfield", user->K_Omega, simCtx->step, "dat"); CHKERRQ(ierr);
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished writing RANS fields.\n");
 
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteSwarmField()

PetscErrorCode WriteSwarmField	(	UserCtx *	user,
		const char *	field_name,
		PetscInt	ti,
		const char *	ext
	)

Writes data from a specific field in a PETSc Swarm to a file.

This function retrieves the Swarm from the UserCtx (i.e., user->swarm) and creates a global PETSc vector from the specified Swarm field. It then calls the existing WriteFieldData() function to handle the actual I/O operation. After writing the data, the function destroys the temporary global vector to avoid memory leaks.

Parameters

[in]	user	Pointer to the UserCtx structure containing simulation context and the PetscSwarm (as user->swarm).
[in]	field_name	Name of the Swarm field to be written (e.g., "my_field").
[in]	ti	Time index used to construct the output file name.
[in]	ext	File extension (e.g., "dat", "bin").

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Note

Compatible with PETSc 3.14.4.

Writes data from a specific field in a PETSc Swarm to a file.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

WriteSwarmField()

Definition at line 1695 of file io.c.

{
  PetscErrorCode ierr;
  Vec            fieldVec;
  DM     swarm;  
 
  PetscFunctionBeginUser; /* PETSc macro indicating start of function */
 
  /* 
   * 1) Retrieve the PetscSwarm from the user context.
   *    Ensure user->swarm is initialized and not NULL.
   */
  swarm = user->swarm;      
 
  /* 
   * 2) Create a global vector from the specified swarm field.
   *    This function is available in PETSc 3.14.4.
   *    It provides a read/write "view" of the swarm field as a global Vec.
   */
  LOG_ALLOW(GLOBAL, LOG_DEBUG,
            "Attempting to create global vector from field: %s\n",
            field_name);
  ierr = DMSwarmCreateGlobalVectorFromField(swarm, field_name, &fieldVec);CHKERRQ(ierr);
 
  /*
   * 3) Use your existing WriteFieldData() to write the global vector to a file.
   *    The field name, time index, and extension are passed along for naming.
   */
  LOG_ALLOW(GLOBAL, LOG_DEBUG,
            "Calling WriteFieldData for field: %s\n",
            field_name);
  ierr = WriteFieldData(user, field_name, fieldVec, ti, ext);CHKERRQ(ierr);
 
  /*
   * 4) Destroy the global vector once the data is successfully written.
   *    This step is crucial for avoiding memory leaks. 
   *    DMSwarmDestroyGlobalVectorFromField() is also available in PETSc 3.14.4.
   */
  LOG_ALLOW(GLOBAL, LOG_DEBUG,
            "Destroying the global vector for field: %s\n",
            field_name);
  ierr = DMSwarmDestroyGlobalVectorFromField(swarm, field_name, &fieldVec);CHKERRQ(ierr);
 
  /* Log and return success. */
  LOG_ALLOW(GLOBAL, LOG_INFO,
            "Successfully wrote swarm data for field: %s\n",
            field_name);
 
  PetscFunctionReturn(0); /* PETSc macro indicating end of function */
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteSwarmIntField()

PetscErrorCode WriteSwarmIntField	(	UserCtx *	user,
		const char *	field_name,
		PetscInt	ti,
		const char *	ext
	)

Writes integer data from a specific PETSc Swarm field to a file.

This function is designed for swarm fields that store integer data (e.g., DMSwarm_CellID), which cannot be converted to a standard PETSc Vec of PetscScalars. It accesses the raw data pointer for the field on each rank using DMSwarmGetField(), writes the local data to a rank-specific binary file, and then restores the field access.

Parameters

[in]	user	Pointer to the UserCtx structure containing the PetscSwarm.
[in]	field_name	Name of the integer Swarm field to be written.
[in]	ti	Time index used to construct the output file name.
[in]	ext	File extension (e.g., "dat", "bin").

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes integer data from a specific PETSc Swarm field to a file.

Local to this translation unit.

Definition at line 1750 of file io.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    Vec            temp_vec;       // Temporary Vec to hold casted data
    PetscInt       nlocal, nglobal,bs,i;
    void           *field_array_void;
    PetscScalar    *scalar_array;  // Pointer to the temporary Vec's scalar data
 
    PetscFunctionBeginUser;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Casting '%s' to Vec for writing.\n", field_name);
 
    // Get the swarm field properties
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    ierr = DMSwarmGetSize(swarm, &nglobal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, field_name, &bs, NULL, &field_array_void); CHKERRQ(ierr);
 
    // Create Temporary parallel Vec wit the CORRECT layout
    ierr = VecCreate(PETSC_COMM_WORLD, &temp_vec); CHKERRQ(ierr);
    ierr = VecSetType(temp_vec, VECMPI); CHKERRQ(ierr);
    ierr = VecSetSizes(temp_vec, nlocal*bs, nglobal*bs); CHKERRQ(ierr);
    ierr = VecSetUp(temp_vec); CHKERRQ(ierr);
 
    // Defining Vector field to mandatory field 'position'
    DMSwarmVectorDefineField(swarm,"position");
              
    ierr = VecGetArray(temp_vec, &scalar_array); CHKERRQ(ierr);
 
    if(strcasecmp(field_name,"DMSwarm_pid") == 0){
        PetscInt64 *int64_array = (PetscInt64 *)field_array_void;
        // Perform the cast from PetscInt64 to PetscScalar
        for (i = 0; i < nlocal*bs; i++) {
            scalar_array[i] = (PetscScalar)int64_array[i];
        }
    }else{
        PetscInt *int_array = (PetscInt *)field_array_void;
        //Perform the cast from PetscInt to PetscScalar
        for (i = 0; i < nlocal*bs; i++) {
            scalar_array[i] = (PetscScalar)int_array[i];
        }
    }
 
    // Restore access to both arrays
    ierr = VecRestoreArray(temp_vec, &scalar_array); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, field_name, &bs, NULL, &field_array_void); CHKERRQ(ierr);
 
    // Call your existing writer with the temporary, populated Vec
    ierr = WriteFieldData(user, field_name, temp_vec, ti, ext); CHKERRQ(ierr);
 
    // Clean up
    ierr = VecDestroy(&temp_vec); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

WriteAllSwarmFields()

PetscErrorCode WriteAllSwarmFields

(

UserCtx *

user

)

Writes a predefined set of PETSc Swarm fields to files.

This function iterates through a hardcoded list of common swarm fields (position, velocity, etc.) and calls the WriteSwarmField() helper function for each one. This provides a straightforward way to output essential particle data at a given simulation step.

This function will only execute if particles are enabled in the simulation (i.e., user->simCtx->np > 0 and user->swarm is not NULL).

Parameters

[in]

user

Pointer to the UserCtx structure containing the simulation context and the PetscSwarm.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Writes a predefined set of PETSc Swarm fields to files.

Local to this translation unit.

Definition at line 1810 of file io.c.

{
    PetscErrorCode ierr;
    SimCtx         *simCtx = user->simCtx;
    
    PetscFunctionBeginUser;
 
    // If no swarm is configured or there are no particles, do nothing and return.
    if (!user->swarm || simCtx->np <= 0) {
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to write swarm fields.\n");
 
    // Ensure the current IO directory is set
    ierr = PetscSNPrintf(simCtx->_io_context_buffer, sizeof(simCtx->_io_context_buffer),
                         "%s/%s", simCtx->output_dir, simCtx->particle_subdir);CHKERRQ(ierr);
    simCtx->current_io_directory = simCtx->_io_context_buffer;
 
    // Write particle position field
    ierr = WriteSwarmField(user, "position", simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write particle velocity field
    ierr = WriteSwarmField(user, "velocity", simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write particle weight field
    ierr = WriteSwarmField(user, "weight", simCtx->step, "dat"); CHKERRQ(ierr);
    
    // Write custom particle field "Psi"
    ierr = WriteSwarmField(user, "Psi", simCtx->step, "dat"); CHKERRQ(ierr);
    
    // Integer fields require special handling
 
    // Write the background mesh cell ID for each particle
    ierr = WriteSwarmIntField(user, "DMSwarm_CellID", simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write the particle location status (e.g., inside or outside the domain)
    ierr = WriteSwarmIntField(user, "DMSwarm_location_status", simCtx->step, "dat"); CHKERRQ(ierr);
 
    // Write the unique particle ID
    ierr = WriteSwarmIntField(user, "DMSwarm_pid", simCtx->step, "dat"); CHKERRQ(ierr);
 
    simCtx->current_io_directory = NULL;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished writing swarm fields.\n");
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadDataFileToArray()

PetscInt ReadDataFileToArray	(	const char *	filename,
		double **	data_out,
		PetscInt *	Nout,
		MPI_Comm	comm
	)

Reads a simple ASCII data file containing one numeric value per line.

This helper performs rank-0 file I/O, broadcasts the parsed result to the rest of the communicator, and returns a replicated array on every rank.

Parameters

filename	Path to the input data file.
data_out	Output pointer to the allocated scalar array.
Nout	Output pointer storing the number of values read.
comm	MPI communicator used for the coordinated read/broadcast sequence.

Returns

Integer value produced by ReadDataFileToArray().

Reads a simple ASCII data file containing one numeric value per line.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

ReadDataFileToArray()

Definition at line 2456 of file io.c.

{
    /* STEP 0: Prepare local variables & log function entry */
    PetscMPIInt    rank, size;
    PetscErrorCode ierr;
    FILE  *fp = NULL;
    PetscInt    N   = 0;            /* number of lines/values read on rank 0 */
    double *array = NULL;      /* pointer to local array on each rank */
    PetscInt    fileExistsFlag = 0; /* 0 = doesn't exist, 1 = does exist */
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Start reading from file: %s\n",
              filename);
 
    /* Basic error checking: data_out, Nout must be non-null. */
    if (!filename || !data_out || !Nout) {
        LOG_ALLOW(GLOBAL, LOG_WARNING,
                  "Null pointer argument provided.\n");
        return 1;
    }
 
    /* Determine rank/size for coordinating I/O. */
    MPI_Comm_rank(comm, &rank);
    MPI_Comm_size(comm, &size);
 
    /* STEP 1: On rank 0, check if file can be opened. */
    if (!rank) {
        fp = fopen(filename, "r");
        if (fp) {
            fileExistsFlag = 1;
            fclose(fp);
        }
    }
 
    /* STEP 2: Broadcast file existence to all ranks. */
    // In ReadDataFileToArray:
    ierr = MPI_Bcast(&fileExistsFlag, 1, MPI_INT, 0, comm); CHKERRQ(ierr);
 
    if (!fileExistsFlag) {
        /* If file does not exist, log & return. */
        if (!rank) {
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "File '%s' not found.\n",
                      filename);
        }
        return 2;
    }
 
    /* STEP 3: Rank 0 re-opens and reads the file, counting lines, etc. */
    if (!rank) {
        fp = fopen(filename, "r");
        if (!fp) {
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "File '%s' could not be opened for reading.\n",
                      filename);
            return 3;
        }
 
        /* (3a) Count lines first. */
        {
            char line[256];
            while (fgets(line, sizeof(line), fp)) {
                N++;
            }
        }
 
        LOG_ALLOW(GLOBAL, LOG_DEBUG,
                  "File '%s' has %d lines.\n",
                  filename, N);
 
        /* (3b) Allocate array on rank 0. */
        array = (double*)malloc(N * sizeof(double));
        if (!array) {
            fclose(fp);
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "malloc failed for array.\n");
            return 4;
        }
 
        /* (3c) Rewind & read values into array. */
        rewind(fp);
        {
            PetscInt i = 0;
            char line[256];
            while (fgets(line, sizeof(line), fp)) {
                double val;
                if (sscanf(line, "%lf", &val) == 1) {
                    array[i++] = val;
                }
            }
        }
        fclose(fp);
 
        LOG_ALLOW(GLOBAL, LOG_INFO,
                  "Successfully read %d values from '%s'.\n",
                  N, filename);
    }
 
    /* STEP 4: Broadcast the integer N to all ranks. */
    ierr = MPI_Bcast(&N, 1, MPI_INT, 0, comm); CHKERRQ(ierr);
 
    /* STEP 5: Each rank allocates an array to receive the broadcast if rank>0. */
    if (rank) {
        array = (double*)malloc(N * sizeof(double));
        if (!array) {
            LOG_ALLOW(GLOBAL, LOG_WARNING,
                      "malloc failed on rank %d.\n",
                      rank);
            return 5;
        }
    }
 
    /* STEP 6: Broadcast the actual data from rank 0 to all. */
    ierr = MPI_Bcast(array, N, MPI_DOUBLE, 0, comm); CHKERRQ(ierr);
 
    /* STEP 7: Assign outputs on all ranks. */
    *data_out = array;
    *Nout     = N;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Done. Provided array of length=%d to all ranks.\n",
              N);
    return 0; /* success */
}

CreateVTKFileFromMetadata()

PetscInt CreateVTKFileFromMetadata	(	const char *	filename,
		const VTKMetaData *	meta,
		MPI_Comm	comm
	)

Creates a VTK file from prepared metadata and field payloads.

This helper dispatches to the structured-grid or polydata writer based on the metadata contents and emits the assembled VTK file on the requested communicator.

Parameters

filename	Path to the output VTK file.
meta	VTK metadata describing the output geometry and field payloads.
comm	MPI communicator used by the write operation.

Returns

Integer value produced by CreateVTKFileFromMetadata().

Creates a VTK file from prepared metadata and field payloads.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

CreateVTKFileFromMetadata()

Definition at line 153 of file vtk_io.c.

{
    PetscMPIInt rank;
    MPI_Comm_rank(comm, &rank);
    PetscErrorCode ierr = 0;
 
    if (!rank) {
        char tmp_filename[PETSC_MAX_PATH_LEN];
        FILE *fp = NULL;
 
        ierr = PetscSNPrintf(tmp_filename, sizeof(tmp_filename), "%s.tmp", filename);CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 writing combined VTK file '%s'.\n", filename);
        fp = fopen(tmp_filename, "wb");
        if (!fp) {
            LOG_ALLOW(GLOBAL, LOG_ERROR, "fopen failed for %s.\n", tmp_filename);
            return PETSC_ERR_FILE_OPEN;
        }
 
        PetscInt boffset = 0;
 
        ierr = WriteVTKFileHeader(fp, meta, &boffset);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
 
        if (meta->coords) {
            ierr = WriteVTKAppendedBlock(fp, meta->coords, 3 * meta->npoints, sizeof(PetscScalar));
            if (ierr) {
                fclose(fp);
                remove(tmp_filename);
                return ierr;
            }
        }
 
        for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
            const VTKFieldInfo *field = &meta->point_data_fields[i];
            if (field->data) {
                ierr = WriteVTKAppendedBlock(fp, field->data, field->num_components * meta->npoints, sizeof(PetscScalar));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        if (meta->fileType == VTK_POLYDATA) {
            if (meta->connectivity) {
                ierr = WriteVTKAppendedBlock(fp, meta->connectivity, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
            if (meta->offsets) {
                ierr = WriteVTKAppendedBlock(fp, meta->offsets, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        ierr = WriteVTKFileFooter(fp, meta);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
        if (fclose(fp) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        if (rename(tmp_filename, filename) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 finished writing VTK file '%s'.\n", filename);
    }
    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

VecToArrayOnRank0()

PetscErrorCode VecToArrayOnRank0	(	Vec	inVec,
		PetscInt *	N,
		double **	arrayOut
	)

Gathers the contents of a distributed PETSc Vec into a single array on rank 0.

Parameters

[in]	inVec	The input (possibly distributed) Vec.
[out]	N	The global size of the vector.
[out]	arrayOut	On rank 0, points to the newly allocated array holding all data. On other ranks, it is set to NULL.

Returns

PetscErrorCode Return 0 on success, nonzero on failure.

Gathers the contents of a distributed PETSc Vec into a single array on rank 0.

Local to this translation unit.

Definition at line 1863 of file io.c.

{
    PetscErrorCode ierr;
    MPI_Comm       comm;
    PetscMPIInt    rank;
    PetscInt       globalSize;
    DM             dm = NULL;
    const char    *dmtype = NULL;
 
    /* For DMDA path */
    Vec            nat = NULL, seqNat = NULL;
    VecScatter     scatNat = NULL;
    const PetscScalar *nar = NULL;
    PetscScalar    *buf = NULL;
 
    /* For generic (no DM) path */
    Vec            seq = NULL;
    VecScatter     scat = NULL;
    const PetscScalar *sar = NULL;
 
    PetscFunctionBeginUser;
 
    ierr = PetscObjectGetComm((PetscObject)inVec, &comm); CHKERRQ(ierr);
    ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr);
    ierr = VecGetSize(inVec, &globalSize); CHKERRQ(ierr);
    *N = globalSize;
    *arrayOut = NULL;
 
    ierr = VecGetDM(inVec, &dm); CHKERRQ(ierr);
    if (dm) { ierr = DMGetType(dm, &dmtype); CHKERRQ(ierr); }
 
    if (dmtype && !strcmp(dmtype, DMDA)) {
        /* --- DMDA path: go to NATURAL ordering, then gather to rank 0 --- */
        ierr = DMDACreateNaturalVector(dm, &nat); CHKERRQ(ierr);
        ierr = DMDAGlobalToNaturalBegin(dm, inVec, INSERT_VALUES, nat); CHKERRQ(ierr);
        ierr = DMDAGlobalToNaturalEnd  (dm, inVec, INSERT_VALUES, nat); CHKERRQ(ierr);
 
        ierr = VecScatterCreateToZero(nat, &scatNat, &seqNat); CHKERRQ(ierr);
        ierr = VecScatterBegin(scatNat, nat, seqNat, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = VecScatterEnd  (scatNat, nat, seqNat, INSERT_VALUES, SCATTER_FORWARD);   CHKERRQ(ierr);
 
        if (rank == 0) {
            PetscInt nseq;
            ierr = VecGetLocalSize(seqNat, &nseq); CHKERRQ(ierr);
            ierr = VecGetArrayRead(seqNat, &nar);  CHKERRQ(ierr);
 
            ierr = PetscMalloc1(nseq, &buf);       CHKERRQ(ierr);
            ierr = PetscMemcpy(buf, nar, (size_t)nseq * sizeof(PetscScalar)); CHKERRQ(ierr);
 
            ierr = VecRestoreArrayRead(seqNat, &nar); CHKERRQ(ierr);
            *arrayOut = (double*)buf; /* hand back as double* to match the helper signature */
        }
 
        ierr = VecScatterDestroy(&scatNat); CHKERRQ(ierr);
        ierr = VecDestroy(&seqNat);        CHKERRQ(ierr);
        ierr = VecDestroy(&nat);           CHKERRQ(ierr);
    } else {
        /* --- No DM attached: plain gather in Vec’s global (parallel) layout order --- */
        ierr = VecScatterCreateToZero(inVec, &scat, &seq); CHKERRQ(ierr);
        ierr = VecScatterBegin(scat, inVec, seq, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = VecScatterEnd  (scat, inVec, seq, INSERT_VALUES, SCATTER_FORWARD);   CHKERRQ(ierr);
 
        if (rank == 0) {
            PetscInt nseq;
            ierr = VecGetLocalSize(seq, &nseq); CHKERRQ(ierr);
            ierr = VecGetArrayRead(seq, &sar);  CHKERRQ(ierr);
 
            ierr = PetscMalloc1(nseq, &buf);    CHKERRQ(ierr);
            ierr = PetscMemcpy(buf, sar, (size_t)nseq * sizeof(PetscScalar)); CHKERRQ(ierr);
 
            ierr = VecRestoreArrayRead(seq, &sar); CHKERRQ(ierr);
            *arrayOut = (double*)buf;
        }
 
        ierr = VecScatterDestroy(&scat); CHKERRQ(ierr);
        ierr = VecDestroy(&seq);        CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SwarmFieldToArrayOnRank0()

PetscErrorCode SwarmFieldToArrayOnRank0	(	DM	swarm,
		const char *	field_name,
		PetscInt *	n_total_particles,
		PetscInt *	n_components,
		void **	gathered_array
	)

Gathers any DMSwarm field from all ranks to a single, contiguous array on rank 0.

This is a generic, type-aware version of SwarmFieldToArrayOnRank0. It is a COLLECTIVE operation.

Parameters

[in]	swarm	The DMSwarm to gather from.
[in]	field_name	The name of the field to gather.
[out]	n_total_particles	(Output on rank 0) Total number of particles in the global swarm.
[out]	n_components	(Output on rank 0) Number of components for the field.
[out]	gathered_array	(Output on rank 0) A newly allocated array containing the full, gathered data. The caller is responsible for freeing this memory and for casting it to the correct type.

Returns

PetscErrorCode

Gathers any DMSwarm field from all ranks to a single, contiguous array on rank 0.

Local to this translation unit.

Definition at line 1948 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank, size;
    PetscInt       nlocal, nglobal, bs;
    void           *local_array_void;
    size_t         element_size = 0;
 
    PetscFunctionBeginUser;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
 
    // All ranks get swarm properties to determine send/receive counts
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    ierr = DMSwarmGetSize(swarm, &nglobal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, field_name, &bs, NULL, &local_array_void); CHKERRQ(ierr);
    
    // Determine the size of one element of the field's data type
    if (strcasecmp(field_name, "DMSwarm_pid") == 0) {
        element_size = sizeof(PetscInt64);
    } else if (strcasecmp(field_name, "DMSwarm_CellID") == 0 || strcasecmp(field_name, "DMSwarm_location_status") == 0) {
        element_size = sizeof(PetscInt);
    } else {
        element_size = sizeof(PetscScalar);
    }
 
    if (rank == 0) {
        *n_total_particles = nglobal;
        *n_components = bs;
        *gathered_array = NULL; // Initialize output
    }
 
    if (size == 1) { // Serial case is a simple copy
        if (rank == 0) {
            ierr = PetscMalloc(nglobal * bs * element_size, gathered_array); CHKERRQ(ierr);
            ierr = PetscMemcpy(*gathered_array, local_array_void, nglobal * bs * element_size); CHKERRQ(ierr);
        }
    } else { // Parallel case: use MPI_Gatherv
        PetscInt *recvcounts = NULL, *displs = NULL;
        if (rank == 0) {
            ierr = PetscMalloc1(size, &recvcounts); CHKERRQ(ierr);
            ierr = PetscMalloc1(size, &displs); CHKERRQ(ierr);
        }
        PetscInt sendcount = nlocal * bs;
        
        // Gather the number of elements (not bytes) from each rank
        ierr = MPI_Gather(&sendcount, 1, MPIU_INT, recvcounts, 1, MPIU_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        if (rank == 0) {
            displs[0] = 0;
            // Convert counts and calculate displacements in terms of BYTES
            for (PetscMPIInt i = 0; i < size; i++) recvcounts[i] *= element_size;
            for (PetscMPIInt i = 1; i < size; i++) displs[i] = displs[i-1] + recvcounts[i-1];
            
            ierr = PetscMalloc(nglobal * bs * element_size, gathered_array); CHKERRQ(ierr);
        }
        
        // Use Gatherv with MPI_BYTE to handle any data type generically
        ierr = MPI_Gatherv(local_array_void, nlocal * bs * element_size, MPI_BYTE,
                           *gathered_array, recvcounts, displs, MPI_BYTE,
                           0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
        if (rank == 0) {
            ierr = PetscFree(recvcounts); CHKERRQ(ierr);
            ierr = PetscFree(displs); CHKERRQ(ierr);
        }
    }
 
    ierr = DMSwarmRestoreField(swarm, field_name, &bs, NULL, &local_array_void); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ReadSwarmField()

PetscErrorCode ReadSwarmField	(	UserCtx *	user,
		const char *	field_name,
		PetscInt	ti,
		const char *	ext
	)

Reads data from a file into a specified field of a PETSc DMSwarm.

This function is the counterpart to WriteSwarmField(). It creates a global PETSc vector that references the specified DMSwarm field, uses ReadFieldData() to read the data from a file, and then destroys the global vector reference.

Parameters

[in]	user	Pointer to the UserCtx structure (containing user->swarm).
[in]	field_name	Name of the DMSwarm field to read into (must be previously declared/allocated).
[in]	ti	Time index used to construct the input file name.
[in]	ext	File extension (e.g., "dat" or "bin").

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Note

Compatible with PETSc 3.14.x.

Reads data from a file into a specified field of a PETSc DMSwarm.

Local to this translation unit.

Definition at line 1286 of file io.c.

{
  PetscErrorCode ierr;
  DM             swarm;
  Vec            fieldVec;
 
  PetscFunctionBegin;
 
  swarm = user->swarm;
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG," ReadSwarmField Begins \n");
 
  /* 2) Create a global vector that references the specified Swarm field. */
  ierr = DMSwarmCreateGlobalVectorFromField(swarm, field_name, &fieldVec);CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL,LOG_DEBUG," Vector created from Field \n");
 
  /* 3) Use the ReadFieldData() function to read data into fieldVec. */
  ierr = ReadFieldData(user, field_name, fieldVec, ti, ext);CHKERRQ(ierr);
 
  /* 4) Destroy the global vector reference. */
  ierr = DMSwarmDestroyGlobalVectorFromField(swarm, field_name, &fieldVec);CHKERRQ(ierr);
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadSwarmIntField()

PetscErrorCode ReadSwarmIntField	(	UserCtx *	user,
		const char *	field_name,
		PetscInt	ti,
		const char *	ext
	)

Reads integer swarm data by using ReadFieldData and casting the result.

This function is the counterpart to WriteSwarmIntField. It reads a file containing floating-point data (that was originally integer) into a temporary Vec and then casts it back to the integer swarm field. It works by:

1. Creating a temporary parallel Vec.
2. Calling the standard ReadFieldData() to populate this Vec.
3. Accessing the local data of both the Vec and the swarm field.
4. Populating the swarm's integer field by casting each PetscScalar back to a PetscInt.
5. Destroying the temporary Vec.

Parameters

[in]	user	Pointer to the UserCtx structure.
[in]	field_name	Name of the integer Swarm field to be read.
[in]	ti	Time index for the input file.
[in]	ext	File extension.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads integer swarm data by using ReadFieldData and casting the result.

Local to this translation unit.

Definition at line 1316 of file io.c.

{
    PetscErrorCode ierr;
    DM             swarm = user->swarm;
    Vec            temp_vec;
    PetscInt       nlocal, nglobal, bs, i;
    const PetscScalar *scalar_array; // Read-only pointer from the temp Vec
    void           *field_array_void;
 
 
    PetscFunctionBeginUser;
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Reading '%s' via temporary Vec.\n", field_name);
 
    // Get the properties of the swarm field to determine the expected layout
    ierr = DMSwarmGetLocalSize(swarm, &nlocal); CHKERRQ(ierr);
    ierr = DMSwarmGetSize(swarm, &nglobal); CHKERRQ(ierr);
    // We get the block size but not the data pointer yet
    ierr = DMSwarmGetField(swarm, field_name, &bs, NULL, NULL); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, field_name, &bs, NULL, NULL); CHKERRQ(ierr);
 
    // Create a temporary Vec with the CORRECT layout to receive the data
    ierr = VecCreate(PETSC_COMM_WORLD, &temp_vec); CHKERRQ(ierr);
    ierr = VecSetType(temp_vec, VECMPI); CHKERRQ(ierr);
    ierr = VecSetSizes(temp_vec, nlocal * bs, nglobal * bs); CHKERRQ(ierr);
    ierr = VecSetBlockSize(temp_vec, bs); CHKERRQ(ierr);
    ierr = VecSetUp(temp_vec); CHKERRQ(ierr);
 
    // Call your existing reader to populate the temporary Vec
    ierr = ReadFieldData(user, field_name, temp_vec, ti, ext); CHKERRQ(ierr);
 
    // Get local pointers
    ierr = VecGetArrayRead(temp_vec, &scalar_array); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, field_name, NULL, NULL, &field_array_void); CHKERRQ(ierr);
    
    // Perform the cast back, using the correct loop size (nlocal * bs)
    if (strcmp(field_name, "DMSwarm_pid") == 0) {
        PetscInt64 *int64_array = (PetscInt64 *)field_array_void;
        for (i = 0; i < nlocal * bs; i++) {
            int64_array[i] = (PetscInt64)scalar_array[i];
        }
    } else {
        PetscInt *int_array = (PetscInt *)field_array_void;
        for (i = 0; i < nlocal * bs; i++) {
            int_array[i] = (PetscInt)scalar_array[i];
        }
    }    
 
    // Restore access
    ierr = DMSwarmRestoreField(swarm, field_name, NULL, NULL, &field_array_void); CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(temp_vec, &scalar_array); CHKERRQ(ierr);
 
    // 6. Clean up
    ierr = VecDestroy(&temp_vec); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadAllSwarmFields()

PetscErrorCode ReadAllSwarmFields	(	UserCtx *	user,
		PetscInt	ti
	)

Reads multiple fields (positions, velocity, CellID, and weight) into a DMSwarm.

This function is analogous to ReadSimulationFields() but targets a DMSwarm. Each Swarm field is read from a separate file using ReadSwarmField().

Parameters

[in,out]	user	Pointer to the UserCtx structure containing the DMSwarm (user->swarm).
[in]	ti	Time index for constructing the file name.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Reads multiple fields (positions, velocity, CellID, and weight) into a DMSwarm.

Local to this translation unit.

Definition at line 1378 of file io.c.

{
  PetscErrorCode ierr;
  PetscInt nGlobal;
  SimCtx *simCtx = user->simCtx;
  const char *source_path = NULL;
  const char *particle_ext = "dat";
 
  PetscFunctionBeginUser;
  ierr = DMSwarmGetSize(user->swarm, &nGlobal); CHKERRQ(ierr);
  LOG_ALLOW(GLOBAL, LOG_INFO, "Reading DMSwarm fields for timestep %d (swarm size is %d).\n", ti, nGlobal);
 
  if (nGlobal == 0) {
      LOG_ALLOW(GLOBAL, LOG_INFO, "Swarm is empty for timestep %d. Nothing to read.\n", ti);
      PetscFunctionReturn(0);
  }
 
    // First, determine the top-level source directory based on the execution mode.
    if (simCtx->exec_mode == EXEC_MODE_SOLVER) {
        source_path = simCtx->restart_dir;
    } else if (simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR) {
        source_path = simCtx->pps->source_dir;
        if (simCtx->pps->particleExt[0] != '\0') {
            particle_ext = simCtx->pps->particleExt;
        }
    } else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "Invalid execution mode for reading simulation fields.");
    }
 
    // Set the current I/O directory context
    ierr = PetscSNPrintf(simCtx->_io_context_buffer, sizeof(simCtx->_io_context_buffer),
                         "%s/%s", source_path, simCtx->particle_subdir); CHKERRQ(ierr);
 
    simCtx->current_io_directory = simCtx->_io_context_buffer;
 
  /* 1) Read positions (REQUIRED) */
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Reading mandatory position field...\n");
  ierr = ReadSwarmField(user, "position", ti, particle_ext);
  if (ierr) {
      SETERRQ(PETSC_COMM_WORLD, ierr, "Failed to read MANDATORY 'position' field for step %d. Cannot continue.", ti);
  }
  LOG_ALLOW(GLOBAL, LOG_INFO, "Successfully read mandatory position field for step %d.\n", ti);
 
  /* 2) Read all OPTIONAL fields using the helper function. */
  /* The helper will print a warning and continue if a file is not found. */
  ierr = ReadOptionalSwarmField(user, "velocity",                "Velocity",                   ti, particle_ext); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_pid",             "Particle ID",                ti, particle_ext); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_CellID",          "Cell ID",                    ti, particle_ext); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "weight",                  "Particle Weight",            ti, particle_ext); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "Psi",                     "Scalar Psi",                 ti, particle_ext); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_location_status", "Migration Status",   ti, particle_ext); CHKERRQ(ierr);
 
  simCtx->current_io_directory = NULL; // Clear the I/O context after reading
 
  LOG_ALLOW(GLOBAL, LOG_INFO, "Finished reading DMSwarm fields for timestep %d.\n", ti);
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ReadPositionsFromFile()

PetscErrorCode ReadPositionsFromFile	(	PetscInt	timeIndex,
		UserCtx *	user,
		double **	coordsArray,
		PetscInt *	Ncoords
	)

Reads coordinate data (for particles) from file into a PETSc Vec, then gathers it to rank 0.

This function uses ReadFieldData to fill a PETSc Vec with coordinate data, then leverages VecToArrayOnRank0 to gather that data into a contiguous array (valid on rank 0 only).

Parameters

[in]	timeIndex	The time index used to construct file names.
[in]	user	Pointer to the user context.
[out]	coordsArray	On rank 0, will point to a newly allocated array holding the coordinates.
[out]	Ncoords	On rank 0, the length of coordsArray. On other ranks, 0.

Returns

PetscErrorCode Returns 0 on success, or non-zero on failures.

Reads coordinate data (for particles) from file into a PETSc Vec, then gathers it to rank 0.

Local to this translation unit.

Definition at line 2588 of file io.c.

{
  PetscFunctionBeginUser;
 
  PetscErrorCode ierr;
  Vec            coordsVec;
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating coords Vec.\n");
  ierr = VecCreate(PETSC_COMM_WORLD, &coordsVec);CHKERRQ(ierr);
  ierr = VecSetFromOptions(coordsVec);CHKERRQ(ierr);
 
  // For example: "position" is the name of the coordinate data
  ierr = ReadFieldData(user, "position", coordsVec, timeIndex, "dat");
  if (ierr) {
    LOG_ALLOW(GLOBAL, LOG_ERROR,
              "Error reading position data (ti=%d).\n",
              timeIndex);
    PetscFunctionReturn(ierr);
  }
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "ReadPositions - Gathering coords Vec to rank 0.\n");
  ierr = VecToArrayOnRank0(coordsVec, Ncoords, coordsArray);CHKERRQ(ierr);
 
  ierr = VecDestroy(&coordsVec);CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG,
            "Successfully gathered coordinates. Ncoords=%d.\n", *Ncoords);
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

ReadFieldDataToRank0()

PetscErrorCode ReadFieldDataToRank0	(	PetscInt	timeIndex,
		const char *	fieldName,
		UserCtx *	user,
		double **	scalarArray,
		PetscInt *	Nscalars
	)

Reads a named field from file into a PETSc Vec, then gathers it to rank 0.

This function wraps ReadFieldData and VecToArrayOnRank0 into a single step. The gathered data is stored in scalarArray on rank 0, with its length in Nscalars.

Parameters

[in]	timeIndex	The time index used to construct file names.
[in]	fieldName	Name of the field to be read (e.g., "velocity").
[in]	user	Pointer to the user context.
[out]	scalarArray	On rank 0, a newly allocated array holding the field data.
[out]	Nscalars	On rank 0, length of scalarArray. On other ranks, 0.

Returns

PetscErrorCode Returns 0 on success, or non-zero on failures.

Reads a named field from file into a PETSc Vec, then gathers it to rank 0.

Local to this translation unit.

Definition at line 2626 of file io.c.

{
  PetscFunctionBeginUser;
 
  PetscErrorCode ierr;
  Vec            fieldVec;
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Creating field Vec.\n");
  ierr = VecCreate(PETSC_COMM_WORLD, &fieldVec);CHKERRQ(ierr);
  ierr = VecSetFromOptions(fieldVec);CHKERRQ(ierr);
 
  ierr = ReadFieldData(user, fieldName, fieldVec, timeIndex, "dat");
  if (ierr) {
    LOG_ALLOW(GLOBAL, LOG_ERROR,
              "Error reading field '%s' (ti=%d).\n",
              fieldName, timeIndex);
    PetscFunctionReturn(ierr);
  }
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG, "Gathering field Vec to rank 0.\n");
  ierr = VecToArrayOnRank0(fieldVec, Nscalars, scalarArray);CHKERRQ(ierr);
 
  ierr = VecDestroy(&fieldVec);CHKERRQ(ierr);
 
  LOG_ALLOW(GLOBAL, LOG_DEBUG,
            "Successfully gathered field '%s'. Nscalars=%d.\n",
            fieldName, *Nscalars);
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

DisplayBanner()

PetscErrorCode DisplayBanner

(

SimCtx *

simCtx

)

Displays a structured banner summarizing the simulation configuration.

This function prints key simulation parameters to standard output. It is intended to be called ONLY by MPI rank 0. It retrieves global domain bounds and block metadata from simCtx.

Parameters

[in]

simCtx

Pointer to the master simulation context.

Returns

PetscErrorCode Returns 0 on success.

Displays a structured banner summarizing the simulation configuration.

Local to this translation unit.

Definition at line 2026 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    Cmpnts         global_min_coords, global_max_coords;
    PetscReal      StartTime;
    PetscInt       StartStep,StepsToRun,total_num_particles;
    PetscMPIInt    num_mpi_procs;
 
    // SimCtx *simCtx = user->simCtx;
    UserCtx *user = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    num_mpi_procs = simCtx->size;
    StartTime = simCtx->StartTime;
    StartStep = simCtx->StartStep;
    StepsToRun = simCtx->StepsToRun;
    total_num_particles = simCtx->np;
    BoundingBox *bboxlist_on_rank0 = simCtx->bboxlist;
    
 
    PetscFunctionBeginUser;
 
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "DisplayBanner - UserCtx pointer is NULL.");
    global_min_coords = user->bbox.min_coords;
    global_max_coords = user->bbox.max_coords;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    if (rank == 0) {
        // If global_domain_bbox is not pre-populated in UserCtx, compute it here from bboxlist_on_rank0
        // This assumes bboxlist_on_rank0 is valid and contains all local bounding boxes on rank 0.
        if (bboxlist_on_rank0 && num_mpi_procs > 0) {
            global_min_coords = bboxlist_on_rank0[0].min_coords;
            global_max_coords = bboxlist_on_rank0[0].max_coords;
            for (PetscMPIInt p = 1; p < num_mpi_procs; ++p) {
                global_min_coords.x = PetscMin(global_min_coords.x, bboxlist_on_rank0[p].min_coords.x);
                global_min_coords.y = PetscMin(global_min_coords.y, bboxlist_on_rank0[p].min_coords.y);
                global_min_coords.z = PetscMin(global_min_coords.z, bboxlist_on_rank0[p].min_coords.z);
                global_max_coords.x = PetscMax(global_max_coords.x, bboxlist_on_rank0[p].max_coords.x);
                global_max_coords.y = PetscMax(global_max_coords.y, bboxlist_on_rank0[p].max_coords.y);
                global_max_coords.z = PetscMax(global_max_coords.z, bboxlist_on_rank0[p].max_coords.z);
            }
            // Optionally store this in user->global_domain_bbox if it's useful elsewhere
            // user->global_domain_bbox.min_coords = global_min_coords;
            // user->global_domain_bbox.max_coords = global_max_coords;
        } else {
            // Fallback or warning if bboxlist is not available for global calculation
            LOG_ALLOW(LOCAL, LOG_WARNING, "(Rank 0) - bboxlist not provided or num_mpi_procs <=0; using user->bbox for domain bounds.\n");
        // global_min_coords = user->bbox.min_coords; // Use local bbox of rank 0 as fallback
        // global_max_coords = user->bbox.max_coords;
        }
 
        ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "=============================================================\n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "                          CASE SUMMARY                       \n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "=============================================================\n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Grid Points     : %d X %d X %d\n", user->IM, user->JM, user->KM); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Cells           : %d X %d X %d\n", user->IM - 1, user->JM - 1, user->KM - 1); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Global Domain Bounds (X)    : %.6f to %.6f\n", (double)global_min_coords.x, (double)global_max_coords.x); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Global Domain Bounds (Y)    : %.6f to %.6f\n", (double)global_min_coords.y, (double)global_max_coords.y); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Global Domain Bounds (Z)    : %.6f to %.6f\n", (double)global_min_coords.z, (double)global_max_coords.z); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Periodic Axes (BC-derived)  : I=%s, J=%s, K=%s\n",
                           simCtx->i_periodic ? "YES" : "NO",
                           simCtx->j_periodic ? "YES" : "NO",
                           simCtx->k_periodic ? "YES" : "NO"); CHKERRQ(ierr);
        for (PetscInt axis = 0; axis < 3; axis++) {
            if (!user->periodic_translation_valid[axis]) continue;
            ierr = PetscPrintf(PETSC_COMM_SELF,
                               " Periodic %c Translation     : (%.6e, %.6e, %.6e)\n",
                               "IJK"[axis],
                               (double)user->periodic_translation[axis].x,
                               (double)user->periodic_translation[axis].y,
                               (double)user->periodic_translation[axis].z); CHKERRQ(ierr);
        }
        if (total_num_particles > 0 &&
            (simCtx->i_periodic || simCtx->j_periodic || simCtx->k_periodic)) {
            ierr = PetscPrintf(PETSC_COMM_SELF,
                               " Particle Periodic Wrapping  : UNSUPPORTED (Eulerian periodicity only)\n"); CHKERRQ(ierr);
        }
        if(strcmp(simCtx->eulerianSource,"load")==0 || strcmp(simCtx->eulerianSource,"solve")==0){
            ierr = PetscPrintf(PETSC_COMM_SELF, "-------------------- Boundary Conditions --------------------\n"); CHKERRQ(ierr);
            const int face_name_width = 17; // Adjusted for longer names (Zeta,Eta,Xi)
            for (PetscInt i_face = 0; i_face < 6; ++i_face) {
                BCFace current_face = (BCFace)i_face;
                // The BCFaceToString will now return the Xi, Eta, Zeta versions
                const char* face_str = BCFaceToString(current_face); 
                const char* bc_type_str = BCTypeToString(user->boundary_faces[current_face].mathematical_type);
                const char* bc_handler_type_str = BCHandlerTypeToString(user->boundary_faces[current_face].handler_type);
                if(user->boundary_faces[current_face].mathematical_type == INLET){
                    if(user->boundary_faces[current_face].handler_type == BC_HANDLER_INLET_CONSTANT_VELOCITY){
                        Cmpnts inlet_velocity = {0.0,0.0,0.0};
                        PetscBool found;
                        ierr = GetBCParamReal(user->boundary_faces[current_face].params,"vx",&inlet_velocity.x,&found); CHKERRQ(ierr);
                        ierr = GetBCParamReal(user->boundary_faces[current_face].params,"vy",&inlet_velocity.y,&found); CHKERRQ(ierr);
                        ierr = GetBCParamReal(user->boundary_faces[current_face].params,"vz",&inlet_velocity.z,&found); CHKERRQ(ierr);
                        ierr = PetscPrintf(PETSC_COMM_SELF, " Face %-*s : %s - %s - [%.4f,%.4f,%.4f]\n", 
                        face_name_width, face_str, bc_type_str, bc_handler_type_str,inlet_velocity.x,inlet_velocity.y,inlet_velocity.z); CHKERRQ(ierr);
                        } else if(user->boundary_faces[current_face].handler_type == BC_HANDLER_INLET_PARABOLIC){
                        PetscReal v_max = 0.0;
                        PetscBool found;
                        ierr = GetBCParamReal(user->boundary_faces[current_face].params,"v_max",&v_max,&found); CHKERRQ(ierr);
                        ierr = PetscPrintf(PETSC_COMM_SELF, " Face %-*s : %s - %s - v_max=%.4f\n",
                        face_name_width, face_str, bc_type_str, bc_handler_type_str, v_max); CHKERRQ(ierr);
                        } else if(user->boundary_faces[current_face].handler_type == BC_HANDLER_INLET_PROFILE_FROM_FILE){
                        const char *source_file = "(missing)";
                        for (BC_Param *param = user->boundary_faces[current_face].params; param; param = param->next) {
                            if (strcasecmp(param->key, "source_file") == 0 && param->value) {
                                source_file = param->value;
                                break;
                            }
                        }
                        ierr = PetscPrintf(PETSC_COMM_SELF, " Face %-*s : %s - %s - source_file=%s\n",
                        face_name_width, face_str, bc_type_str, bc_handler_type_str, source_file); CHKERRQ(ierr);
                        }
                    } else if(user->boundary_faces[current_face].handler_type == BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX){
                        PetscReal flux;
                        PetscBool trimflag,foundflux,foundtrimflag;
                        ierr = GetBCParamReal(user->boundary_faces[current_face].params,"target_flux",&flux,&foundflux); CHKERRQ(ierr);
                        ierr = GetBCParamBool(user->boundary_faces[current_face].params,"apply_trim",&trimflag,&foundtrimflag); CHKERRQ(ierr);
                        ierr = PetscPrintf(PETSC_COMM_SELF, " Face %-*s : %s - %s - [%.4f] - %s\n", 
                        face_name_width, face_str, bc_type_str, bc_handler_type_str,flux,trimflag?"Trim Ucont":"No Trim"); CHKERRQ(ierr);
                    } else{    
                        ierr = PetscPrintf(PETSC_COMM_SELF, " Face %-*s : %s - %s\n", 
                                face_name_width, face_str, bc_type_str,bc_handler_type_str); CHKERRQ(ierr);
                    }
            }
        }
        ierr = PetscPrintf(PETSC_COMM_SELF, "-------------------------------------------------------------\n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Run Mode                   : %s\n", simCtx->OnlySetup ? "SETUP ONLY" : "Full Simulation"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Start Time                  : %.4f\n", (double)StartTime); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Timestep Size               : %.4f\n", (double)simCtx->dt); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Starting Step               : %d\n", StartStep); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Total Steps to Run          : %d\n", StepsToRun); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, " Ending Step                : %d\n", StartStep + StepsToRun); CHKERRQ(ierr);
        if (simCtx->tiout > 0) {
            ierr = PetscPrintf(PETSC_COMM_SELF, " Field/Restart Cadence      : every %d step(s)\n", simCtx->tiout); CHKERRQ(ierr);
        } else {
            ierr = PetscPrintf(PETSC_COMM_SELF, " Field/Restart Cadence      : DISABLED\n"); CHKERRQ(ierr);
        }
        ierr = PetscPrintf(PETSC_COMM_SELF, " Immersed Boundary          : %s\n", simCtx->immersed ? "ENABLED" : "DISABLED"); CHKERRQ(ierr);
        if (simCtx->walltimeGuardEnabled) {
            ierr = PetscPrintf(
                PETSC_COMM_SELF,
                " Runtime Walltime Guard     : %s (warmup=%d, multiplier=%.2f, min=%.1f s, alpha=%.2f)\n",
                simCtx->walltimeGuardActive ? "ENABLED" : "CONFIGURED BUT INACTIVE",
                simCtx->walltimeGuardWarmupSteps,
                (double)simCtx->walltimeGuardMultiplier,
                (double)simCtx->walltimeGuardMinSeconds,
                (double)simCtx->walltimeGuardEstimatorAlpha
            ); CHKERRQ(ierr);
        } else {
            ierr = PetscPrintf(PETSC_COMM_SELF, " Runtime Walltime Guard     : DISABLED\n"); CHKERRQ(ierr);
        }
        ierr = PetscPrintf(PETSC_COMM_SELF, " Number of MPI Processes     : %d\n", num_mpi_procs); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD," Number of Particles         : %d\n", total_num_particles); CHKERRQ(ierr);
        if (simCtx->np > 0) {
            const char *particle_init_str = ParticleInitializationToString(simCtx->ParticleInitialization);
 
            if (simCtx->particleConsoleOutputFreq > 0) {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Console Cadence   : every %d step(s)\n", simCtx->particleConsoleOutputFreq); CHKERRQ(ierr);
            } else {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Console Cadence   : DISABLED\n"); CHKERRQ(ierr);
            }
            ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Log Row Sampling  : every %d particle(s)\n", simCtx->LoggingFrequency); CHKERRQ(ierr);
            if (simCtx->StartStep > 0) {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Restart Mode      : %s\n", simCtx->particleRestartMode); CHKERRQ(ierr);
            }
            ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Initialization Mode: %s\n", particle_init_str); CHKERRQ(ierr);
            ierr = PetscPrintf(PETSC_COMM_SELF, " Interpolation Method       : %s\n",
                simCtx->interpolationMethod == INTERP_TRILINEAR ? "Trilinear (direct cell-center)" : "CornerAveraged (legacy)"); CHKERRQ(ierr);
            if (simCtx->ParticleInitialization == PARTICLE_INIT_SURFACE_RANDOM ||
                simCtx->ParticleInitialization == PARTICLE_INIT_SURFACE_EDGES) {
                if (user->inletFaceDefined) {
                    ierr = PetscPrintf(PETSC_COMM_SELF, " Particles Initialized At    : %s (Enum Val: %d)\n", BCFaceToString(user->identifiedInletBCFace), user->identifiedInletBCFace); CHKERRQ(ierr);
                } else {
                    ierr = PetscPrintf(PETSC_COMM_SELF, " Particles Initialized At    : --- (No INLET face identified)\n"); CHKERRQ(ierr);
                }
            }
        }
        if(strcmp(simCtx->eulerianSource,"solve")==0 || strcmp(simCtx->eulerianSource,"load")==0){
            ierr = PetscPrintf(PETSC_COMM_WORLD," Reynolds Number             : %le\n", simCtx->ren); CHKERRQ(ierr);
            //ierr = PetscPrintf(PETSC_COMM_WORLD," Von-Neumann Number          : %le\n", simCtx->vnn); CHKERRQ(ierr);
            if(strcmp(simCtx->eulerianSource,"solve")==0){
                //ierr = PetscPrintf(PETSC_COMM_WORLD," Stanton Number              : %le\n", simCtx->st); CHKERRQ(ierr);
                ierr = PetscPrintf(PETSC_COMM_WORLD," Momentum Equation Solver    : %s\n", MomentumSolverTypeToString(simCtx->mom_solver_type)); CHKERRQ(ierr);
                ierr = PetscPrintf(PETSC_COMM_WORLD," Initial Pseudo-CFL (Courant): %le\n", simCtx->pseudo_cfl); CHKERRQ(ierr);
                ierr = PetscPrintf(PETSC_COMM_WORLD," Large Eddy Simulation Model : %s\n", LESModelToString(simCtx->les)); CHKERRQ(ierr);
            }
            if (strcmp(simCtx->eulerianSource, "load") == 0) {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Eulerian State Source       : load (%s)\n",
                                   simCtx->restart_dir); CHKERRQ(ierr);
            } else if (simCtx->StartStep > 0) {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Eulerian State Source       : restart step %d (%s)\n",
                                   simCtx->StartStep, simCtx->restart_dir); CHKERRQ(ierr);
            } else {
                const char* field_init_str = InitialConditionModeToString(simCtx->initialConditionMode);
                ierr = PetscPrintf(PETSC_COMM_SELF, " Eulerian State Source       : initial condition (%s)\n",
                                   field_init_str); CHKERRQ(ierr);
            }
            if (strcmp(simCtx->eulerianSource, "solve") == 0 && simCtx->StartStep == 0 &&
                simCtx->initialConditionMode == IC_MODE_CONSTANT_CARTESIAN) {
                    ierr = PetscPrintf(PETSC_COMM_SELF,
                        " Constant Velocity (Cart.)   : x=%.4f  y=%.4f  z=%.4f\n",
                        (double)simCtx->InitialConstantContra.x,
                        (double)simCtx->InitialConstantContra.y,
                        (double)simCtx->InitialConstantContra.z); CHKERRQ(ierr);
            } else if (strcmp(simCtx->eulerianSource, "solve") == 0 && simCtx->StartStep == 0 &&
                       simCtx->initialConditionMode == IC_MODE_CONSTANT_STREAMWISE) {
                    ierr = PetscPrintf(PETSC_COMM_SELF,
                        " Constant Velocity (Curv.)   : speed=%.4f  direction=%s\n",
                        (double)simCtx->icVelocityPhysical,
                        FlowDirectionToString(simCtx->flowDirection)); CHKERRQ(ierr);
            } else if (strcmp(simCtx->eulerianSource, "solve") == 0 && simCtx->StartStep == 0 &&
                       simCtx->initialConditionMode == IC_MODE_POISEUILLE) {
                ierr = PetscPrintf(PETSC_COMM_SELF,
                    " Poiseuille Peak Velocity    : speed=%.4f  direction=%s\n",
                    (double)simCtx->icVelocityPhysical,
                    FlowDirectionToString(simCtx->flowDirection)); CHKERRQ(ierr);
            } else if (strcmp(simCtx->eulerianSource, "solve") == 0 && simCtx->StartStep == 0 &&
                       simCtx->initialConditionMode == IC_MODE_FILE) {
                ierr = PetscPrintf(PETSC_COMM_SELF,
                    " Initial Velocity File       : field=%s directory=%s\n",
                    simCtx->initialConditionField == IC_FIELD_UCAT ? "Ucat" : "Ucont",
                    simCtx->initialConditionDirectory); CHKERRQ(ierr);
            }
        } else if(strcmp(simCtx->eulerianSource,"analytical")==0){
            ierr = PetscPrintf(PETSC_COMM_WORLD," Analytical Solution Type : %s\n", simCtx->AnalyticalSolutionType); CHKERRQ(ierr);
        }
        ierr = PetscPrintf(PETSC_COMM_SELF, "=============================================================\n"); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_SELF, "\n"); CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

StringToBCFace()

PetscErrorCode StringToBCFace	(	const char *	str,
		BCFace *	face_out
	)

Converts a face-token string (e.g., "-Xi", "+Eta") to the internal BCFace enum.

Parameters

[in]	str	Input token from configuration.
[out]	face_out	Parsed enum value on success.

Returns

PetscErrorCode 0 on success, non-zero for invalid tokens or null pointers.

Converts a face-token string (e.g., "-Xi", "+Eta") to the internal BCFace enum.

Local to this translation unit.

Definition at line 317 of file io.c.

                                                                 {
    if      (strcasecmp(str, "-Xi")   == 0) *face_out = BC_FACE_NEG_X;
    else if (strcasecmp(str, "+Xi")   == 0) *face_out = BC_FACE_POS_X;
    else if (strcasecmp(str, "-Eta")  == 0) *face_out = BC_FACE_NEG_Y;
    else if (strcasecmp(str, "+Eta")  == 0) *face_out = BC_FACE_POS_Y;
    else if (strcasecmp(str, "-Zeta") == 0) *face_out = BC_FACE_NEG_Z;
    else if (strcasecmp(str, "+Zeta") == 0) *face_out = BC_FACE_POS_Z;
    else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown face specifier: %s", str);
    return 0;
}

Here is the caller graph for this function:

StringToBCType()

PetscErrorCode StringToBCType	(	const char *	str,
		BCType *	type_out
	)

Converts a mathematical BC type string (e.g., "PERIODIC", "WALL") to BCType.

Parameters

[in]	str	Input token from configuration.
[out]	type_out	Parsed enum value on success.

Returns

PetscErrorCode 0 on success, non-zero for invalid tokens or null pointers.

Converts a mathematical BC type string (e.g., "PERIODIC", "WALL") to BCType.

Local to this translation unit.

Definition at line 332 of file io.c.

                                                                 {
    if      (strcasecmp(str, "WALL")      == 0) *type_out = WALL;
    else if (strcasecmp(str, "SYMMETRY")  == 0) *type_out = SYMMETRY;
    else if (strcasecmp(str, "INLET")     == 0) *type_out = INLET;
    else if (strcasecmp(str, "OUTLET")    == 0) *type_out = OUTLET;
    else if (strcasecmp(str, "PERIODIC")  == 0) *type_out = PERIODIC;
    // ... add other BCTypes here ...
    else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown BC Type string: %s", str);
    return 0;
}

Here is the caller graph for this function:

StringToBCHandlerType()

PetscErrorCode StringToBCHandlerType	(	const char *	str,
		BCHandlerType *	handler_out
	)

Converts a BC handler token (implementation strategy) to BCHandlerType.

Parameters

[in]	str	Input handler token from configuration.
[out]	handler_out	Parsed enum value on success.

Returns

PetscErrorCode 0 on success, non-zero for invalid tokens or null pointers.

Converts a BC handler token (implementation strategy) to BCHandlerType.

Local to this translation unit.

Definition at line 347 of file io.c.

                                                                                  {
    if      (strcasecmp(str, "noslip")              == 0) *handler_out = BC_HANDLER_WALL_NOSLIP;
    else if (strcasecmp(str, "constant_velocity")   == 0) *handler_out = BC_HANDLER_INLET_CONSTANT_VELOCITY;
    else if (strcasecmp(str, "conservation")        == 0) *handler_out = BC_HANDLER_OUTLET_CONSERVATION;
    else if (strcasecmp(str, "parabolic")           == 0) *handler_out = BC_HANDLER_INLET_PARABOLIC;
    else if (strcasecmp(str, "prescribed_flow")     == 0) *handler_out = BC_HANDLER_INLET_PROFILE_FROM_FILE;
    else if (strcasecmp(str,"geometric")            == 0) *handler_out = BC_HANDLER_PERIODIC_GEOMETRIC;
    else if (strcasecmp(str,"constant_flux")        == 0) *handler_out = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
    else if (strcasecmp(str,"initial_flux")         == 0) *handler_out = BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX;
    // ... add other BCHandlerTypes here ...
    else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown BC Handler string: %s", str);
    return 0;
}

Here is the caller graph for this function:

ValidateBCHandlerForBCType()

PetscErrorCode ValidateBCHandlerForBCType	(	BCType	type,
		BCHandlerType	handler
	)

Validates that a selected handler is compatible with a mathematical BC type.

Parameters

[in]	type	Mathematical BC type (e.g., WALL, PERIODIC).
[in]	handler	Selected handler implementation enum.

Returns

PetscErrorCode 0 if compatible, non-zero if the combination is invalid.

Validates that a selected handler is compatible with a mathematical BC type.

Local to this translation unit.

Definition at line 365 of file io.c.

                                                                              {
    switch (type) {
        case OUTLET:
        if(handler != BC_HANDLER_OUTLET_CONSERVATION) return PETSC_ERR_ARG_WRONG;
        break;
        case WALL:
            if (handler != BC_HANDLER_WALL_NOSLIP && handler != BC_HANDLER_WALL_MOVING) return PETSC_ERR_ARG_WRONG;
            break;
        case INLET:
            if (handler != BC_HANDLER_INLET_CONSTANT_VELOCITY &&
                handler != BC_HANDLER_INLET_PARABOLIC &&
                handler != BC_HANDLER_INLET_PROFILE_FROM_FILE) return PETSC_ERR_ARG_WRONG;
            break;
        case PERIODIC:
            if(handler != BC_HANDLER_PERIODIC_GEOMETRIC && handler != BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX && handler != BC_HANDLER_PERIODIC_DRIVEN_INITIAL_FLUX) return PETSC_ERR_ARG_WRONG;
        // ... add other validation cases here ...
        default: break;
    }
    return 0; // Combination is valid
}

Here is the caller graph for this function:

FreeBC_ParamList()

void FreeBC_ParamList

(

BC_Param *

head

)

Frees an entire linked list of boundary-condition parameters.

Parameters

[in,out]

head

Head pointer of the BC_Param list to destroy.

Frees an entire linked list of boundary-condition parameters.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

FreeBC_ParamList()

Definition at line 302 of file io.c.

                                      {
    BC_Param *current = head;
    while (current != NULL) {
        BC_Param *next = current->next;
        PetscFree(current->key);
        PetscFree(current->value);
        PetscFree(current);
        current = next;
    }
}

Here is the caller graph for this function:

GetBCParamReal()

PetscErrorCode GetBCParamReal	(	BC_Param *	params,
		const char *	key,
		PetscReal *	value_out,
		PetscBool *	found
	)

Searches a BC_Param linked list for a key and returns its value as a double.

Parameters

	params	The head of the BC_Param linked list.
	key	The key to search for (case-insensitive).
[out]	value_out	The found value, converted to a PetscReal.
[out]	found	Set to PETSC_TRUE if the key was found, PETSC_FALSE otherwise.

Returns

0 on success.

Searches a BC_Param linked list for a key and returns its value as a double.

Local to this translation unit.

Definition at line 390 of file io.c.

                                                                                                         {
    *found = PETSC_FALSE;
    *value_out = 0.0;
    if (!key) return 0; // No key to search for
 
    BC_Param *current = params;
    while (current) {
        if (strcasecmp(current->key, key) == 0) {
            *value_out = atof(current->value);
            *found = PETSC_TRUE;
            return 0; // Found it, we're done
        }
        current = current->next;
    }
    return 0; // It's not an error to not find the key.
}

Here is the caller graph for this function:

GetBCParamBool()

PetscErrorCode GetBCParamBool	(	BC_Param *	params,
		const char *	key,
		PetscBool *	value_out,
		PetscBool *	found
	)

Searches a BC_Param linked list for a key and returns its value as a bool.

Parameters

	params	The head of the BC_Param linked list.
	key	The key to search for (case-insensitive).
[out]	value_out	The found value, converted to a PetscBool.
[out]	found	Set to PETSC_TRUE if the key was found, PETSC_FALSE otherwise.

Returns

0 on success.

Searches a BC_Param linked list for a key and returns its value as a bool.

Local to this translation unit.

Definition at line 411 of file io.c.

                                                                                                         {
    *found = PETSC_FALSE;
    *value_out = PETSC_FALSE;
    if (!key) return 0; // No key to search for
 
    BC_Param *current = params;
    while (current) {
        if (strcasecmp(current->key, key) == 0) {
                        // Key was found.
            *found = PETSC_TRUE;
            
            // Check the value string. Default to FALSE if the value is NULL or doesn't match a "true" string.
            if (current->value && 
               (strcasecmp(current->value, "true") == 0 ||
                strcmp(current->value, "1") == 0         ||
                strcasecmp(current->value, "yes") == 0))
            {
                *value_out = PETSC_TRUE;
            } else {
                *value_out = PETSC_FALSE;
            }
            return 0; // Found it, we're done
        }
        current = current->next;
    }
    return 0; // It's not an error to not find the key.
}

Here is the caller graph for this function:

ParseAllBoundaryConditions()

PetscErrorCode ParseAllBoundaryConditions	(	UserCtx *	user,
		const char *	bcs_input_filename
	)

Parses the boundary conditions file to configure the type, handler, and any associated parameters for all 6 global faces of the domain.

This function performs the following steps:

1. On MPI rank 0, it reads the specified configuration file line-by-line.
2. It parses each line for <Face> <Type> <Handler> [param=value]... format.
3. It validates the parsed strings and stores the configuration, including a linked list of parameters, in a temporary array.
4. It then serializes this configuration and broadcasts it to all other MPI ranks.
5. All ranks (including rank 0) then deserialize the broadcasted data to populate their local user->boundary_faces array identically.
6. It also sets the particle inlet lookup fields in UserCtx.

Parameters

[in,out]	user	The main UserCtx struct where the final configuration for all ranks will be stored.
[in]	bcs_input_filename	The path to the boundary conditions configuration file.

Returns

PetscErrorCode 0 on success, error code on failure.

Parses the boundary conditions file to configure the type, handler, and any associated parameters for all 6 global faces of the domain.

Local to this translation unit.

Definition at line 450 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
 
    // Temporary storage for rank 0 to build the configuration before broadcasting.
    BoundaryFaceConfig configs_rank0[6];
    
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    if (rank == 0) {
        FILE *file;
        char line_buffer[1024];
 
        // Initialize the temporary config array with safe defaults on rank 0.
        for (int i = 0; i < 6; i++) {
            configs_rank0[i].face_id = (BCFace)i;
            configs_rank0[i].mathematical_type = WALL;
            configs_rank0[i].handler_type = BC_HANDLER_WALL_NOSLIP;
            configs_rank0[i].params = NULL;
            configs_rank0[i].handler = NULL; // Handler object is not created here.
        }
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing BC configuration from '%s' on rank 0... \n", bcs_input_filename);
        file = fopen(bcs_input_filename, "r");
        if (!file) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open BCs file '%s'.", bcs_input_filename);
 
        while (fgets(line_buffer, sizeof(line_buffer), file)) {
            char *current_pos = line_buffer;
            while (isspace((unsigned char)*current_pos)) current_pos++; // Skip leading whitespace
            if (*current_pos == '#' || *current_pos == '\0' || *current_pos == '\n' || *current_pos == '\r') continue;
 
            char *face_str = strtok(current_pos, " \t\n\r");
            char *type_str = strtok(NULL, " \t\n\r");
            char *handler_str = strtok(NULL, " \t\n\r");
 
            if (!face_str || !type_str || !handler_str) {
                LOG_ALLOW(GLOBAL, LOG_WARNING, "Malformed line in bcs.dat, skipping: %s \n", line_buffer);
                continue;
            }
 
            BCFace      face_enum;
            BCType      type_enum;
            BCHandlerType handler_enum;
            const char* handler_name_for_log;
 
            // --- Convert strings to enums and validate ---
            ierr = StringToBCFace(face_str, &face_enum); CHKERRQ(ierr);
            ierr = StringToBCType(type_str, &type_enum); CHKERRQ(ierr);
            ierr = StringToBCHandlerType(handler_str, &handler_enum); CHKERRQ(ierr);
            ierr = ValidateBCHandlerForBCType(type_enum, handler_enum);
            if (ierr) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Validation failed: Handler '%s' is not valid for Type '%s' on Face '%s'.\n", handler_str, type_str, face_str);
 
            // Store the core types for the corresponding face
            configs_rank0[face_enum].mathematical_type = type_enum;
            configs_rank0[face_enum].handler_type = handler_enum;
            handler_name_for_log = BCHandlerTypeToString(handler_enum); // Assumes this utility exists
            LOG_ALLOW(GLOBAL, LOG_DEBUG, "  Parsed Face '%s': Type=%s, Handler=%s \n", face_str, type_str, handler_name_for_log);
 
            // --- Parse optional key=value parameters for this face ---
            FreeBC_ParamList(configs_rank0[face_enum].params); // Clear any previous (default) params
            configs_rank0[face_enum].params = NULL;
            BC_Param **param_next_ptr = &configs_rank0[face_enum].params; // Pointer to the 'next' pointer to build the list
 
            char* token;
            while ((token = strtok(NULL, " \t\n\r")) != NULL) {
                char* equals_ptr = strchr(token, '=');
                if (!equals_ptr) {
                    LOG_ALLOW(GLOBAL, LOG_WARNING, "Malformed parameter '%s' on face '%s', skipping. \n", token, face_str);
                    continue;
                }
                
                *equals_ptr = '\0'; // Temporarily split the string at '=' to separate key and value
                char* key_str = token;
                char* value_str = equals_ptr + 1;
 
                BC_Param *new_param;
                ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
                ierr = PetscStrallocpy(key_str, &new_param->key); CHKERRQ(ierr);
                ierr = PetscStrallocpy(value_str, &new_param->value); CHKERRQ(ierr);
                new_param->next = NULL;
 
                *param_next_ptr = new_param;
                param_next_ptr = &new_param->next;
                LOG_ALLOW(GLOBAL, LOG_TRACE, "    - Found param: [%s] = [%s] \n", new_param->key, new_param->value);
            }
        }
        fclose(file);
    }
 
    // =========================================================================
    //               BROADCASTING THE CONFIGURATION FROM RANK 0
    // =========================================================================
    // This is a critical step to ensure all processes have the same configuration.
    
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d broadcasting/receiving BC configuration.\n", rank);
 
    for (int i = 0; i < 6; i++) {
        // --- Broadcast simple enums ---
        if (rank == 0) {
            user->boundary_faces[i] = configs_rank0[i]; // Rank 0 populates its final struct
        }
        ierr = MPI_Bcast(&user->boundary_faces[i].mathematical_type, 1, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
        ierr = MPI_Bcast(&user->boundary_faces[i].handler_type, 1, MPI_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
        
        // --- Serialize and Broadcast the parameter linked list ---
        PetscInt n_params = 0;
        if (rank == 0) { // On rank 0, count the number of parameters to send
            for (BC_Param *p = user->boundary_faces[i].params; p; p = p->next) n_params++;
        }
        ierr = MPI_Bcast(&n_params, 1, MPI_INT, 0, PETSC_COMM_WORLD);CHKERRQ(ierr);
        
        if (rank != 0) { // Non-root ranks need to receive and build the list
            FreeBC_ParamList(user->boundary_faces[i].params); // Ensure list is empty before building
            user->boundary_faces[i].params = NULL;
        }
 
        BC_Param **param_next_ptr = &user->boundary_faces[i].params;
 
        for (int j = 0; j < n_params; j++) {
            char key_buf[256] = {0}, val_buf[256] = {0};
            if (rank == 0) {
                // On rank 0, navigate to the j-th param and copy its data to buffers
                BC_Param *p = user->boundary_faces[i].params;
                for (int k = 0; k < j; k++) p = p->next;
                strncpy(key_buf, p->key, 255);
                strncpy(val_buf, p->value, 255);
            }
 
            ierr = MPI_Bcast(key_buf, 256, MPI_CHAR, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
            ierr = MPI_Bcast(val_buf, 256, MPI_CHAR, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
            if (rank != 0) {
                // On non-root ranks, deserialize: create a new node and append it
                BC_Param *new_param;
                ierr = PetscMalloc1(1, &new_param); CHKERRQ(ierr);
                ierr = PetscStrallocpy(key_buf, &new_param->key); CHKERRQ(ierr);
                ierr = PetscStrallocpy(val_buf, &new_param->value); CHKERRQ(ierr);
                new_param->next = NULL;
                *param_next_ptr = new_param;
                param_next_ptr = &new_param->next;
            } else {
                 // On rank 0, just advance the pointer for the next iteration
                 param_next_ptr = &((*param_next_ptr)->next);
            }
        }
        user->boundary_faces[i].face_id = (BCFace)i; // Ensure face_id is set on all ranks
    }
    
    // --- Set particle inlet lookup fields used by the particle system ---
    user->inletFaceDefined = PETSC_FALSE;
    for (int i=0; i<6; i++) {
        
        if (user->boundary_faces[i].mathematical_type == INLET && !user->inletFaceDefined) {
            user->inletFaceDefined = PETSC_TRUE;
            user->identifiedInletBCFace = (BCFace)i;
            LOG_ALLOW(GLOBAL, LOG_INFO, "Inlet face for particle initialization identified as Face %d.\n", i);
            break; // Found the first one, stop looking
        }
    }
 
    
    if (rank == 0) {
        // Rank 0 can now free the linked lists it created for the temporary storage.
        // As written, user->boundary_faces was populated directly on rank 0, so no extra free is needed.
        // for(int i=0; i<6; i++) FreeBC_ParamList(configs_rank0[i].params); // This would be needed if we used configs_rank0 exclusively
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DeterminePeriodicity()

PetscErrorCode DeterminePeriodicity

(

SimCtx *

simCtx

)

Scans all block-specific boundary condition files to determine a globally consistent periodicity for each dimension, reusing the core type parser.

This is a lightweight pre-parser intended to be called before DMDA creation. It ensures that the periodicity setting is consistent across all blocks, which is a physical requirement for the domain.

1. It collectively verifies that the mandatory BCS file for each block exists.
2. On MPI rank 0, it then iterates through the files.
3. For each line, it attempts to convert the type string to a BCType enum using the standard StringToBCType helper.
4. If the conversion is successful AND the type is PERIODIC, it flags the corresponding face.
5. If the conversion fails (e.g., for "WALL", "INLET", etc.), the error is cleared and the line is simply ignored, as it's not relevant to periodicity.
6. It validates consistency (e.g., -Xi and +Xi match) and ensures all block files specify the same global periodicity.
7. It broadcasts the final three flags (as integers 0 or 1) to all MPI ranks.
8. All ranks update the i_periodic, j_periodic, and k_periodic fields in their SimCtx.

Parameters

[in,out]

simCtx

The master SimCtx struct, containing the bcs_files list and where the final periodicity flags will be stored.

Returns

PetscErrorCode 0 on success, error code on failure.

Scans all block-specific boundary condition files to determine a globally consistent periodicity for each dimension, reusing the core type parser.

Local to this translation unit.

Definition at line 637 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscInt       periodic_flags[3] = {0, 0, 0}; // Index 0:I, 1:J, 2:K
 
    PetscFunctionBeginUser;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
 
    // --- Part 1: Collectively verify all BCS files exist before proceeding ---
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        const char *bcs_filename = simCtx->bcs_files[bi];
        if (!bcs_filename) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_NULL, "BCS filename for block %d is not set in SimCtx.", bi);
        char desc_buf[256];
        PetscBool file_exists;
        snprintf(desc_buf, sizeof(desc_buf), "BCS file for block %d", bi);
        ierr = VerifyPathExistence(bcs_filename, PETSC_FALSE, PETSC_FALSE, desc_buf, &file_exists); CHKERRQ(ierr);
    }
 
    // --- Part 2: Rank 0 does the parsing, since we know all files exist ---
    if (rank == 0) {
        PetscBool global_is_periodic[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};
        PetscBool is_set = PETSC_FALSE;
 
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            const char *bcs_filename = simCtx->bcs_files[bi];
            FILE *file = fopen(bcs_filename, "r");
 
            PetscBool face_is_periodic[6] = {PETSC_FALSE};
            char line_buffer[1024];
 
            while (fgets(line_buffer, sizeof(line_buffer), file)) {
                char *current_pos = line_buffer;
                while (isspace((unsigned char)*current_pos)) current_pos++;
                if (*current_pos == '#' || *current_pos == '\0' || *current_pos == '\n') continue;
 
                // --- Tokenize the line exactly like the main parser ---
                char *face_str = strtok(current_pos, " \t\n\r");
                char *type_str = strtok(NULL, " \t\n\r");
 
                // If the line doesn't have at least two tokens, we can't determine the type.
                if (!face_str || !type_str) continue;
 
                // --- Perform a direct, non-erroring check on the mathematical type string ---
                if (strcasecmp(type_str, "PERIODIC") == 0) {
                    BCFace face_enum;
                    // A malformed face string on a periodic line IS a fatal error.
                    ierr = StringToBCFace(face_str, &face_enum); CHKERRQ(ierr);
                    face_is_periodic[face_enum] = PETSC_TRUE;
                }
                // Any other type_str (e.g., "WALL", "INLET") is correctly and silently ignored.
            }
            fclose(file);
 
            // --- Validate consistency within this file ---
            if (face_is_periodic[BC_FACE_NEG_X] != face_is_periodic[BC_FACE_POS_X])
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Inconsistent X-periodicity in file '%s' for block %d. Both -Xi and +Xi must be periodic or neither.", bcs_filename, bi);
            if (face_is_periodic[BC_FACE_NEG_Y] != face_is_periodic[BC_FACE_POS_Y])
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Inconsistent Y-periodicity in file '%s' for block %d. Both -Eta and +Eta must be periodic or neither.", bcs_filename, bi);
            if (face_is_periodic[BC_FACE_NEG_Z] != face_is_periodic[BC_FACE_POS_Z])
                SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Inconsistent Z-periodicity in file '%s' for block %d. Both -Zeta and +Zeta must be periodic or neither.", bcs_filename, bi);
 
            PetscBool local_is_periodic[3] = {face_is_periodic[BC_FACE_NEG_X], face_is_periodic[BC_FACE_NEG_Y], face_is_periodic[BC_FACE_NEG_Z]};
 
            // --- Validate consistency across block files ---
            if (!is_set) {
                global_is_periodic[0] = local_is_periodic[0];
                global_is_periodic[1] = local_is_periodic[1];
                global_is_periodic[2] = local_is_periodic[2];
                is_set = PETSC_TRUE;
            } else {
                if (global_is_periodic[0] != local_is_periodic[0] || global_is_periodic[1] != local_is_periodic[1] || global_is_periodic[2] != local_is_periodic[2]) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
                            "Periodicity mismatch between blocks. Block 0 requires (I:%d, J:%d, K:%d), but block %d (file '%s') has (I:%d, J:%d, K:%d).",
                            (int)global_is_periodic[0], (int)global_is_periodic[1], (int)global_is_periodic[2],
                            bi, bcs_filename,
                            (int)local_is_periodic[0], (int)local_is_periodic[1], (int)local_is_periodic[2]);
                }
            }
        } // end loop over blocks
 
        periodic_flags[0] = (global_is_periodic[0]) ? 1 : 0;
        periodic_flags[1] = (global_is_periodic[1]) ? 1 : 0;
        periodic_flags[2] = (global_is_periodic[2]) ? 1 : 0;
 
        LOG_ALLOW(GLOBAL, LOG_INFO, "Global periodicity determined: I-periodic=%d, J-periodic=%d, K-periodic=%d\n",
                  periodic_flags[0], periodic_flags[1], periodic_flags[2]);
    }
 
    // --- Part 3: Broadcast the final flags from rank 0 to all other ranks ---
    ierr = MPI_Bcast(periodic_flags, 3, MPIU_INT, 0, PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    // --- All ranks now update their SimCtx ---
    simCtx->i_periodic = periodic_flags[0];
    simCtx->j_periodic = periodic_flags[1];
    simCtx->k_periodic = periodic_flags[2];
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

TrimWhitespace()

void TrimWhitespace

(

char *

str

)

Helper function to trim leading/trailing whitespace from a string.

Parameters

str	The string to trim in-place.

Helper function to trim leading/trailing whitespace from a string.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

TrimWhitespace()

Definition at line 38 of file io.c.

                               {
    if (!str) return;
 
    char *start = str;
    // Find the first non-whitespace character
    while (isspace((unsigned char)*start)) {
        start++;
    }
 
    // Find the end of the string
    char *end = str + strlen(str) - 1;
    // Move backwards from the end to find the last non-whitespace character
    while (end > start && isspace((unsigned char)*end)) {
        end--;
    }
 
    // Null-terminate after the last non-whitespace character
    *(end + 1) = '\0';
 
    // If there was leading whitespace, shift the string to the left
    if (str != start) {
        memmove(str, start, (end - start) + 2); // +2 to include the new null terminator
    }
}

Here is the caller graph for this function:

ParsePostProcessingSettings()

PetscErrorCode ParsePostProcessingSettings

(

SimCtx *

simCtx

)

Initializes post-processing settings from a config file and command-line overrides.

This function establishes the configuration for a post-processing run by:

1. Setting hardcoded default values in the PostProcessParams struct.
2. Reading a configuration file to override the defaults.
3. Parsing command-line options (-startTime, -endTime, etc.) which can override both the defaults and the file settings.

Parameters

simCtx

The pointer to the simulation context that contains the postprocessing file and struct.

Returns

PetscErrorCode

Initializes post-processing settings from a config file and command-line overrides.

Local to this translation unit.

Definition at line 2264 of file io.c.

{
    FILE *file;
    char line[1024];
    PetscBool startTimeSet, endTimeSet, timeStepSet;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (!simCtx || !simCtx->pps) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_NULL, "SimCtx or its pps member is NULL in ParsePostProcessingSettings.");
    }
 
    char *configFile = simCtx->PostprocessingControlFile;
    PostProcessParams *pps = simCtx->pps;
 
 
    // --- 1. Set Sane Defaults First ---
    pps->startTime = 0;
    pps->endTime = 0;
    pps->timeStep = 1;
    pps->outputParticles = PETSC_FALSE;
    pps->particle_output_freq = simCtx->LoggingFrequency; // Default to logging frequency;
    strcpy(pps->process_pipeline, "");
    strcpy(pps->output_fields_instantaneous, "Ucat,P");
    strcpy(pps->output_fields_averaged, "");
    strcpy(pps->output_prefix, "Field");
    strcpy(pps->particle_output_prefix,"Particle");
    strcpy(pps->particle_fields,"velocity,CellID,weight,pid");
    strcpy(pps->particle_pipeline,"");
    strncpy(pps->statistics_pipeline,      "",       MAX_PIPELINE_LENGTH - 1);
    strncpy(pps->statistics_output_prefix, "Stats",  MAX_FILENAME_LENGTH - 1);
    strcpy(pps->particleExt,"dat"); // The input file format for particles.
    strcpy(pps->eulerianExt,"dat"); // The input file format for Eulerian fields.
    pps->reference[0] = pps->reference[1] = pps->reference[2] = 1;
    strncpy(pps->source_dir, simCtx->output_dir, sizeof(pps->source_dir) - 1);
    pps->source_dir[sizeof(pps->source_dir) - 1] = '\0'; // Ensure null-termination
 
    // --- 2. Parse the Configuration File (overrides defaults) ---
    file = fopen(configFile, "r");
    if (file) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Parsing post-processing config file: %s\n", configFile);
        while (fgets(line, sizeof(line), file)) {
            char *key, *value, *comment;
            comment = strchr(line, '#'); if (comment) *comment = '\0';
            TrimWhitespace(line); if (strlen(line) == 0) continue;
            key = strtok(line, "="); value = strtok(NULL, "=");
            if (key && value) {
                TrimWhitespace(key); TrimWhitespace(value);
                if (strcmp(key, "startTime") == 0) pps->startTime = atoi(value);
                else if (strcmp(key, "endTime") == 0) pps->endTime = atoi(value);
                else if (strcmp(key, "timeStep") == 0) pps->timeStep = atoi(value);
                else if (strcmp(key, "output_particles") == 0) {
                    if (strcasecmp(value, "true") == 0) pps->outputParticles = PETSC_TRUE;
                }
                else if (strcasecmp(key, "process_pipeline") == 0) {
                    strncpy(pps->process_pipeline, value, MAX_PIPELINE_LENGTH - 1);
                    pps->process_pipeline[MAX_PIPELINE_LENGTH - 1] = '\0'; // Ensure null-termination
                } else if (strcasecmp(key, "output_fields_instantaneous") == 0) {
                    strncpy(pps->output_fields_instantaneous, value, MAX_FIELD_LIST_LENGTH - 1);
                    pps->output_fields_instantaneous[MAX_FIELD_LIST_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "output_fields_averaged") == 0) {
                    strncpy(pps->output_fields_averaged, value, MAX_FIELD_LIST_LENGTH - 1);
                    pps->output_fields_averaged[MAX_FIELD_LIST_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "output_prefix") == 0) {
                    strncpy(pps->output_prefix, value, MAX_FILENAME_LENGTH - 1);
                    pps->output_prefix[MAX_FILENAME_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "particle_output_prefix") == 0) {
                    strncpy(pps->particle_output_prefix, value, MAX_FILENAME_LENGTH - 1);
                    pps->particle_output_prefix[MAX_FILENAME_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "particle_fields_instantaneous") == 0) {
                    strncpy(pps->particle_fields, value, MAX_FIELD_LIST_LENGTH - 1);
                    pps->particle_fields[MAX_FIELD_LIST_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "particle_pipeline") == 0) {
                    strncpy(pps->particle_pipeline, value, MAX_PIPELINE_LENGTH - 1);
                    pps->particle_pipeline[MAX_PIPELINE_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "particle_output_freq") == 0) {
                    pps->particle_output_freq = atoi(value);
                } else if (strcasecmp(key, "statistics_pipeline") == 0) {
                    strncpy(pps->statistics_pipeline, value, MAX_PIPELINE_LENGTH - 1);
                    pps->statistics_pipeline[MAX_PIPELINE_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "statistics_output_prefix") == 0) {
                    strncpy(pps->statistics_output_prefix, value, MAX_FILENAME_LENGTH - 1);
                    pps->statistics_output_prefix[MAX_FILENAME_LENGTH - 1] = '\0';
                } else if (strcasecmp(key, "particleExt") == 0) {
                    strncpy(pps->particleExt, value, sizeof(pps->particleExt) - 1);
                    pps->particleExt[sizeof(pps->particleExt) - 1] = '\0';
                } else if (strcasecmp(key, "eulerianExt") == 0) {
                    strncpy(pps->eulerianExt, value, sizeof(pps->eulerianExt) - 1);
                    pps->eulerianExt[sizeof(pps->eulerianExt) - 1] = '\0';
                } else if (strcmp(key, "reference_ip") == 0) {pps->reference[0] = atoi(value);
                } else if (strcmp(key, "reference_jp") == 0) {pps->reference[1] = atoi(value);
                } else if (strcmp(key, "reference_kp") == 0) {pps->reference[2] = atoi(value);
                } else if (strcasecmp(key, "source_directory") == 0) {
                    strncpy(pps->source_dir, value, sizeof(pps->source_dir) - 1);
                    pps->source_dir[sizeof(pps->source_dir) - 1] = '\0';
                } else {
                    LOG_ALLOW(GLOBAL, LOG_WARNING, "Unknown key '%s' in post-processing config file. Ignoring.\n", key);
                }
                // Add parsing for pipeline, fields, etc. in later phases
            }
        }
        fclose(file);
    } else {
        LOG_ALLOW(GLOBAL, LOG_WARNING, "Could not open post-processing config file '%s'. Using defaults and command-line overrides.\n", configFile);
    }
 
    // --- 3. Parse Command-Line Options (overrides file settings and defaults) ---
    PetscOptionsGetInt(NULL, NULL, "-startTime", &pps->startTime, &startTimeSet);
    PetscOptionsGetInt(NULL, NULL, "-endTime", &pps->endTime, &endTimeSet);
    PetscOptionsGetInt(NULL, NULL, "-timeStep", &pps->timeStep, &timeStepSet);
    PetscOptionsGetBool(NULL, NULL, "-output_particles", &pps->outputParticles, NULL);
 
    if(pps->endTime==-1){
        pps->endTime = simCtx->StartStep + simCtx->StepsToRun; // Total steps if endTime is set to -1.
    }
 
    // If only startTime is given on command line, run for a single step
    if (startTimeSet && !endTimeSet) {
        pps->endTime = pps->startTime;
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Post-processing configured to run from t=%d to t=%d with step %d. Particle output: %s.\n",
              pps->startTime, pps->endTime, pps->timeStep, pps->outputParticles ? "TRUE" : "FALSE");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Process Pipeline: %s\n", pps->process_pipeline);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Instantaneous Output Fields: %s\n", pps->output_fields_instantaneous);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Output Prefix: %s\n", pps->output_prefix);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Output Prefix: %s\n", pps->particle_output_prefix);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Fields: %s\n", pps->particle_fields);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Pipeline: %s\n", pps->particle_pipeline);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Particle Output Frequency: %d\n", pps->particle_output_freq);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Averaged Output Fields (reserved): %s\n", pps->output_fields_averaged);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Post input extensions: Eulerian='.%s', Particle='.%s'\n", pps->eulerianExt, pps->particleExt);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ParseScalingInformation()

PetscErrorCode ParseScalingInformation

(

SimCtx *

simCtx

)

Parses physical scaling parameters from command-line options.

This function reads the reference length, velocity, and density from the PETSc options database (provided via -scaling_L_ref, etc.). It populates the simCtx->scaling struct and calculates the derived reference pressure. It sets default values of 1.0 for a fully non-dimensional case if the options are not provided.

Parameters

[in,out]

simCtx

The simulation context whose 'scaling' member will be populated.

Returns

PetscErrorCode

Parses physical scaling parameters from command-line options.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/io.h.

See also

ParseScalingInformation()

Definition at line 2412 of file io.c.

{
    PetscErrorCode ierr;
    PetscBool      flg;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx is NULL in ParseScalingInformation");
 
    // --- 1. Set default values to 1.0 ---
    // This represents a purely non-dimensional run if no scaling is provided.
    simCtx->scaling.L_ref   = 1.0;
    simCtx->scaling.U_ref   = 1.0;
    simCtx->scaling.rho_ref = 1.0;
    
    // --- 2. Read overrides from the command line / control file ---
    ierr = PetscOptionsGetReal(NULL, NULL, "-scaling_L_ref", &simCtx->scaling.L_ref, &flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-scaling_U_ref", &simCtx->scaling.U_ref, &flg); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL, NULL, "-scaling_rho_ref", &simCtx->scaling.rho_ref, &flg); CHKERRQ(ierr);
    
    // --- 3. Calculate derived scaling factors ---
    // Check for division by zero to be safe, though U_ref should be positive.
    if (simCtx->scaling.U_ref <= 0.0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Reference velocity U_ref must be positive. Got %g", (double)simCtx->scaling.U_ref);
    }
    simCtx->scaling.P_ref = simCtx->scaling.rho_ref * simCtx->scaling.U_ref * simCtx->scaling.U_ref;
 
    // --- 4. Log the final, effective scales for verification ---
    LOG(GLOBAL, LOG_INFO, "---------------- Physical Scales Initialized -----------------\n");
    LOG(GLOBAL, LOG_INFO, "  L_ref: %.4f, U_ref: %.4f, rho_ref: %.4f, P_ref: %.4f\n",
        simCtx->scaling.L_ref, simCtx->scaling.U_ref, simCtx->scaling.rho_ref, simCtx->scaling.P_ref);
    LOG(GLOBAL, LOG_INFO, "--------------------------------------------------------------\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function: