AnalyticalSolution.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscdmcomposite.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "setup.h"
#include "interpolation.h"

Include dependency graph for AnalyticalSolution.h:

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

Go to the source code of this file.

Macros
#define	SetLocalCartesianField(fieldName, fieldValue, coor, FieldInitialization, Constant)
	Generic macro to select the appropriate local Cartesian field setter based on field type.

Functions
PetscErrorCode	SetLocalCartesianField_Scalar (const char *fieldName, PetscReal *scalarField, Cmpnts *coor, PetscInt FieldInitialization, PetscReal Constant)
	Sets the local Cartesian scalar field based on input coordinates.
PetscErrorCode	SetLocalCartesianField_Vector (const char *fieldName, Cmpnts *vecField, Cmpnts *coor, PetscInt FieldInitialization, PetscReal Constant)
	Sets the local Cartesian vector field based on input coordinates.
PetscErrorCode	SetAnalyticalCartesianField (UserCtx *user, const char *fieldName)
	Sets an analytical Cartesian field (scalar or vector) for cell centers based on a field name.
PetscErrorCode	SetAnalyticalContravariantField (UserCtx *user, const char *fieldName)
	Sets an analytical contravariant field (Ucont) at nodes, where each component is evaluated using the physical coordinates of its respective face center.
PetscErrorCode	SetAnalyticalSolution (Vec tempVec, PetscInt FieldInitialization)
	Applies the analytical solution to the position vector.
PetscErrorCode	ApplyAnalyticalBC (UserCtx *user, const char *fieldName)
	Applies analytical boundary conditions to a specified global vector or scalar field.
PetscErrorCode	SetInitialInteriorField (UserCtx *user, const char *fieldName)
	Sets the initial values for the INTERIOR of a specified Eulerian field.

Macro Definition Documentation

SetLocalCartesianField

#define SetLocalCartesianField	(		fieldName,
			fieldValue,
			coor,
			FieldInitialization,
			Constant
	)

Value:

    _Generic((fieldValue), \
    Cmpnts*:    SetLocalCartesianField_Vector, \
    PetscReal*: SetLocalCartesianField_Scalar \
    )(fieldName, fieldValue, coor,FieldInitialization,Constant)

Generic macro to select the appropriate local Cartesian field setter based on field type.

This macro dispatches to:

• SetLocalCartesianField_Vector if the field value pointer is of type Cmpnts*
• SetLocalCartesianField_Scalar if the field value pointer is of type PetscReal*

Parameters

fieldName	Field name for logging.
fieldValue	Pointer to the cell's field value (either Cmpnts* or PetscReal*).
coor	Pointer to the coordinate (Cmpnts*) to be used in the computation.

Definition at line 33 of file AnalyticalSolution.h.

           :    SetLocalCartesianField_Vector, \
    PetscReal*: SetLocalCartesianField_Scalar \
    )(fieldName, fieldValue, coor,FieldInitialization,Constant)

Function Documentation

SetLocalCartesianField_Scalar()

PetscErrorCode SetLocalCartesianField_Scalar	(	const char *	fieldName,
		PetscReal *	scalarField,
		Cmpnts *	coor,
		PetscInt	FieldInitialization,
		PetscReal	Constant
	)

Sets the local Cartesian scalar field based on input coordinates.

This function computes the scalar field value by combining the sine of the input coordinate values. In this example, the scalar field is computed as the sum of the sine functions of the x, y, and z coordinates.

Parameters

[in]	fieldName	Pointer to a string representing the field name (for logging purposes).
[in,out]	scalarField	Pointer to the PetscReal where the computed scalar field value will be stored.
[in]	coor	Pointer to the Cmpnts structure containing the input coordinate values.

Returns

PetscErrorCode Returns 0 on successful execution, non-zero on failure.

SetLocalCartesianField_Vector()

PetscErrorCode SetLocalCartesianField_Vector	(	const char *	fieldName,
		Cmpnts *	vecField,
		Cmpnts *	coor,
		PetscInt	FieldInitialization,
		PetscReal	Constant
	)

Sets the local Cartesian vector field based on input coordinates.

This function computes the vector components by applying the sine function to the input coordinate values along the x, y, and z directions.

Parameters

[in]	fieldName	Pointer to a string representing the field name (for logging purposes).
[in,out]	vecField	Pointer to the Cmpnts structure where the computed vector field will be stored.
[in]	coor	Pointer to the Cmpnts structure containing the input coordinate values.

Returns

PetscErrorCode Returns 0 on successful execution, non-zero on failure.

SetAnalyticalCartesianField()

PetscErrorCode SetAnalyticalCartesianField	(	UserCtx *	user,
		const char *	fieldName
	)

Sets an analytical Cartesian field (scalar or vector) for cell centers based on a field name.

This function looks up the field within the user context by comparing the provided field name against the supported names:

• Vector fields: "Ucat" and "Ucont"
• Scalar fields: "P" and "nvert"

If the field is found, the function verifies that the DM block size matches the expected value (3 for vector fields and 1 for scalar fields), retrieves local DM information for both the coordinate DM (da) and the cell-centered DM (fda), interpolates the corner-based coordinates (from da) to cell centers, and then updates the field using the generic helper macro SetLocalCartesianField. Interior cells are updated with the interpolated coordinates, and boundary cells are updated using the original coordinate data.

If the field name is not found in the user context, the function throws an error.

Parameters

[in]	user	Pointer to the UserCtx structure containing: da: DM for coordinate (corner) data. fda: DM for cell-centered data. Ucat: vector field (Cmpnts ***) Ucont: vector field (Cmpnts ***) P: scalar field (PetscReal ***) nvert: scalar field (PetscReal ***)
[in]	fieldName	Name of the field to update.

Returns

PetscErrorCode 0 on success, non-zero on failure.

SetAnalyticalContravariantField()

PetscErrorCode SetAnalyticalContravariantField	(	UserCtx *	user,
		const char *	fieldName
	)

Sets an analytical contravariant field (Ucont) at nodes, where each component is evaluated using the physical coordinates of its respective face center.

This function:

1. Interpolates nodal physical coordinates to the centers of i-faces, j-faces, and k-faces for the cells relevant to the current rank, storing them in temporary local arrays.
2. For each node (i,j,k) owned by the current rank (where Ucont components are stored): a. Retrieves the pre-calculated physical coordinate of the i-face center. b. Evaluates the analytical function for U^1 using this i-face center coordinate. c. Retrieves the pre-calculated physical coordinate of the j-face center. d. Evaluates the analytical function for U^2 using this j-face center coordinate. e. Retrieves the pre-calculated physical coordinate of the k-face center. f. Evaluates the analytical function for U^3 using this k-face center coordinate. g. Stores these (U^1, U^2, U^3) into user->Ucont[k][j][i].

Parameters

[in]	user	Pointer to the UserCtx structure.
[in]	fieldName	Name of the field to update (currently only "Ucont").

Returns

PetscErrorCode 0 on success, non-zero on failure.

SetAnalyticalSolution()

PetscErrorCode SetAnalyticalSolution	(	Vec	tempVec,
		PetscInt	FieldInitialization
	)

Applies the analytical solution to the position vector.

This function updates each entry in the provided PETSc vector by computing its sine, thereby replacing each position with sin(position).

Parameters

tempVec

The PETSc Vec containing particle positions which will be used to store the velocities.

Returns

PetscErrorCode Returns 0 on success.

ApplyAnalyticalBC()

PetscErrorCode ApplyAnalyticalBC	(	UserCtx *	user,
		const char *	fieldName
	)

Applies analytical boundary conditions to a specified global vector or scalar field.

This function acts as a dispatcher based on the fieldName. It determines the appropriate DM and Vec from the UserCtx and calls field-specific helper routines (e.g., ApplyAnalyticalBC_Vector, ApplyAnalyticalBC_Scalar) to set boundary values according to a predefined analytical function (e.g., sin(coord)).

This should be called AFTER setting interior values and BEFORE updating local ghosts.

Parameters

[in]	user	Pointer to the UserCtx containing DMs, coordinates, and field Vecs.
[in]	fieldName	Name of the field to apply BCs to ("Ucat", "Ucont", "P", "Nvert").

Returns

PetscErrorCode 0 on success, non-zero on failure.

SetInitialInteriorField()

PetscErrorCode SetInitialInteriorField	(	UserCtx *	user,
		const char *	fieldName
	)

Sets the initial values for the INTERIOR of a specified Eulerian field.

This function initializes the interior nodes of Ucont based on a profile selected by user->FieldInitialization. It explicitly skips any node that lies on a global boundary, as those values are set by the Boundary System's Initialize methods.

The initialization is directional, aligned with the primary INLET face that was identified by the parser. This ensures the initial flow is physically meaningful.

Supported user->FieldInitialization profiles for "Ucont":

• 0: Zero Velocity. All interior components of Ucont are set to 0.
• 1: Constant Normal Velocity. The contravariant velocity component normal to the inlet direction is set such that the physical velocity normal to those grid planes is a constant uin. Other contravariant components are zero.
• 2: Poiseuille Normal Velocity. The contravariant component normal to the inlet direction is set with a parabolic profile.

Parameters

user	The main UserCtx struct, containing all simulation data and configuration.
fieldName	A string ("Ucont" or "P") identifying which field to initialize.

Returns

PetscErrorCode 0 on success.

Definition at line 30 of file initialcondition.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    SimCtx *simCtx = user->simCtx;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Setting initial INTERIOR field for '%s' with profile %d.\n", fieldName, simCtx->FieldInitialization);
 
    // This function currently only implements logic for Ucont.
    if (strcmp(fieldName, "Ucont") != 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Skipping SetInitialInteriorField for non-Ucont field '%s'.\n", fieldName);
        PetscFunctionReturn(0);
    }
 
    // --- 1. Get references to required data and PETSc arrays ---
    DM            fieldDM = user->fda;
    Vec           fieldVec = user->Ucont;
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(fieldDM, &info); CHKERRQ(ierr);
 
    Vec      localCoor;
    Cmpnts ***coor_arr;
    ierr = DMGetCoordinatesLocal(user->da, &localCoor); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, localCoor, &coor_arr); CHKERRQ(ierr);
    
    Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
   
    // --- 2. Compute Cell-Center Coordinates (only if needed by the selected profile) ---
    Cmpnts ***cent_coor = NULL;
    PetscInt xs_cell=0, xm_cell=0, ys_cell=0, ym_cell=0, zs_cell=0, zm_cell=0;
    
    if (simCtx->FieldInitialization == 2) { // Profile 2 (Poiseuille) requires cell centers.
      /*
        ierr = GetOwnedCellRange(&info, 0, &xs_cell, &xm_cell); CHKERRQ(ierr);
        ierr = GetOwnedCellRange(&info, 1, &ys_cell, &ym_cell); CHKERRQ(ierr);
        ierr = GetOwnedCellRange(&info, 2, &zs_cell, &zm_cell); CHKERRQ(ierr);
      */
      ierr = GetGhostedCellRange(&info,&user->neighbors,0, &xs_cell, &xm_cell); CHKERRQ(ierr);
      ierr = GetGhostedCellRange(&info,&user->neighbors,1, &ys_cell, &ym_cell); CHKERRQ(ierr);
      ierr = GetGhostedCellRange(&info,&user->neighbors,2, &zs_cell, &zm_cell); CHKERRQ(ierr);
      
        if (xm_cell > 0 && ym_cell > 0 && zm_cell > 0) {
            ierr = Allocate3DArray(&cent_coor, zm_cell, ym_cell, xm_cell); CHKERRQ(ierr);
            ierr = InterpolateFieldFromCornerToCenter(coor_arr, cent_coor, user); CHKERRQ(ierr);
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Computed temporary cell-center coordinates for Poiseuille profile.\n");
        }
    }
    
    // --- 3. Loop Over Owned Nodes and Apply Initial Condition to Interior ---
    Cmpnts ***ucont_arr;
    ierr = DMDAVecGetArray(fieldDM, fieldVec, &ucont_arr); CHKERRQ(ierr);
    
    PetscInt i, j, k;
    const PetscInt mx = info.mx, my = info.my, mz = info.mz; // Global node dimensions
    const PetscInt xs = info.xs, xe = info.xs + info.xm;
    const PetscInt ys = info.ys, ye = info.ys + info.ym;
    const PetscInt zs = info.zs, ze = info.zs + info.zm;
    
    const Cmpnts uin = simCtx->InitialConstantContra; // Max/average velocity from user options
    // Flow into a negative face (e.g., -Zeta at k=0) is in the positive physical direction (+z).
    // Flow into a positive face (e.g., +Zeta at k=mz-1) is in the negative physical direction (-z).
 
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"InitialConstantContra = (%.3f, %.3f, %.3f)\n",(double)uin.x, (double)uin.y, (double)uin.z);
    
    const PetscReal flow_direction_sign = (user->identifiedInletBCFace % 2 == 0) ? 1.0 : -1.0;
        
    for (k = zs; k < ze; k++) {
        for (j = ys; j < ye; j++) {
            for (i = xs; i < xe; i++) {
                
                // The crucial check to ensure we only modify interior nodes.
                const PetscBool is_interior = (i > 0 && i < mx - 1 &&
                                               j > 0 && j < my - 1 &&
                                               k > 0 && k < mz - 1);
 
                if (is_interior) {
                    Cmpnts ucont_val = {0.0, 0.0, 0.0}; // Default to zero velocity
                    PetscReal normal_velocity_mag = 0.0;
 
                    // Step A: Determine the magnitude of the desired physical normal velocity.
                    switch (simCtx->FieldInitialization) {
                        case 0: // Zero initial velocity
                            normal_velocity_mag = 0.0;
                            break;
                case 1: /* Constant Normal Velocity */
                    if (user->identifiedInletBCFace == BC_FACE_NEG_X ||
                user->identifiedInletBCFace == BC_FACE_POS_X) {
                normal_velocity_mag = uin.x;
                    } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
               normal_velocity_mag = uin.y;
                    } else {
               normal_velocity_mag = uin.z;
                    }
              break; 
                        case 2: // Poiseuille Normal Velocity
                            {
                                                  // This profile assumes flow is aligned with the k-index direction.
                                // It uses grid indices (i,j) to define the cross-section, which works for bent geometries.
                                PetscReal u0 = 0.0;
                                if (user->identifiedInletBCFace <= BC_FACE_POS_X) u0 = uin.x;
                                else if (user->identifiedInletBCFace <= BC_FACE_POS_Y) u0 = uin.y;
                                else u0 = uin.z; // Assumes Z-like inlet direction
 
                                // Define channel geometry in "index space" based on global grid dimensions
                                // We subtract 2.0 because the interior runs from index 1 to mx-2 (or my-2).
                                const PetscReal i_width  = (PetscReal)(mx - 2);
                                const PetscReal j_width  = (PetscReal)(my - 2);
                                const PetscReal i_center = 1.0 + i_width / 2.0;
                                const PetscReal j_center = 1.0 + j_width / 2.0;
 
                                // Create normalized coordinates for the current point (i,j), ranging from -1 to 1
                                const PetscReal i_norm = (i - i_center) / (i_width / 2.0);
                                const PetscReal j_norm = (j - j_center) / (j_width / 2.0);
 
                                // Apply the parabolic profile for a rectangular/square channel
                                // V(i,j) = V_max * (1 - i_norm^2) * (1 - j_norm^2)
                                const PetscReal profile_i = 1.0 - i_norm * i_norm;
                                const PetscReal profile_j = 1.0 - j_norm * j_norm;
                                normal_velocity_mag = u0 * profile_i * profile_j;
 
                                // Clamp to zero for any points outside the channel (or due to minor float errors)
                                if (normal_velocity_mag < 0.0) {
                                    normal_velocity_mag = 0.0;
                                }
                            }
                            break;
                /*
                                PetscReal r_sq = 0.0;
                                const PetscInt i_local = i - xs_cell, j_local = j - ys_cell, k_local = k - zs_cell;
                                if (cent_coor && i_local >= 0 && i_local < xm_cell && j_local >= 0 && j_local < ym_cell && k_local >= 0 && k_local < zm_cell) {
                                    const Cmpnts* center = &cent_coor[k_local][j_local][i_local];
                                    if (user->identifiedInletBCFace <= BC_FACE_POS_X) r_sq = center->y * center->y + center->z * center->z;
                                    else if (user->identifiedInletBCFace <= BC_FACE_POS_Y) r_sq = center->x * center->x + center->z * center->z;
                                    else r_sq = center->x * center->x + center->y * center->y;
                    // pick the correct contravariant component for center‐line speed 
                    PetscReal u0;
                    if (user->identifiedInletBCFace == BC_FACE_NEG_X ||
                    user->identifiedInletBCFace == BC_FACE_POS_X) {
                      u0 = uin.x;
                    } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y ||
                           user->identifiedInletBCFace == BC_FACE_POS_Y) {
                      u0 = uin.y;
                    } else {
                      u0 = uin.z;
                    }
                    // now form the parabolic profile as before 
                                    normal_velocity_mag = 2.0 * u0 * (1.0 - 4.0 * r_sq);
                                }
                            }
                            break;
                */
                        default:
                            LOG_ALLOW(LOCAL, LOG_WARNING, "Unrecognized FieldInitialization profile %d. Defaulting to zero.\n", simCtx->FieldInitialization);
                            normal_velocity_mag = 0.0;
                            break;
                    }
 
                    // Step B: Apply direction sign and set the correct contravariant component.
                    // The contravariant component U^n = v_n * Area_n, where v_n is the physical normal velocity.
                    if (normal_velocity_mag != 0.0) {
               const PetscReal signed_normal_vel = normal_velocity_mag * flow_direction_sign*user->GridOrientation;
                        
                        if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
                            const PetscReal area_i = sqrt(csi_arr[k][j][i].x * csi_arr[k][j][i].x + csi_arr[k][j][i].y * csi_arr[k][j][i].y + csi_arr[k][j][i].z * csi_arr[k][j][i].z);
                
                            ucont_val.x = signed_normal_vel * area_i;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_DEBUG,k,50," ucont_val.x = %.6f (signed_normal_vel=%.3f × area=%.4f)\n",ucont_val.x, signed_normal_vel, area_i);
                        } 
                        else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
                            const PetscReal area_j = sqrt(eta_arr[k][j][i].x * eta_arr[k][j][i].x + eta_arr[k][j][i].y * eta_arr[k][j][i].y + eta_arr[k][j][i].z * eta_arr[k][j][i].z);
                 
                            ucont_val.y = signed_normal_vel * area_j;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_DEBUG,k,50," ucont_val.y = %.6f (signed_normal_vel=%.3f × area=%.4f)\n",ucont_val.y, signed_normal_vel, area_j);
                        } 
                        else { // Z-inlet
                            const PetscReal area_k = sqrt(zet_arr[k][j][i].x * zet_arr[k][j][i].x + zet_arr[k][j][i].y * zet_arr[k][j][i].y + zet_arr[k][j][i].z * zet_arr[k][j][i].z);
 
                            ucont_val.z = signed_normal_vel * area_k;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_DEBUG,k,50," i,j,k,ucont_val.z = %d, %d, %d, %.6f (signed_normal_vel=%.3f × area=%.4f)\n",i,j,k,ucont_val.z, signed_normal_vel, area_k);
                        }
                    }
                    ucont_arr[k][j][i] = ucont_val;
                } // end if(is_interior)
            }
        }
    }
    ierr = DMDAVecRestoreArray(fieldDM, fieldVec, &ucont_arr); CHKERRQ(ierr);
 
    // --- 5. Cleanup: Restore arrays and free temporary memory ---
    ierr = DMDAVecRestoreArrayRead(user->fda, localCoor, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
    
    if (cent_coor) {
        ierr = Deallocate3DArray(cent_coor, zm_cell, ym_cell); CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function: