Code Repository
ENV["QT_QPA_PLATFORM"] = "xcb"
using WaterLily
using StaticArrays
using Plots
using CUDA
norm(x) = √sum(abs2, x)
# Define SDF struct for triangle, circles, and one blob
struct TriangleCircleBlobSDF{T}
centers::NTuple{4, SVector{2, T}}
radius::NTuple{4, T}
blob_center::SVector{2, T}
blob_radius::T
amplitudes::SVector{6, T}
end
function (sdf::TriangleCircleBlobSDF)(x, t)
L = 32
r_tri, k = L / 3, sqrt(3.0)
p = x .- SA[3L, L]
x1 = abs(p[1]) - r_tri
y1 = p[2] + r_tri / k
p2 = SA[clamp(x1, -2r_tri, 2r_tri), y1]
if x1 + k * y1 > 0.0
p2 = SA[clamp(x1 - 0.5 * k * y1, -2r_tri, 2r_tri), -k * x1 - y1] ./ 2.0
end
d_min = -norm(p2) * sign(p2[2])
dx = x .- sdf.blob_center
r = norm(dx)
θ = atan(dx[2], dx[1])
r_effective = sdf.blob_radius
for (k, a_k) in enumerate(sdf.amplitudes)
r_effective += a_k * sin(k * θ)
end
d_blob = r - r_effective
d_min = min(d_min, d_blob)
for i in 1:4
d = norm(x .- sdf.centers[i]) - sdf.radius[i]
d_min = min(d_min, d)
end
return d_min
end
function make_sim(; L=2^5, U=1, Re=250, mem=CuArray)
ν = U * L / Re
size = (12L, 4L)
centers = (
SA[6.0L, 1.0L],
SA[2.0L, 2.0L],
SA[3.0L, 3.0L],
SA[2.0L, 3.0L]
)
radius = (L/2.5, L/4, L/4, L/3)
blob_center = SA[8.0L, 2.0L]
blob_radius = L/4 + (L/2) * rand()^0.5
amplitudes = @SVector [0.2 * blob_radius * (2rand() - 1) for _ in 1:6]
sdf = TriangleCircleBlobSDF(centers, radius, blob_center, blob_radius, amplitudes)
body = AutoBody(sdf)
return Simulation(size, (U, 0), L; U=U, ν=ν, body=body, T=Float64, mem=mem)
end
sim = make_sim()
sim_gif!(sim; duration=10, clims=(-10, 10), plotbody=true, dpi=300)
#pragma once
#include "step-69.h"
#include "particle_tracking.h"
#include
#include
#include
#include
namespace Step69
{
template
ParticleTracking::ParticleTracking(
const MPI_Comm mpi_communicator,
TimerOutput &computing_timer,
const OfflineData &offline_data,
const std::string &subsection /*= "ParticleTracking"*/)
: ParameterAcceptor(subsection)
, mpi_communicator(mpi_communicator)
, computing_timer(computing_timer)
, offline_data(&offline_data)
, pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{}
template
void ParticleTracking::generate_particles()
{
const auto &discretization = offline_data->discretization;
const auto &triangulation = discretization->triangulation;
const auto &mapping = discretization->mapping;
particle_handler.initialize(triangulation, mapping, 1 + dim);
Point center;
center[0] = c_x;
center[1] = c_y;
if (dim == 3)
center[2] = 0.0;
const double outer_radius = r_2;
const double inner_radius = r_1;
unsigned int particle_insertion_refinement = 3;
parallel::distributed::Triangulation particle_triangulation(
MPI_COMM_WORLD);
GridGenerator::hyper_shell(
particle_triangulation, center, inner_radius, outer_radius, 6);
particle_triangulation.refine_global(particle_insertion_refinement);
const auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
triangulation, IteratorFilters::LocallyOwnedCell());
const auto global_bounding_boxes =
Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
std::vector> properties(
particle_triangulation.n_locally_owned_active_cells(),
std::vector(dim + 1, 0.));
Particles::Generators::quadrature_points(particle_triangulation,
QMidpoint(),
global_bounding_boxes,
particle_handler,
mapping,
properties);
}
template
void ParticleTracking::euler_step_interpolated(const vector_type &U,
const double dt)
{
const auto get_oriented_normal = [](const auto &face) {
const auto tangential = face->vertex(0) - face->vertex(1);
Tensor<1, dim> normal;
normal[0] = -tangential[1];
normal[1] = tangential[0];
normal /= normal.norm();
return normal;
};
const auto &dof_handler = offline_data->dof_handler;
const auto &finite_element = dof_handler.get_fe();
const double particle_density = dens_graphite;
const double gas_viscosity = mu_H2;
const double diameter = d_0;
double average_velocity = 0;
int counter = 0;
double drag_coeff;
double drag;
Vector local_density_values(finite_element.dofs_per_cell);
std::array, dim> local_velocity_values;
double gas_density = 0.;
for (auto &it : local_velocity_values)
it.reinit(finite_element.dofs_per_cell);
auto particle = particle_handler.begin();
while (particle != particle_handler.end())
{
const auto cell = particle->get_surrounding_cell();
const auto dh_cell =
typename DoFHandler::cell_iterator(*cell, &dof_handler);
dh_cell->get_dof_values(/*density:*/ U[0], local_density_values);
for (int d = 0; d < dim; ++d)
{
dh_cell->get_dof_values(/*momentum component*/ U[1 + d],
local_velocity_values[d]);
for (unsigned int i = 0; i < finite_element.dofs_per_cell; ++i)
local_velocity_values[d][i] /= local_density_values[i];
}
for (unsigned int i = 0; i < finite_element.dofs_per_cell; ++i)
gas_density += local_density_values[i];
const auto pic = particle_handler.particles_in_cell(cell);
Assert(pic.begin() == particle, ExcInternalError());
for (auto &p : pic)
{
const Point reference_location = p.get_reference_location();
Tensor<1, dim> particle_velocity;
Tensor<1, dim> particle_velocity_old;
for (unsigned int j = 0; j < finite_element.dofs_per_cell; ++j)
for (int d = 0; d < dim; ++d)
particle_velocity[d] +=
finite_element.shape_value(j, reference_location) *
local_velocity_values[d][j];
// Drag Force Effect for a Flow
/*********************************/
ArrayView properties_old = p.get_properties();
for (int d = 0; d < dim; ++d) // drag force effect on velocity
particle_velocity_old[d] = 0; //properties_old[d];
double Reynolds =
gas_density *
sqrt(particle_velocity[0] * particle_velocity[0] +
particle_velocity[1] * particle_velocity[1]) *
diameter / gas_viscosity;
drag_coeff = 24.0 / Reynolds +4*pow(Reynolds,0.333); // In Stokes flow
drag = drag_coeff * 3. / 4. * 1 / diameter * gas_density /
particle_density *
sqrt((particle_velocity[0] - particle_velocity_old[0]) *
(particle_velocity[0] - particle_velocity_old[0]) +
(particle_velocity[1] - particle_velocity_old[1]) *
(particle_velocity[1] - particle_velocity_old[1]));
/**************************************************************************/
/* Update particle location: */
Point particle_location = particle->get_location();
for (int d = 0; d < dim; ++d)
{
particle_location[d] += particle_velocity[d] * dt;}
p.set_location(particle_location);
/*CHECKING BOUNDARY CONDITION*/
for (const auto i: dh_cell->face_indices())
{
const auto face = dh_cell->face(i);
if (!face->at_boundary())
continue;
if (face->boundary_id() != Boundaries::free_slip)
continue;
const auto center = face->center();
const auto normal = get_oriented_normal(face);
const double normal_distance =
(particle_location - center) * normal;
if (normal_distance >= 0.)
{
particle_location = particle_location -
(2 * normal_distance) * normal;
p.set_location(particle_location);
}
}
ArrayView properties = p.get_properties();
for (int d = 0; d < dim; ++d) // drag force effect on velocity
properties[d] =
particle_velocity[d] +
dt * drag * (particle_velocity[d] - particle_velocity_old[d]);
/* The drag force and drag heating should interpolate with gas*/
auto drag_force[d] = 1/6*PI*diameter**3*particle_density*
drag*(properties[d]-particle_velocity[d]);
auto drag_heating[d] = 1/6*PI*diameter**3*particle_density*
drag*(properties[d]-particle_velocity[d])**2;
for (int d = 0; d < dim; ++d)
average_velocity += sqrt(properties[0] * properties[0] +
properties[1] * properties[1]);
properties[dim] =
Utilities::MPI::this_mpi_process(mpi_communicator);
++particle;
counter += 1;
}
}
average_velocity = average_velocity / counter;
std::string output_directory = "./";
const std::string file_name("data.dat");
particle_handler.sort_particles_into_subdomains_and_cells();
}
template
void ParticleTracking::output_particles(const unsigned int it)
{
Particles::DataOut particle_output;
Tensor<1, dim> particle_velocity;
std::vector solution_names(dim, "velocity");
solution_names.push_back("process_id");
std::vector
data_component_interpretation(
dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation.push_back(
DataComponentInterpretation::component_is_scalar);
particle_output.build_patches(particle_handler,
solution_names,
data_component_interpretation);
std::string output_directory = "./";
const std::string file_name("interpolated-particles");
pcout << "Writing particle output file: " << file_name << '-' << it
<< std::endl;
particle_output.write_vtu_with_pvtu_record(
output_directory, file_name, it, mpi_communicator, 6);
}
} // namespace Step69