Adding outflow to BreakingDam3D example case

[stonepreston]

I am attempting to change the outlet of the breaking dam to something that would allow outflow instead of bounce back, and could really use some help. The simplest outlet BC I thought could work would be ZeroGradient, however the simulation diverges when a zero gradient BC is applied to the outlet. (Note: id prefer to use the 2d case as its simpler but the zero gradient BC is only implemented for 3d it seems, so thats why Ive gone with the 3d example) Am I applying something incorrectly? I tried to change as little as possible from the default example

My steps are as follows:

1. Edit prepareGeometry to add a material number for the outlet cells:

void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
                      SuperGeometry<T,3>& superGeometry, const FreeSurfaceAppHelper& helper) {

  // ...
  // Untouched example code above
  // ...

  superGeometry.rename( 0,2 );
  superGeometry.rename( 2,1, {1,1,1} );

  Vector<T,3> outletExtend(1.*converter.getConversionFactorLength(), helper.area[1], helper.area[2]);
  Vector<T,3> outletOrigin;
  outletOrigin[0] = helper.area[0] - 1*converter.getConversionFactorLength();
  IndicatorCuboid3D<T> outlet( outletExtend, outletOrigin );
  superGeometry.rename( 2,5,1, outlet ); 
  
  // ...
  // Untouched example code above
  // ...
}

2. Initialize free surface fields for the new outlet material number (taking care to apply the gravity force to this material number as well as material number 1

void prepareBreakingDam(UnitConverter<T,DESCRIPTOR> const& converter,
                     SuperLattice<T, DESCRIPTOR>& sLattice,
                     SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper)
{
  // ...
  // Untouched example code above
  // ...

  auto outlet = superGeometry.getMaterialIndicator(5);
  sLattice.defineField<FreeSurface::CELL_TYPE>(outlet, zero);
  sLattice.defineField<FreeSurface::MASS>(outlet, zero);
  sLattice.defineField<FreeSurface::EPSILON>(outlet, zero);

  for (int i: {0,2}) {
    //sLattice.defineField<FreeSurface::MASS>(superGeometry, i, one);
    sLattice.defineField<FreeSurface::EPSILON>(superGeometry, i, one);
    sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, i, four);
  }

  T force_factor = 1./ converter.getConversionFactorForce() * converter.getConversionFactorMass();
  AnalyticalConst3D<T,T> force_a{helper.gravity_force[0] * force_factor, helper.gravity_force[1] * force_factor, helper.gravity_force[2] * force_factor};
  sLattice.defineField<descriptors::FORCE>(superGeometry.getMaterialIndicator({1,5}), force_a);

}

3. Prepare the lattice and apply the zerogradient BC:

void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
                     SuperLattice<T, DESCRIPTOR>& sLattice,
                     SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper) {
  // ...
  // Untouched example code above
  // ...

  // Material=0 -->do nothing
  sLattice.defineDynamics<NoDynamics<T,DESCRIPTOR>>(superGeometry, 0);
  // Material=1,5 -->bulk dynamics
  auto bulkIndicator = superGeometry.getMaterialIndicator({1, 5});
  sLattice.defineDynamics<SmagorinskyForcedBGKdynamics<T,DESCRIPTOR>>(bulkIndicator);
  // Material=2 -->no-slip boundary
  sLattice.defineDynamics<BounceBack<T,DESCRIPTOR>>( superGeometry, 2);
  //setSlipBoundary<T,DESCRIPTOR>(sLattice, superGeometry, 2);
  // Material 5 -> zero gradient boundary
  setZeroGradientBoundary<T,DESCRIPTOR>(sLattice, superGeometry, 5);

}

4. Set the initial values (basically just add the new material number to the list of materials

void setInitialValues(SuperLattice<T, DESCRIPTOR>& sLattice, SuperGeometry<T,3>& sGeometry, T lattice_length, UnitConverter<T,DESCRIPTOR> const& converter){
  OstreamManager clout( std::cout,"setInitialValues" );

  AnalyticalConst3D<T,T> u{0., 0.};
  AnalyticalConst3D<T,T> one(1.);

  sLattice.defineRhoU( sGeometry.getMaterialIndicator({0,1,2,5}), one, u );
  for (int i: {0,1,2,5}) {
    sLattice.iniEquilibrium( sGeometry, i, one, u );
  }

  // Set up free surface communicator stages
  FreeSurface::initialize(sLattice);
  // Make the lattice ready for simulation
  sLattice.initialize();
}

Am I missing something? I thought this would be relatively simple to implement but that is proving to not be the case and I could use some help.

I will include the entire file below:

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2021 Claudius Holeksa, Maximilian Schecher
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

#include "olb3D.h"
#include "olb3D.hh"

#include <vector>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <memory>
#include <fstream>
#include <sstream>

using namespace olb;
using namespace olb::descriptors;

using T = float;
using DESCRIPTOR = D3Q27<descriptors::FORCE, FreeSurface::MASS, FreeSurface::EPSILON, FreeSurface::CELL_TYPE, FreeSurface::CELL_FLAGS, FreeSurface::TEMP_MASS_EXCHANGE, FreeSurface::PREVIOUS_VELOCITY>;

/*
 * Helper since there are a lot of values to set and giving a reference to this object is easier than
 * defining the function calls accordingly
 */
struct FreeSurfaceAppHelper {
  std::array<T,3> area{{7.31,0.42,0.42}};
  std::array<T,3> gravity_force = {{0.,0., -9.81}};

  T char_phys_length = 0.03;
  T char_phys_vel = 0.1;
  bool has_surface_tension = true;
  T surface_tension_coefficient = 0.0661;

};

template<typename T, typename DESCRIPTOR>
class FreeSurfaceBreakingDam3D final : public AnalyticalF3D<T,T> {
private:
  T lattice_size;
  // This value doesn't depend on the dimension of the descriptor. It's always 3
  std::array<T,3> cell_values;
  std::array<T,3> area;
public:
  FreeSurfaceBreakingDam3D(T lattice_size_, const std::array<T,3>& cell_vals, const std::array<T,3>& area_):AnalyticalF3D<T,T>{1}, lattice_size{lattice_size_}, cell_values{cell_vals}, area{area_}{}

  // This sets which parts of the area is going to be the fluid
  bool operator()(T output[], const T x[]) override {
    output[0] = cell_values[0];

    if (x[2] <= area[2] * 0.95 &&
        x[1] <= area[1] * 0.95 &&
        x[0] <= area[0] * 0.45){
      output[0] = cell_values[2];
    } else if (x[2] - lattice_size * 1.1 <= area[2] * 0.95 &&
               x[1] - lattice_size * 1.1 <= area[1] * 0.95 &&
               x[0] - lattice_size * 1.1 <= area[0] * 0.45) {
      output[0] = cell_values[1];
    }

    return true;
  }

};

void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
                      SuperGeometry<T,3>& superGeometry, const FreeSurfaceAppHelper& helper) {

  OstreamManager clout( std::cout,"prepareGeometry" );
  clout << "Prepare Geometry ..." << std::endl;

  superGeometry.rename( 0,2 );
  superGeometry.rename( 2,1, {1,1,1} );

  Vector<T,3> outletExtend(1.*converter.getConversionFactorLength(), helper.area[1], helper.area[2]);
  Vector<T,3> outletOrigin;
  outletOrigin[0] = helper.area[0] - 1*converter.getConversionFactorLength();
  IndicatorCuboid3D<T> outlet( outletExtend, outletOrigin );
  superGeometry.rename( 2,5,1, outlet ); 

  superGeometry.clean();
  superGeometry.innerClean();
  superGeometry.checkForErrors();

  superGeometry.print();

  clout << "Prepare Geometry ... OK" << std::endl;
}

void prepareBreakingDam(UnitConverter<T,DESCRIPTOR> const& converter,
                     SuperLattice<T, DESCRIPTOR>& sLattice,
                     SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper)
{
  AnalyticalConst3D<T,T> zero( 0. );
  AnalyticalConst3D<T,T> one( 1. );
  AnalyticalConst3D<T,T> two( 2. );
  AnalyticalConst3D<T,T> four( 4. );
  FreeSurfaceBreakingDam3D<T,DESCRIPTOR> cells_analytical{ lattice_size, {0., 1., 2.}, helper.area};
  FreeSurfaceBreakingDam3D<T,DESCRIPTOR> mass_analytical{ lattice_size, {0., 0.5, 1.}, helper.area};

  AnalyticalConst3D<T,T> force_zero{0., 0., 0.};

  for (int i: {0,1,2}) {
    sLattice.defineField<FreeSurface::MASS>(superGeometry, i, zero);
    sLattice.defineField<FreeSurface::EPSILON>(superGeometry, i, zero);
    sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, i, zero);
    sLattice.defineField<FreeSurface::CELL_FLAGS>(superGeometry, i, zero);
    sLattice.defineField<descriptors::FORCE>(superGeometry, i, force_zero);
  }

  sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, 1, cells_analytical);

  sLattice.defineField<FreeSurface::MASS>(superGeometry, 1, mass_analytical);
  sLattice.defineField<FreeSurface::EPSILON>(superGeometry, 1, mass_analytical);

  auto outlet = superGeometry.getMaterialIndicator(5);
  sLattice.defineField<FreeSurface::CELL_TYPE>(outlet, zero);
  sLattice.defineField<FreeSurface::MASS>(outlet, zero);
  sLattice.defineField<FreeSurface::EPSILON>(outlet, zero);

  for (int i: {0,2}) {
    //sLattice.defineField<FreeSurface::MASS>(superGeometry, i, one);
    sLattice.defineField<FreeSurface::EPSILON>(superGeometry, i, one);
    sLattice.defineField<FreeSurface::CELL_TYPE>(superGeometry, i, four);
  }

  T force_factor = 1./ converter.getConversionFactorForce() * converter.getConversionFactorMass();
  AnalyticalConst3D<T,T> force_a{helper.gravity_force[0] * force_factor, helper.gravity_force[1] * force_factor, helper.gravity_force[2] * force_factor};
  sLattice.defineField<descriptors::FORCE>(superGeometry.getMaterialIndicator({1,5}), force_a);

}

void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
                     SuperLattice<T, DESCRIPTOR>& sLattice,
                     SuperGeometry<T,3>& superGeometry, T lattice_size, const FreeSurfaceAppHelper& helper) {

  OstreamManager clout( std::cout,"prepareLattice" );
  clout << "Prepare Lattice ..." << std::endl;

  // Material=0 -->do nothing
  sLattice.defineDynamics<NoDynamics<T,DESCRIPTOR>>(superGeometry, 0);
  // Material=1 -->bulk dynamics
  auto bulkIndicator = superGeometry.getMaterialIndicator({1, 5});
  sLattice.defineDynamics<SmagorinskyForcedBGKdynamics<T,DESCRIPTOR>>(bulkIndicator);
  // Material=2 -->no-slip boundary
  sLattice.defineDynamics<BounceBack<T,DESCRIPTOR>>( superGeometry, 2);
  //setSlipBoundary<T,DESCRIPTOR>(sLattice, superGeometry, 2);

  setZeroGradientBoundary<T,DESCRIPTOR>(sLattice, superGeometry, 5);

  sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());
  sLattice.setParameter<collision::LES::Smagorinsky>(T(0.2));

  prepareBreakingDam(converter, sLattice, superGeometry, lattice_size, helper);

  clout << "Prepare Lattice ... OK" << std::endl;
}

void setInitialValues(SuperLattice<T, DESCRIPTOR>& sLattice, SuperGeometry<T,3>& sGeometry, T lattice_length, UnitConverter<T,DESCRIPTOR> const& converter){
  OstreamManager clout( std::cout,"setInitialValues" );

  AnalyticalConst3D<T,T> u{0., 0.};
  AnalyticalConst3D<T,T> one(1.);

  sLattice.defineRhoU( sGeometry.getMaterialIndicator({0,1,2,5}), one, u );
  for (int i: {0,1,2,5}) {
    sLattice.iniEquilibrium( sGeometry, i, one, u );
  }

  // Set up free surface communicator stages
  FreeSurface::initialize(sLattice);
  // Make the lattice ready for simulation
  sLattice.initialize();
}

void getResults( SuperLattice<T,DESCRIPTOR>& sLattice,
                 UnitConverter<T,DESCRIPTOR> const& converter, int iT,
                 SuperGeometry<T,3>& superGeometry, util::Timer<T>& timer)
{
  const int vtmIter  = 100;
  const int statIter = 100;

  if ( iT==0 ) {
    // Writes the geometry, cuboid no. and rank no. as vti file for visualization
    SuperVTMwriter3D<T> vtmWriter( "breakingDam3d" );
    SuperLatticeGeometry3D<T, DESCRIPTOR> geometry( sLattice, superGeometry );
    SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid( sLattice );
    SuperLatticeRank3D<T, DESCRIPTOR> rank( sLattice );
    vtmWriter.write( geometry );
    vtmWriter.write( cuboid );
    vtmWriter.write( rank );

    vtmWriter.createMasterFile();
  }

  // Writes output in ParaView
  if ( iT%vtmIter==0 ) {
    sLattice.setProcessingContext(ProcessingContext::Evaluation);

    SuperVTMwriter3D<T> vtmWriter( "breakingDam3d" );
    SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity( sLattice, converter );
    SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure( sLattice, converter );
    SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::EPSILON> epsilon( sLattice );
    SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::CELL_TYPE> cells( sLattice );
    SuperLatticeExternalScalarField3D<T, DESCRIPTOR, FreeSurface::MASS> mass( sLattice );
    SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
    epsilon.getName() = "epsilon";
    cells.getName() = "cell_type";
    mass.getName() = "mass";
    vtmWriter.addFunctor( velocity );
    vtmWriter.addFunctor( pressure );
    vtmWriter.addFunctor( epsilon );
    vtmWriter.addFunctor( cells );
    vtmWriter.addFunctor( mass );
    vtmWriter.addFunctor( geometry );

    vtmWriter.write(iT);
  }

  // Writes output on the console
  if ( iT%statIter==0 ) {
    // Timer console output
    timer.update( iT );
    timer.printStep();

    // Lattice statistics console output
    sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) );
  }
}

namespace {
FreeSurfaceAppHelper free_surface_config {
  {7.31,0.42,0.42}, //Area
  {0.,0.,-9.81}, // Gravity force
  7.31, // char_phys_length
  0.1, // char_phys_vel
  true, // has_surface_tension
  0.05 // surface tension coefficient
};

class FreeSurfaceConfig {
public:
  T viscosity = 1e-4; // viscosity of water: 1e-4
  //T viscosity = 122000e-4; // viscosity of honey: 122000e-4 = 12.2
  T density = 1e3;  // density of water: 1e3
  //T density = 1.358e3;  // density of honey: 1.358e3 = 1358
  T physTime = 40.;
  T latticeRelaxationTime = .501;
  int N = 500;

  // Anti jitter value
  T transitionThreshold = 1e-3;
  // When to remove lonely cells
  T lonelyThreshold = 1.0;
};

}

int main(int argc, char **argv)
{
  olbInit(&argc, &argv, false, false);

  FreeSurfaceConfig c;
  OstreamManager clerr( std::cerr, "main" );
  OstreamManager clout( std::cout, "main" );

  singleton::directories().setOutputDir("./tmp/");

  FreeSurfaceAppHelper& helper = free_surface_config;

  UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
    int {c.N},     // resolution: number of voxels per charPhysL
    (T)   c.latticeRelaxationTime,   // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
    (T)   helper.char_phys_length,     // charPhysLength: reference length of simulation geometry
    (T)   helper.char_phys_vel,     // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
    (T)   c.viscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
    (T)   c.density     // physDensity: physical density in __kg / m^3__
  );

  // Prints the converter log as console output
  converter.print();
  // Writes the converter log in a file
  converter.write("free surface");

  T lattice_size = helper.char_phys_length / c.N;

  T force_conversion_factor = 1./converter.getConversionFactorForce()*converter.getConversionFactorMass();

  // Convert kg / s^2
  // Basically it is multiplied with s^2 / kg = s^2 * m^3 / (kg * m^2 * m) = 1. / (velocity_factor^2 * density * length_factor)
  T surface_tension_coefficient_factor = std::pow(converter.getConversionFactorTime(),2)/ (c.density * std::pow(converter.getConversionFactorLength(),3));

  clout<<"Surface: "<<surface_tension_coefficient_factor * helper.surface_tension_coefficient<<std::endl;
  clout<<"Lattice Size: "<<converter.getPhysDeltaX()<<std::endl;

  // === 2nd Step: Prepare Geometry ===
  Vector<T,3> extend( helper.area[0], helper.area[1], helper.area[2] );
  Vector<T,3> origin;
  IndicatorCuboid3D<T> cuboid( extend, origin );

  // Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
  const int noOfCuboids = singleton::mpi().getSize();
#else
  const int noOfCuboids = 4;
#endif
  CuboidGeometry3D<T> cuboidGeometry( cuboid, converter.getConversionFactorLength(), noOfCuboids );

  HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );
  SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer, 2 );

  prepareGeometry( converter, superGeometry, helper );

  SuperLattice<T, DESCRIPTOR> sLattice( superGeometry );

  clout<<"Overlap: "<<sLattice.getOverlap()<<std::endl;

  prepareLattice( converter, sLattice, superGeometry, lattice_size, helper);

  FreeSurface3DSetup<T,DESCRIPTOR> free_surface_setup{sLattice};

  free_surface_setup.addPostProcessor();

  // Set variables from freeSurfaceHelpers.h
  sLattice.setParameter<FreeSurface::DROP_ISOLATED_CELLS>(true);
  sLattice.setParameter<FreeSurface::TRANSITION>(c.transitionThreshold);
  sLattice.setParameter<FreeSurface::LONELY_THRESHOLD>(c.lonelyThreshold);
  sLattice.setParameter<FreeSurface::HAS_SURFACE_TENSION>(helper.has_surface_tension);
  sLattice.setParameter<FreeSurface::SURFACE_TENSION_PARAMETER>(surface_tension_coefficient_factor * helper.surface_tension_coefficient);
  sLattice.setParameter<FreeSurface::FORCE_CONVERSION_FACTOR>(force_conversion_factor);
  sLattice.setParameter<FreeSurface::LATTICE_SIZE>(converter.getPhysDeltaX());

  // === 4th Step: Main Loop with Timer ===
  clout << "starting simulation..." << std::endl;
  util::Timer<T> timer( converter.getLatticeTime( c.physTime ), superGeometry.getStatistics().getNvoxel() );
  timer.start();
  setInitialValues(sLattice, superGeometry, lattice_size, converter);

  for ( std::size_t iT = 0; iT < converter.getLatticeTime( c.physTime ); ++iT ) {
    getResults( sLattice, converter, iT, superGeometry, timer );
    sLattice.collideAndStream();
  }

  timer.stop();
  timer.printSummary();

  return 0;
}

• This topic was modified 1 month ago by stonepreston.

[mathias]

Such detailed support we cannot offer in the framework of the forum. You may consider coming to the next spring school or look into the options here https://www.openlb.net/consortium/.

[Author]

Posts