Skip to content

[ConvectionDiffusionApplication] Adding crosswind stabilization to Eulerian Convection–Diffusion element#14272

Open
Marco1410 wants to merge 4 commits intomasterfrom
marco/convection_diffusion
Open

[ConvectionDiffusionApplication] Adding crosswind stabilization to Eulerian Convection–Diffusion element#14272
Marco1410 wants to merge 4 commits intomasterfrom
marco/convection_diffusion

Conversation

@Marco1410
Copy link
Copy Markdown
Contributor

📝 Description
This pull request introduces a crosswind stabilization term in the EulerianConvectionDiffusionElement of Kratos.

The implementation follows the discontinuity-capturing approach proposed in A discontinuity‑capturing crosswind‑dissipation for the finite element solution of the convection‑diffusion equation by Ramon Codina. The method adds a nonlinear diffusion in the direction orthogonal to the flow in order to reduce the spurious oscillations that may still appear when using streamline-based stabilization techniques.

The stabilization term can be activated through the element settings.
By default CROSS_WIND_STABILIZATION_FACTOR is zero.

 "transient_parameters" : {
        "dynamic_tau" : 1.0,
        "theta"       : 0.5,
        "cross_wind_stabilization_factor" : 0.0
   }

Additionally, a unit test has been included to verify the correct behaviour of the implementation.

To illustrate the effect of the stabilization, two example simulations are shown below:

  • Without crosswind stabilization:
    Spurious oscillations appear near steep gradients.
    Without CrossWind

  • With crosswind stabilization:
    The oscillations are significantly reduced.
    WithCrossWind

🆕 Changelog

  • eulerian_conv_diff.cpp
  • convection_diffusion_transient_solver.py
  • test_eulerian_conv_diff.cpp

@Marco1410 Marco1410 changed the title [CONVECTION DIFFUSION] Adding crosswind stabilization to Eulerian Convection–Diffusion element [ConvectionDiffusionApplication] Adding crosswind stabilization to Eulerian Convection–Diffusion element Mar 11, 2026
Copy link
Copy Markdown
Member

@rubenzorrilla rubenzorrilla left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not fully convinced about the elemental cross wind factor calculation. Aside of this, the others are minors.

const double res = -time_derivative -inner_prod(vel_gauss, grad_phi_halfstep);

const double disc_capturing_coeff = 0.5*C*h*fabs(res/norm_grad);
BoundedMatrix<double,TDim,TDim> D = disc_capturing_coeff*( IdentityMatrix(TDim));
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are allocating this for each Gauss pt.

//cross-wind term
if(norm_grad > 1e-3 && norm_vel > 1e-9)
{
const double C = rCurrentProcessInfo[CROSS_WIND_STABILIZATION_FACTOR];
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Access once outside the Gauss pt. loop. Besides avoiding accessing the ProcessInfo many times, you can use it to check if the cross wind stabilization is active and avoid all these operations.

const double time_derivative = Variables.dt_inv*(inner_prod(N,Variables.phi)-inner_prod(N,Variables.phi_old));
const double res = -time_derivative -inner_prod(vel_gauss, grad_phi_halfstep);

const double disc_capturing_coeff = 0.5*C*h*fabs(res/norm_grad);
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const double disc_capturing_coeff = 0.5*C*h*fabs(res/norm_grad);
const double disc_capturing_coeff = 0.5*C*h*std::abs(res/norm_grad);

const double disc_capturing_coeff = 0.5*C*h*fabs(res/norm_grad);
BoundedMatrix<double,TDim,TDim> D = disc_capturing_coeff*( IdentityMatrix(TDim));
const double norm_vel_squared = norm_vel*norm_vel;
D += (std::max( disc_capturing_coeff - tau*norm_vel_squared , 0.0) - disc_capturing_coeff)/(norm_vel_squared) * outer_prod(vel_gauss,vel_gauss);
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure about this? Shouldn't it involve the viscosity besides the tau?

Comment on lines +310 to +311
KRATOS_EXPECT_VECTOR_NEAR(RHS, expected_RHS, 1.0e-4)
KRATOS_EXPECT_MATRIX_NEAR(LHS, expected_LHS, 1.0e-4)
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make it more strict. In order to output the values, you can use std::setprecision (you can find some examples in the FluidDynamicsApplication).

"dynamic_tau": 1.0,
"theta" : 1.0
"theta" : 1.0,
"cross_wind_stabilization_factor": 0.0
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't it automatically add it when validating the defaults?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants