Skip to content

Commit

Permalink
Remove filtered velocity until further testing
Browse files Browse the repository at this point in the history
  • Loading branch information
@aprilnovak
aprilnovak committed Oct 12, 2023
1 parent 4b50c0f commit a313c7d
Show file tree
Hide file tree
Showing 3 changed files with 2 additions and 117 deletions.
9 changes: 0 additions & 9 deletions include/base/NekInterface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -684,15 +684,6 @@ struct usrwrkIndices
/// z-velocity of moving boundary (for mesh blending solver)
int mesh_velocity_z;

/// filtered x-velocity of moving boundary
int filtered_velocity_x;

/// filtered y-velocity of moving boundary
int filtered_velocity_y;

/// filtered z-velocity of moving boundary
int filtered_velocity_z;

/// boundary velocity (for separate domain coupling)
int boundary_velocity;

Expand Down
14 changes: 1 addition & 13 deletions include/base/NekRSProblem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,23 +146,14 @@ class NekRSProblem : public NekRSProblemBase
*/
void calculateMeshVelocity(int e, const field::NekWriteEnum & field);

// TODO: Update doxygen entry for velocity integral and prev_disp_update
/**
* Compute the mass flux weighted integral of a given integrand over a set of boundary IDs
* @param[in] boundary_id nekRS boundary IDs for which to perform the integral
* @param[in] integrand field to integrate and weight by mass flux
* @param[in] pp_mesh which NekRS mesh to operate on
* @return mass flux weighted area average of a field
*/
void calculateFilteredVelocity(const std::vector<int> & boundary_id);

/**
* Calculate mesh velocity for NekRS's elasticity solver using current and previous displacement
* values and write it to nrs->usrwrk, from where it can be accessed in nekRS's .oudf file.
* @param[in] e Boundary element that the displacement values belong to
* @param[in] field NekWriteEnum mesh_velocity_x/y/z field
*/
void prev_disp_update(int e, const field::NekWriteEnum & field);

/// Whether a heat source will be applied to NekRS from MOOSE
const bool & _has_heat_source;

Expand Down Expand Up @@ -261,7 +252,4 @@ class NekRSProblem : public NekRSProblemBase

/// volumetric heat source variable read from by nekRS
unsigned int _heat_source_var;

/// whether or not to calculate filtered velocity for FSI problems
const bool & _calc_filtered_velocity;
};
96 changes: 1 addition & 95 deletions src/base/NekRSProblem.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,11 +68,6 @@ NekRSProblem::validParams()
"Whether to conserve the heat flux by individual sideset (as opposed to lumping all sidesets "
"together). Setting this option to true requires syntax changes in the input file to use "
"vector postprocessors, and places restrictions on how the sidesets are set up.");
params.addParam<bool>(
"calculate_filtered_velocity",
false,
"Whether to conserve determine the filtered velocity for FSI calculations. This can be "
"useful if there is significant mesh compression in the solid.");
return params;
}

Expand All @@ -81,8 +76,7 @@ NekRSProblem::NekRSProblem(const InputParameters & params)
_has_heat_source(getParam<bool>("has_heat_source")),
_conserve_flux_by_sideset(getParam<bool>("conserve_flux_by_sideset")),
_abs_tol(getParam<Real>("normalization_abs_tol")),
_rel_tol(getParam<Real>("normalization_rel_tol")),
_calc_filtered_velocity(getParam<bool>("calculate_filtered_velocity"))
_rel_tol(getParam<Real>("normalization_rel_tol"))
{
nekrs::setAbsoluteTol(getParam<Real>("normalization_abs_tol"));
nekrs::setRelativeTol(getParam<Real>("normalization_rel_tol"));
Expand Down Expand Up @@ -118,19 +112,7 @@ NekRSProblem::NekRSProblem(const InputParameters & params)
_usrwrk_indices.push_back("mesh_velocity_x");
_usrwrk_indices.push_back("mesh_velocity_y");
_usrwrk_indices.push_back("mesh_velocity_z");

if (_calc_filtered_velocity)
{
indices.filtered_velocity_x = start++ * nekrs::scalarFieldOffset();
indices.filtered_velocity_y = start++ * nekrs::scalarFieldOffset();
indices.filtered_velocity_z = start++ * nekrs::scalarFieldOffset();
_usrwrk_indices.push_back("filtered_velocity_x");
_usrwrk_indices.push_back("filtered_velocity_y");
_usrwrk_indices.push_back("filtered_velocity_z");
}
}
else
checkUnusedParam(params, "calculate_filtered_velocity", "not using the blending mesh solver");

_minimum_scratch_size_for_coupling = _usrwrk_indices.size() - _first_reserved_usrwrk_slot;

Expand Down Expand Up @@ -373,9 +355,6 @@ NekRSProblem::sendBoundaryDeformationToNek()
if (apply_new_mesh_velocity)
writeBoundarySolution(e, field::mesh_velocity_z, _mesh_velocity_elem);
}

if (_calc_filtered_velocity)
calculateFilteredVelocity(*_boundary);
}
else
{
Expand Down Expand Up @@ -872,79 +851,6 @@ NekRSProblem::calculateMeshVelocity(int e, const field::NekWriteEnum & field)
_nek_mesh->updateDisplacement(e, displacement, disp_field);
}

void
NekRSProblem::calculateFilteredVelocity(const std::vector<int> & boundary_id)
{
// TODO:VERIFY THAT THIS IS WORKING CORRECTLY
mesh_t * mesh = nekrs::flowMesh(); // NOTE: assuming to only use fluid mesh here, this might not
// be entirely correct.
nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();

dfloat * sgeo = nekrs::getSgeo();

double integral = 0.0;

for (int i = 0; i < mesh->Nelements; ++i)
{
for (int j = 0; j < mesh->Nfaces; ++j)
{
int face_id = mesh->EToB[i * mesh->Nfaces + j];

if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
{
int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
for (int v = 0; v < mesh->Nfp; ++v)
{
int vol_id = mesh->vmapM[offset + v];
int surf_offset = mesh->Nsgeo * (offset + v);

double normal_velocity =
nrs->usrwrk[indices.mesh_velocity_x + vol_id] * sgeo[surf_offset + NXID] +
nrs->usrwrk[indices.mesh_velocity_y + vol_id] * sgeo[surf_offset + NYID] +
nrs->usrwrk[indices.mesh_velocity_z + vol_id] * sgeo[surf_offset + NZID];

integral += normal_velocity * sgeo[surf_offset + WSJID];
}
}
}
}

// sum across all processes
double total_integral;
MPI_Allreduce(&integral, &total_integral, 1, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);

// std::cout << "TOTAL VELOCITY INTEGRAL IS " << total_integral << std::endl; //UNCOMMENT FOR
// DEBUGGING

for (int i = 0; i < mesh->Nelements; ++i)
{
for (int j = 0; j < mesh->Nfaces; ++j)
{
int face_id = mesh->EToB[i * mesh->Nfaces + j];

if (std::find(boundary_id.begin(), boundary_id.end(), face_id) != boundary_id.end())
{
int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
for (int v = 0; v < mesh->Nfp; ++v)
{
int vol_id = mesh->vmapM[offset + v];
int surf_offset = mesh->Nsgeo * (offset + v);

nrs->usrwrk[indices.filtered_velocity_x + vol_id] =
nrs->usrwrk[indices.mesh_velocity_x + vol_id] -
total_integral * sgeo[surf_offset + NXID];
nrs->usrwrk[indices.filtered_velocity_y + vol_id] =
nrs->usrwrk[indices.mesh_velocity_y + vol_id] -
total_integral * sgeo[surf_offset + NYID];
nrs->usrwrk[indices.filtered_velocity_z + vol_id] =
nrs->usrwrk[indices.mesh_velocity_z + vol_id] -
total_integral * sgeo[surf_offset + NZID];
}
}
}
}
}

// This function updates the values stores in previous_displacement. This was added to make sure
// that I only update at the start of a timestep rather than at the start of each picard iteration.
void
Expand Down

0 comments on commit a313c7d

Please sign in to comment.