Where is the constraint on the passive scalars applied?

[LuisFLongo]

In :
https://github.com/PrincetonUniversity/athena/wiki/Passive-Scalars

It is stated that:
"Variable floors are applied to these quantities during interface reconstruction and equation of state conversions to ensure that 1) the primitive r remains in [0, 1]"

I was wondering where this enforcement is made, and how can I make sure that this constraint is not applied anymore.

[felker]

athena/src/eos/eos_scalars.cpp

Lines 42 to 55 in a94f18f

	//for (int n=0; n<NSCALARS; ++n) {
	Real& s_n = s(n,k,j,i);
	Real& r_n = r(n,k,j,i);
	// apply passive scalars floor to conserved variable first, then transform:
	// (multi-D fluxes may have caused it to drop below floor)
	s_n = (s_n < scalar_floor_ * d) ? scalar_floor_ * d : s_n;
	r_n = s_n/d;
	// TODO(felker): continue to monitor the acceptability of this absolute 0. floor
	// (may create very large global conservation violations, e.g. the first few
	// cycles of the slotted cylinder test)
	//r_n = (r_n < scalar_floor_) ? scalar_floor_ : r_n;
	//s_n = r_n * d;
	}

The variable floor is mainly applied here. I can't recall at the moment where the max(r) <= 1.0 is explicitly enforced, if at all. It might only be an implicit enforcement through the reconstruction method being monotonicity-preserving. Hence why you see variable floors applied again to PPM-reconstructed states, since that reconstruction method has weaker guarantees compared to the second-order scheme:

athena/src/scalars/calculate_scalar_fluxes.cpp

Lines 74 to 90 in a94f18f

	for (int k=kl; k<=ku; ++k) {
	for (int j=jl; j<=ju; ++j) {
	// reconstruct L/R states
	if (order == 1) {
	pmb->precon->DonorCellX1(k, j, is-1, ie+1, r, rl_, rr_);
	} else if (order == 2) {
	pmb->precon->PiecewiseLinearX1(k, j, is-1, ie+1, r, rl_, rr_);
	} else {
	pmb->precon->PiecewiseParabolicX1(k, j, is-1, ie+1, r, rl_, rr_);
	for (int n=0; n<NSCALARS; ++n) {
	#pragma omp simd
	for (int i=is; i<=ie+1; ++i) {
	pmb->peos->ApplyPassiveScalarFloors(rl_, n, k, j, i);
	pmb->peos->ApplyPassiveScalarFloors(rr_, n, k, j, i);
	}
	}
	}

although I think that second-order reconstruction, despite being monotonicity-preserving in 1D, could possibly lead to r > 1.0 in 2D and 3D in certain cases of strongly converging (possibly grid aligned?) shocks?

[msbc]

At some point, I had convinced myself that this limit is not actually enforced in the code. I have stored thermodynamic quantities in the pasive scalers and it seems to work with positive real numbers.

[felker]

my next step was to ask you and @munan, glad to see your memory matches my recollection