54 using T =
typename CELL::value_t;
55 using DESCRIPTOR =
typename CELL::descriptor_t;
59 const bool drop_isolated_cells = vars.template get<FreeSurface::DROP_ISOLATED_CELLS>();
60 const T transition = vars.template get<FreeSurface::TRANSITION>();
61 const T lonely_threshold = vars.template get<FreeSurface::LONELY_THRESHOLD>();
62 const bool has_surface_tension = vars.template get<FreeSurface::HAS_SURFACE_TENSION>();
63 const T surface_tension_parameter = vars.template get<FreeSurface::SURFACE_TENSION_PARAMETER>();
74 T mass_tmp = cell.template getField<FreeSurface::MASS>();
78 for (
int iPop = 1; iPop < DESCRIPTOR::q; ++iPop){
79 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop, 0),
80 descriptors::c<DESCRIPTOR>(iPop, 1)});
82 int iPop_op = descriptors::opposite<DESCRIPTOR>(iPop);
91 mass_tmp += cell[iPop_op] - cellC[iPop];
100 mass_flow = -cellC[iPop];
102 mass_flow = cell[iPop_op];
104 mass_flow = cell[iPop_op] - cellC[iPop];
107 mass_flow = -cellC[iPop];
112 mass_flow = cell[iPop_op];
114 mass_flow = -cellC[iPop];
116 mass_flow = cell[iPop_op] - cellC[iPop];
119 mass_flow = cell[iPop_op];
123 mass_flow = cell[iPop_op];
125 mass_flow = -cellC[iPop];
127 mass_flow = cell[iPop_op] - cellC[iPop];
134 mass_tmp += mass_flow * 0.5 * (getClampedEpsilon(cell) + getClampedEpsilon(cellC));
138 cell.template setField<FreeSurface::MASS>(mass_tmp);
149 if(has_surface_tension){
152 curvature = calculateSurfaceTensionCurvature(cell);
158 T gas_pressure = 1. - 6. * surface_tension_parameter * curvature;
160 for(
int iPop=1; iPop < DESCRIPTOR::q; iPop++) {
161 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop, 0),
162 descriptors::c<DESCRIPTOR>(iPop, 1)});
163 int iPop_op = descriptors::opposite<DESCRIPTOR>(iPop);
171 for(
size_t u_i = 0; u_i < DESCRIPTOR::d; ++u_i){
178 dfs[iPop_op] = cell[iPop_op];
182 for(
int iPop=1; iPop<DESCRIPTOR::q; iPop++) {
183 cell[iPop] = dfs[iPop];
190 T rho = cell.computeRho();
193 if ( mass_tmp < -transition * rho || (mass_tmp < lonely_threshold * rho && !neighbour_info.
has_fluid_neighbours) ){
196 else if ( mass_tmp > (1. + transition)*rho || ( mass_tmp > (1-lonely_threshold) * rho && !neighbour_info.
has_gas_neighbours) ){
234 T u_avg[DESCRIPTOR::d] = {0., 0.};
238 for(
int iPop=1; iPop<DESCRIPTOR::q; iPop++) {
239 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop,0),
240 descriptors::c<DESCRIPTOR>(iPop,1)});
244 T u_tmp[DESCRIPTOR::d] = {0., 0.};
246 cellC.computeRhoU(rho_tmp, u_tmp);
248 for(
size_t i = 0; i < DESCRIPTOR::d; ++i){
249 u_avg[i] += u_tmp[i];
255 rho_avg /=
static_cast<T
>(ctr);
256 for(
size_t i = 0; i < DESCRIPTOR::d; ++i){
257 u_avg[i] /=
static_cast<T
>(ctr);
261 cell.iniEquilibrium(rho_avg, u_avg);
328 T rho = cell.computeRho();
329 T mass = cell.template getField<FreeSurface::MASS>( );
332 auto normal = computeParkerYoungInterfaceNormal(cell);
340 cell.template setField<FreeSurface::MASS>( 0. );
341 normal = {-normal[0], -normal[1]};
343 mass_excess = mass - rho;
344 cell.template setField<FreeSurface::MASS>( rho );
349 std::array<T,DESCRIPTOR::q> products;
352 std::size_t product_total = 0;
354 for(
int iPop=1; iPop<DESCRIPTOR::q; iPop++) {
355 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop,0),
356 descriptors::c<DESCRIPTOR>(iPop,1)});
366 products[iPop] = (normal[0] * descriptors::c<DESCRIPTOR>(iPop, 0) + normal[1] * descriptors::c<DESCRIPTOR>(iPop,1));
367 if(products[iPop] <= 0){
371 product_sum += products[iPop];
377 mass_excess_vector[0] = 0.;
379 if (product_total > 0) {
380 T product_fraction = 1. / product_total;
381 for(
int iPop=1; iPop < DESCRIPTOR::q; iPop++) {
382 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop,0),
383 descriptors::c<DESCRIPTOR>(iPop,1)});
388 mass_excess_vector[iPop] = mass_excess * product_fraction;
390 mass_excess_vector[iPop] = 0.;
393 cell.template setField<FreeSurface::TEMP_MASS_EXCHANGE>( mass_excess_vector );
395 mass_excess_vector[0] = mass_excess;
396 for(
int iPop=1; iPop < DESCRIPTOR::q; iPop++) {
397 mass_excess_vector[iPop] = 0.;
399 cell.template setField<FreeSurface::TEMP_MASS_EXCHANGE>( mass_excess_vector );
437 cell.template setField<FreeSurface::EPSILON>( 1. );
438 T mass_tmp = cell.template getField<FreeSurface::MASS>();
439 mass_tmp += cell.template getField<FreeSurface::TEMP_MASS_EXCHANGE>()[0];
440 cell.template setField<FreeSurface::MASS>(mass_tmp);
448 cell.template setField<FreeSurface::EPSILON>( 0. );
449 T mass_tmp = cell.template getField<FreeSurface::MASS>();
450 mass_tmp += cell.template getField<FreeSurface::TEMP_MASS_EXCHANGE>()[0];
451 cell.template setField<FreeSurface::MASS>(mass_tmp);
470 T collected_excess = 0.;
471 for(
int iPop = 1; iPop < DESCRIPTOR::q; ++iPop){
472 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop,0),
473 descriptors::c<DESCRIPTOR>(iPop,1)});
474 auto tempMassExchange = cellC.template getFieldPointer<FreeSurface::TEMP_MASS_EXCHANGE>();
476 if (hasCellFlags(cellC,
479 int iPop_op = descriptors::opposite<DESCRIPTOR>(iPop);
480 collected_excess += tempMassExchange[iPop_op];
484 T mass_tmp = cell.template getField<FreeSurface::MASS>();
486 mass_tmp += collected_excess;
488 T u_tmp[DESCRIPTOR::d];
489 cell.computeRhoU(rho, u_tmp);
493 cell.template setField<FreeSurface::EPSILON>( mass_tmp / rho );
494 cell.template setField<FreeSurface::MASS>(mass_tmp);
495 cell.template setField<FreeSurface::PREVIOUS_VELOCITY>(u_vel);
500 T collected_excess = 0.;
502 for(
int iPop = 1; iPop < DESCRIPTOR::q; ++iPop){
503 auto cellC = cell.neighbor({descriptors::c<DESCRIPTOR>(iPop,0),
504 descriptors::c<DESCRIPTOR>(iPop,1)});
505 auto tempMassExchange = cellC.template getFieldPointer<FreeSurface::TEMP_MASS_EXCHANGE>();
506 if (hasCellFlags(cellC,
509 int iPop_op = descriptors::opposite<DESCRIPTOR>(iPop);
510 collected_excess += tempMassExchange[iPop_op];
514 T mass_tmp = cell.template getField<FreeSurface::MASS>();
515 mass_tmp += collected_excess;
516 cell.template setField<FreeSurface::MASS>(mass_tmp);
534 sLattice.template addPostProcessor<FreeSurface::Stage0>(
536 sLattice.template addPostProcessor<FreeSurface::Stage1>(
538 sLattice.template addPostProcessor<FreeSurface::Stage2>(
540 sLattice.template addPostProcessor<FreeSurface::Stage3>(
542 sLattice.template addPostProcessor<FreeSurface::Stage4>(
547 communicator.requestOverlap(2);
548 communicator.template requestField<FreeSurface::EPSILON>();
549 communicator.template requestField<FreeSurface::CELL_TYPE>();
550 communicator.template requestField<descriptors::POPULATION>();
551 communicator.exchangeRequests();
557 communicator.requestOverlap(2);
558 communicator.template requestField<FreeSurface::CELL_FLAGS>();
559 communicator.template requestField<descriptors::POPULATION>();
560 communicator.exchangeRequests();
566 communicator.requestOverlap(2);
567 communicator.template requestField<FreeSurface::CELL_FLAGS>();
568 communicator.exchangeRequests();
574 communicator.requestOverlap(2);
575 communicator.template requestField<FreeSurface::CELL_FLAGS>();
576 communicator.exchangeRequests();
582 communicator.requestOverlap(2);
583 communicator.template requestField<FreeSurface::TEMP_MASS_EXCHANGE>();
584 communicator.exchangeRequests();
587 sLattice.template addCustomTask<stage::PostStream>([&]() {
Free Surface Processor 7 Finishes up left over cell conversions and prepares the state for the next s...
Free Surface Processor 6 Calculates mass excess from the cell type conversions and distributes them t...
void apply(CELL &cell, PARAMETERS ¶meters) any_platform
Free Surface Processor 5 ToGas Converts cells to interface from fluid if a neighbouring cell was conv...