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	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.
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)
PetscInt	CreateVTKFileFromMetadata (const char *filename, const VTKMetaData *meta, MPI_Comm comm)
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 string representation of a face to a BCFace enum.
PetscErrorCode	StringToBCType (const char *str, BCType *type_out)
	Converts a string representation of a BC type to a BCType enum.
PetscErrorCode	StringToBCHandlerType (const char *str, BCHandlerType *handler_out)
	Converts a string representation of a handler to a BCHandlerType enum.
PetscErrorCode	ValidateBCHandlerForBCType (BCType type, BCHandlerType handler)
	Validates that a specific handler is compatible with a general BC type.
void	FreeBC_ParamList (BC_Param *head)
	Frees the memory allocated for a linked list of BC_Param structs.
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.

This function is responsible for reading all per-block array options related to grid geometry, such as dimensions (-im), domain bounds (-xMins, -xMaxs), and stretching ratios (-rxs). It reads the entire array for each option, then uses the block index stored in user->_this to select the correct value and populate the fields of the provided UserCtx struct.

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.

Definition at line 77 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:

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.

This function uses a static-variable pattern to ensure the grid file header is read only once, collectively, by all processes. The first time any process calls this function, it triggers a collective operation where rank 0 reads the file and broadcasts the dimensions for all blocks. This data is stored in static "cached" arrays.

On every call (including the first), the function retrieves the dimensions for the specific block (identified by user->_this) from the cached arrays and populates the IM, JM, and KM fields of the provided UserCtx.

Parameters

user	Pointer to the UserCtx for a specific block.

Returns

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

Definition at line 177 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);
 
            fscanf(fd, "%d\n", &g_nblk_from_file);
            // ---- Read first token; decide whether it is the header or nblk ----
            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) {
                // Header is present – read nblk from the next line
                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);
            } else {
                // No header – the token we just read is actually nblk
                g_nblk_from_file = (PetscInt)strtol(firstTok, NULL, 10);
            }       
            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

Definition at line 762 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:

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]	The	timestep at which the simulation field data needs to be read.

Returns

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

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

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

Parameters

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

Returns

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

Definition at line 1231 of file io.c.

{
    PetscErrorCode ierr;
 
    SimCtx *simCtx = user->simCtx;
    const char *source_path = NULL;
 
    if(simCtx->exec_mode == EXEC_MODE_POSTPROCESSOR){
        source_path = simCtx->pps->source_dir;
    } 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, "dat"); CHKERRQ(ierr);
 
    // Read contravariant velocity field
    ierr = ReadFieldData(user, "vfield", user->Ucont, ti, "dat"); CHKERRQ(ierr);
 
    // Read pressure field
    ierr = ReadFieldData(user, "pfield", user->P, ti, "dat"); CHKERRQ(ierr);
 
    // Read node state field (nvert)
    ierr = ReadFieldData(user, "nvfield", user->Nvert, ti, "dat"); 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, "dat"); CHKERRQ(ierr);
    
    ierr = ReadOptionalField(user, "psifield", "Scalar Psi Field", user->Psi, ti, "dat"); 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.

Definition at line 995 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]	The	timestep at which the simulation field data needs to be read.

Returns

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

Reads statistical fields used for time-averaged simulations.

This function reads data for fields such as Ucat_sum, Ucat_cross_sum, Ucat_square_sum, and P_sum, used for statistical analysis during simulation.

Parameters

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

Returns

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

Definition at line 1311 of file io.c.

{
    PetscErrorCode ierr;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read statistical fields.\n");
 
    ierr = ReadOptionalField(user, "su0", "Velocity Sum",       user->Ucat_sum,        ti, "dat"); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "su1", "Velocity Cross Sum", user->Ucat_cross_sum,  ti, "dat"); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "su2", "Velocity Square Sum",user->Ucat_square_sum, ti, "dat"); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "sp",  "Pressure Sum",       user->P_sum,           ti, "dat"); 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.

This function reads LES-related fields such as Cs (Smagorinsky constant) into their respective PETSc vectors.

Parameters

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

Returns

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

Definition at line 1338 of file io.c.

{
    PetscErrorCode ierr;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read LES fields.\n");
 
    ierr = ReadOptionalField(user, "Nu_t", "Turbulent Viscosity", user->Nu_t, ti, "dat"); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "cs", "Smagorinsky Constant (Cs)", user->CS, ti, "dat"); CHKERRQ(ierr);
 
    DMGlobalToLocalBegin(user->da, user->CS, INSERT_VALUES, user->lCs);
    DMGlobalToLocalEnd(user->da, user->CS, INSERT_VALUES, user->lCs);
    
    DMGlobalToLocalBegin(user->da, user->Nu_t, INSERT_VALUES, user->lNu_t);
    DMGlobalToLocalEnd(user->da, user->Nu_t, INSERT_VALUES, user->lNu_t);
 
    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.

This function reads RANS-related fields such as K_Omega into their respective PETSc vectors.

Parameters

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

Returns

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

Definition at line 1369 of file io.c.

{
    PetscErrorCode ierr;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting to read RANS fields.\n");
 
    ierr = ReadOptionalField(user, "kfield", "K-Omega RANS", user->K_Omega, ti, "dat"); CHKERRQ(ierr);
    ierr = ReadOptionalField(user, "Nu_t", "Turbulent Viscosity", user->Nu_t, ti, "dat"); CHKERRQ(ierr);
 
    VecCopy(user->K_Omega, user->K_Omega_o);
 
    DMGlobalToLocalBegin(user->fda2, user->K_Omega, INSERT_VALUES, user->lK_Omega);
    DMGlobalToLocalEnd(user->fda2, user->K_Omega, INSERT_VALUES, user->lK_Omega);
 
    DMGlobalToLocalBegin(user->da, user->Nu_t, INSERT_VALUES, user->lNu_t);
    DMGlobalToLocalEnd(user->da, user->Nu_t, INSERT_VALUES, user->lNu_t);
 
    DMGlobalToLocalBegin(user->fda2, user->K_Omega_o, INSERT_VALUES, user->lK_Omega_o);
    DMGlobalToLocalEnd(user->fda2, user->K_Omega_o, INSERT_VALUES, user->lK_Omega_o);
 
    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.

This function is now parallel-safe.

• In PARALLEL: All processes send their local data to Rank 0. Rank 0 assembles the data into a temporary sequential vector and writes it to a single file.
• In SERIAL: It performs a direct, simple write.

This ensures the output file is always in a simple, portable format.

Parameters

[in]	user	Pointer to the UserCtx structure.
[in]	field_name	Name of the field (e.g., "position").
[in]	field_vec	The parallel PETSc vector containing the 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.

Dump a distributed PETSc Vec to the single sequential file format used by our restart / post-processing tools.

The companion of ReadFieldData(): it always produces one file (e.g. results/ufield00006_0.dat) regardless of how many MPI ranks are running.

Behaviour

# MPI ranks	Strategy
1	Direct VecView() into the file.
> 1	VecScatterCreateToZero() gathers the distributed Vec onto rank 0.
	Rank 0 writes the sequential Vec; all other ranks allocate no storage.

The routine never alters or destroys the parallel Vec passed in; the gather buffer is created and freed internally.

Parameters

[in]	user	Simulation context (used only for logging).
[in]	field_name	Logical field name (forms part of the filename).
[in]	field_vec	Distributed PETSc Vec to write.
[in]	ti	Timestep index used in the filename.
[in]	ext	File extension ("dat" in our workflow).

Returns

0 on success or a PETSc error code.

Definition at line 1636 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->output_dir){
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Output 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.

Definition at line 1761 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.

Definition at line 1824 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.

Definition at line 1852 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.

Definition at line 1886 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.

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 and newer.

Definition at line 1920 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.

This function provides a bridge to write integer-based swarm fields (like DMSwarm_CellID) using the existing Vec-based I/O routine (WriteFieldData). It works by:

1. Creating a temporary parallel Vec with the same layout as other swarm fields.
2. Accessing the local integer data from the swarm field.
3. Populating the temporary Vec by casting each integer to a PetscScalar.
4. Calling the standard WriteFieldData() function with the temporary Vec.
5. Destroying the temporary Vec.

Parameters

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

Returns

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

Definition at line 1989 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.

Definition at line 2061 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
	)

Definition at line 2643 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
	)

Definition at line 128 of file vtk_io.c.

{
    PetscMPIInt rank;
    MPI_Comm_rank(comm, &rank);
    PetscErrorCode ierr = 0;
 
    if (!rank) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 writing combined VTK file '%s'.\n", filename);
        FILE *fp = fopen(filename, "wb");
        if (!fp) {
            LOG_ALLOW(GLOBAL, LOG_ERROR, "fopen failed for %s.\n", filename);
            return PETSC_ERR_FILE_OPEN;
        }
 
        PetscInt boffset = 0;
        
        ierr = WriteVTKFileHeader(fp, meta, &boffset);
        if(ierr) { fclose(fp); return ierr; }
 
        if (meta->coords) {
            ierr = WriteVTKAppendedBlock(fp, meta->coords, 3 * meta->npoints, sizeof(PetscScalar));
            if(ierr) { fclose(fp); return ierr; }
        }
        /*
        // ======== DEBUG: Dump first few values of Ucat_nodal if present ========
        if (!rank) {
            for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
                const VTKFieldInfo* f = &meta->point_data_fields[i];
                if (strcasecmp(f->name, "Ucat_nodal") == 0) {
                    const PetscScalar *a = f->data;
                    const PetscInt npts = meta->npoints;
                    const PetscInt nshow = (npts < 5) ? npts : 5;
 
                    LOG_ALLOW(GLOBAL, LOG_INFO,
                            "DBG (writer) Ucat_nodal: first %d of %" PetscInt_FMT " tuples:\n",
                            (int)nshow, npts);
                    for (PetscInt t = 0; t < nshow; ++t) {
                        LOG_ALLOW(GLOBAL, LOG_INFO,
                                "  Ucat_nodal[%3" PetscInt_FMT "] = (%g, %g, %g)\n",
                                t, (double)a[3*t+0], (double)a[3*t+1], (double)a[3*t+2]);
                    }
                }
            }
        }
        */
        // ======== END DEBUG ========
               
 
        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); return ierr; }
            }
        }
        if (meta->fileType == VTK_POLYDATA) {
            if (meta->connectivity) {
                ierr = WriteVTKAppendedBlock(fp, meta->connectivity, meta->npoints, sizeof(PetscInt));
                if(ierr) { fclose(fp); return ierr; }
            }
            if (meta->offsets) {
                ierr = WriteVTKAppendedBlock(fp, meta->offsets, meta->npoints, sizeof(PetscInt));
                if(ierr) { fclose(fp); return ierr; }
            }
        }
        ierr = WriteVTKFileFooter(fp, meta);
        if(ierr) { fclose(fp); return ierr; }
        
        fclose(fp);
        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.

If the Vec has a DMDA attached, the gather is performed in DMDA "natural" ordering.

Parameters

[in]	inVec	The PETSc Vec (may be distributed).
[out]	N	Global length of the Vec (includes dof).
[out]	arrayOut	On rank 0, newly allocated buffer with the gathered values (PetscScalar-sized). On other ranks, set to NULL.

Definition at line 2119 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* for drop-in compatibility */
        }
 
        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

Definition at line 2214 of file io.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank, size;
    PetscInt       nlocal, nglobal, bs;
    void           *local_array_void;
    size_t         element_size = 0;
    MPI_Datatype   mpi_type = MPI_BYTE; // We'll send raw bytes
 
    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.

Definition at line 1411 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.

Definition at line 1456 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.

This function reads all necessary and optional fields for a DMSwarm for a given timestep. It assumes the swarm has already been resized to match the particle count in the input files.

The 'position' field is considered MANDATORY. If its file is missing or corrupt, the function will return a fatal error.

All other fields (velocity, CellID, weight, etc.) are OPTIONAL. If their corresponding files are not found, a warning is logged, and the function continues without error. If an optional file IS found but is corrupt or has a size mismatch, it is treated as a fatal error.

Parameters

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

Returns

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

Definition at line 1534 of file io.c.

{
  PetscErrorCode ierr;
  PetscInt nGlobal;
  SimCtx *simCtx = user->simCtx;
  const char *source_path = NULL;
 
  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;
    } 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, "dat");
  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, "dat"); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_pid",             "Particle ID",                ti, "dat"); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_CellID",          "Cell ID",                    ti, "dat"); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "weight",                  "Particle Weight",            ti, "dat"); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "Psi",                     "Scalar Psi",                 ti, "dat"); CHKERRQ(ierr);
  ierr = ReadOptionalSwarmField(user, "DMSwarm_location_status", "Migration Status",   ti, "dat"); 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.

Definition at line 2785 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.

Definition at line 2833 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 from user->global_domain_bbox and boundary conditions for all faces from user->face_bc_types.

Parameters

[in]	user	Pointer to UserCtx structure.
[in]	bboxlist	(If rank 0 needed to compute global_domain_bbox here, otherwise NULL)

Returns

PetscErrorCode Returns 0 on success.

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 from user->global_domain_bbox and boundary conditions for all faces from user->face_bc_types.

Parameters

[in]	user	Pointer to UserCtx structure.
[in]	StartTime	Initial simulation time.
[in]	StartStep	Starting timestep index.
[in]	StepsToRun	Total number of timesteps to run.
[in]	num_mpi_procs	Total number of MPI processes.
[in]	total_num_particles	Total number of particles.
[in]	bboxlist	(If rank 0 needed to compute global_domain_bbox here, otherwise NULL)

Returns

PetscErrorCode Returns 0 on success.

Definition at line 2307 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.");
    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);
        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_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, " 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, " 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(strcmp(simCtx->eulerianSource,"solve")==0 || strcmp(simCtx->eulerianSource,"load")==0){
            const char* field_init_str = FieldInitializationToString(simCtx->FieldInitialization);
            const char* particle_init_str = ParticleInitializationToString(simCtx->ParticleInitialization);
            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);
            ierr = PetscPrintf(PETSC_COMM_SELF, " Particle Initialization Mode: %s\n", particle_init_str); 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    : %le\n", simCtx->pseudo_cfl); CHKERRQ(ierr);
                ierr = PetscPrintf(PETSC_COMM_WORLD," Large Eddy Simulation Model : %s\n", LESModelToString(simCtx->les)); 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);
                }
            }
 
                ierr = PetscPrintf(PETSC_COMM_SELF, " Field Initialization Mode   : %s\n", field_init_str); CHKERRQ(ierr);
            if (simCtx->FieldInitialization == 1) {
                ierr = PetscPrintf(PETSC_COMM_SELF, " Constant Velocity           : 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,"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 string representation of a face to a BCFace enum.

Parameters

	str	The input string (e.g., "-Xi", "+Zeta"). Case-insensitive.
[out]	face_out	The resulting BCFace enum.

Returns

0 on success.

Definition at line 280 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 string representation of a BC type to a BCType enum.

Parameters

	str	The input string (e.g., "WALL", "INLET"). Case-insensitive.
[out]	type_out	The resulting BCType enum.

Returns

0 on success.

Definition at line 297 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 string representation of a handler to a BCHandlerType enum.

Parameters

	str	The input string (e.g., "noslip", "constant_velocity"). Case-insensitive.
[out]	handler_out	The resulting BCHandlerType enum.

Returns

0 on success.

Definition at line 314 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,"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 specific handler is compatible with a general BC type.

Parameters

type	The general BCType.
handler	The specific BCHandlerType.

Returns

0 if compatible, error code otherwise.

Definition at line 333 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) 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 the memory allocated for a linked list of BC_Param structs.

Parameters

head	A pointer to the head of the linked list to be freed.

Definition at line 263 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.

Definition at line 360 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.

Definition at line 385 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 legacy fields in UserCtx for compatibility with other modules.

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.

Definition at line 439 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 legacy fields for compatibility with 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.

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. It parses each line to extract the first (face) and second (type) tokens, mimicking the main BC parser's tokenization logic.
4. It performs a direct, case-insensitive string comparison on the "type" token to check if it is "PERIODIC". Other types are silently ignored.
5. It validates consistency (e.g., -Xi and +Xi match) and ensures all block files specify the same global periodicity.
6. It broadcasts the final three flags (as integers 0 or 1) to all MPI ranks.
7. 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.

Definition at line 645 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.

Parameters

str	The string to be trimmed.

Definition at line 36 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

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

configFile	The path to the configuration file to parse.
pps	Pointer to the PostProcessParams struct to be populated.

Returns

PetscErrorCode

Definition at line 2458 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,"");
    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_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 (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);
 
    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

Definition at line 2593 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: