Momentum source terms added to the contravariant RHS. More...

#include "variables.h"
#include "logging.h"
#include "Metric.h"

Include dependency graph for BodyForces.h:

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

Functions
PetscErrorCode	ComputeDrivenChannelFlowSource (UserCtx *user, Vec Rct)
	Applies a momentum source term to drive flow in a periodic channel or pipe.

PetscErrorCode	LogDrivenFlowDiagnostics (UserCtx *user)
	Appends one row of driven-flow controller state to `driven_flow.csv`.

Detailed Description

Momentum source terms added to the contravariant RHS.

Contract for a body force

Every body force here is invoked from ComputeBodyForces() (src/rhs.c), which is the single extension point. A new force should:

take (UserCtx *user, Vec Rct),
detect for itself whether it is active and return early if not,
accumulate into Rct with +=, never assign.

RULE: DO NOT ADVANCE PER-TIMESTEP STATE HERE

ComputeBodyForces() is called from ComputeRHS(), which runs once per residual evaluation - that is, once per Jameson RK stage under the Picard solver and once per Newton residual evaluation (including every finite difference probe) under Newton-Krylov. It is emphatically NOT called once per physical timestep.

Any force that carries state across calls - a filter, a ramp, a moving average, an integral controller term - must therefore gate its update on simCtx->step and reuse the resolved value for the rest of that step:

if (simCtx->myForceStep != simCtx->step) {
    ... advance the state ...
    simCtx->myForceStep = simCtx->step;
}

Advancing it unconditionally makes the applied force depend on how many residual evaluations preceded it. That is history dependence, and it breaks two things at once: MomentumNewtonKrylov_FormResidual() requires F(X) to be a deterministic function of the trial vector alone, and the Picard shadow-Jacobian estimate assumes body forces are a constant forcing with zero velocity Jacobian.

This is not hypothetical. The driven-flow smoothing EMA in ComputeDrivenChannelFlowSource() had exactly this defect: the applied force walked 0.5, 0.75, 0.875 ... of the way toward its target across evaluations within a single timestep. tests/smoke/run_driven_periodic_regression.sh asserts the force is piecewise constant per step; extend it when adding a stateful force.

Definition in file BodyForces.h.

Function Documentation

◆ ComputeDrivenChannelFlowSource()

PetscErrorCode ComputeDrivenChannelFlowSource	(	UserCtx *	user,
		Vec	Rct
	)

Applies a momentum source term to drive flow in a periodic channel or pipe.

This function is the "engine" of the driven flow control system. It operates by:

Introspecting the boundary condition handlers to see if a DRIVEN_ flow handler is active on any face. This determines if a driven flow is enabled and in which direction ('X', 'Y', or 'Z').
If a driven flow is active, it reads the bulkVelocityCorrection value that was computed by the handler's PreStep method and stored in the SimCtx.
It translates this velocity correction into a momentum source term.
It adds this source term to the appropriate component of the contravariant RHS vector (Rct) for all fluid cells in the domain.

If no driven flow handler is found, this function does nothing.

Parameters

user	The UserCtx containing the simulation state for a single block.
Rct	The PETSc Vec for the contravariant RHS, which will be modified in-place.

Returns: PetscErrorCode 0 on success.

Applies a momentum source term to drive flow in a periodic channel or pipe.

Local to this translation unit.

Definition at line 37 of file BodyForces.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
    PetscFunctionBeginUser;
 
    // --- Step 1: Discover if and where a driven flow is active ---
    const char drivenDirection = DrivenFlowDirection(user); // ' ' when none is configured
 
    // --- Step 2: Early exit if no driven flow is configured ---
    if (drivenDirection == ' ') {
        PetscFunctionReturn(0);
    }
 
    // --- Step 3: Get the control signal and exit if no correction is needed ---
    PetscReal bulkVelocityCorrection = simCtx->bulkVelocityCorrection;
    if (PetscAbsReal(bulkVelocityCorrection) < 1.0e-12) {
        PetscFunctionReturn(0);
    }
 
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Applying driven flow momentum source in '%c' direction.\n",
              simCtx->rank, user->_this, drivenDirection);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "  - Received Bulk Velocity Correction: %le\n", bulkVelocityCorrection);
 
    // --- Step 4: Setup for calculation ---
    DMDALocalInfo info = user->info;
    PetscInt i, j, k;
    PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
    PetscInt lys = (info.ys == 0) ? 1 : info.ys;
    PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
    PetscInt lxe = (info.xs + info.xm == info.mx) ? info.mx - 1 : info.xs + info.xm;
    PetscInt lye = (info.ys + info.ym == info.my) ? info.my - 1 : info.ys + info.ym;
    PetscInt lze = (info.zs + info.zm == info.mz) ? info.mz - 1 : info.zs + info.zm;
 
    Cmpnts ***rct, ***csi, ***eta, ***zet;
    PetscReal ***nvert;
    ierr = DMDAVecGetArray(user->fda, Rct, &rct); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    // Calculate the driving force magnitude for the current timestep, smoothed
    // with the value from the previous step for stability.
    //
    // ONCE PER PHYSICAL STEP, NOT ONCE PER CALL. This function runs from
    // ComputeRHS, which executes once per Jameson RK stage under the Picard
    // solver and once per residual evaluation under Newton-Krylov. The smoothing
    // below carries state in simCtx across calls, so advancing it every call
    // would walk the applied force toward its target within a single timestep
    // (0.5, then 0.75, then 0.875 ... of the way there). The force would then
    // depend on how many residual evaluations preceded it - history dependence
    // that MomentumNewtonKrylov_FormResidual() explicitly forbids, and that
    // breaks the constant-forcing assumption behind the Picard shadow-Jacobian
    // estimate. Resolve it once for the step and reuse it thereafter.
    const PetscReal forceScalingFactor  = simCtx->forceScalingFactor;
    PetscReal drivingForceMagnitude;
 
    if (simCtx->drivingForceStep != simCtx->step) {
        const PetscReal targetForce = (bulkVelocityCorrection / simCtx->dt / 1.0 * COEF_TIME_ACCURACY); // replaced simCtx->st with 1.0.
        drivingForceMagnitude = (simCtx->drivingForceMagnitude * 0.5) + (targetForce * 0.5);
        simCtx->drivingForceMagnitude = drivingForceMagnitude;
        simCtx->drivingForceStep = simCtx->step;
    } else {
        drivingForceMagnitude = simCtx->drivingForceMagnitude;
    }
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "  - Previous driving force:            %le\n", simCtx->drivingForceMagnitude);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "  - New smoothed driving force:        %le\n", drivingForceMagnitude);
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "  - Force scaling factor:              %f\n",  simCtx->forceScalingFactor);
 
    PetscBool hasLoggedApplication = PETSC_FALSE; // Flag to log details only once per rank.
    // --- Step 5: Apply the momentum source to the correct RHS component ---
    for (k = lzs; k < lze; k++) {
        for (j = lys; j < lye; j++) {
            for (i = lxs; i < lxe; i++) {
                if (nvert[k][j][i] < 0.1) { // Apply only to fluid cells
                    PetscReal faceArea = 0.0;
                    PetscReal momentumSource = 0.0;
 
                    switch (drivenDirection) {
                        case 'X':
                            faceArea = sqrt(csi[k][j][i].x * csi[k][j][i].x + csi[k][j][i].y * csi[k][j][i].y + csi[k][j][i].z * csi[k][j][i].z);
                            momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
                            rct[k][j][i].x += momentumSource;
 
                            // Log details for the very first point where force is applied on this rank.
                            if (!hasLoggedApplication) {
                                LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
                                hasLoggedApplication = PETSC_TRUE;
                            }
                                break;
                        case 'Y':
                            faceArea = sqrt(eta[k][j][i].x * eta[k][j][i].x + eta[k][j][i].y * eta[k][j][i].y + eta[k][j][i].z * eta[k][j][i].z);
                            momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
                            rct[k][j][i].y += momentumSource;
                            
                            // Log details for the very first point where force is applied on this rank.
                            if (!hasLoggedApplication) {
                                LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
                                hasLoggedApplication = PETSC_TRUE;
                            }
                            break;
                        case 'Z':
                            faceArea = sqrt(zet[k][j][i].x * zet[k][j][i].x + zet[k][j][i].y * zet[k][j][i].y + zet[k][j][i].z * zet[k][j][i].z);
                            momentumSource = drivingForceMagnitude * forceScalingFactor * faceArea;
                            rct[k][j][i].z += momentumSource;
 
                            // Log details for the very first point where force is applied on this rank.
                            if (!hasLoggedApplication) {
                                LOG_ALLOW(LOCAL, LOG_DEBUG,"Body Force %le added at (%d,%d,%d)\n",momentumSource, k, j, i);
                                hasLoggedApplication = PETSC_TRUE;
                            }
                            break;
                    }
                }
            }
        }
    }
 
    // --- Step 6: Restore arrays ---
    ierr = DMDAVecRestoreArray(user->fda, Rct, &rct); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, (const Cmpnts***)&csi); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, (const Cmpnts***)&eta); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, (const Cmpnts***)&zet); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->da, user->lNvert, (const PetscReal***)&nvert); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ LogDrivenFlowDiagnostics()

PetscErrorCode LogDrivenFlowDiagnostics ( UserCtx * user )

Appends one row of driven-flow controller state to driven_flow.csv.

Active whenever a constant_flux or initial_flux handler drives the flow, in any geometry. Call once per physical step, after the momentum solve, so the row holds the force that step actually applied. Rank 0 writes; block 0 only, since the controller state lives in SimCtx.

Columns: target_flux and measured_flux (the planar-average volumetric flux the correction was computed from this step), cross_section_area, bulk_velocity (their ratio), bulk_velocity_correction, and driving_acceleration, the momentum source per unit mass along the driven axis in solver units: the smoothed magnitude times -driven_flow_scaling_factor, or zero on a step where the correction was negligible and no source was added. In a statistically stationary driven channel its time average is the mean wall stress over the half-height, u_tau^2 / h; in a duct, the mean wall stress times the wetted perimeter over the cross-section area.