Temperature blow up in NS-ADE coupling

[Anonymous]

Hi,

I am modifying the laminar 3D cylinder example code such that I can set the cylinder wall at constant temperature which is different from fluid temperature to observe convection.

I am not able to understand the reason why my temperature solution appears to blow up.

Following is the full code: I appreciate your help.

Best, P Shenai

PS: this is my first NS-ADE-LBM code.

#ifndef CYLINDER_3D_H
#define CYLINDER_3D_H

#include <olb.h>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;

using T = FLOATING_POINT_TYPE;

using NSDESCRIPTOR = D3Q19<FORCE>;
using TDESCRIPTOR =  D3Q7<VELOCITY>;

#define BOUZIDI

// Parameters for the simulation setup
const T Re = 20.;       // Reynolds number
const T maxPhysT = 16.; // max. simulation time in s, SI unit
const T Twall = 283.15;
const T Tfluid = 273.15;
const T viscosity =1.505e-5;			// kinematic viscosity 
const T density = 1.2041;			// Fluid density (kg/m^3)
const T conductivity = 0.0255;			// Fluid thermal conductivity (W/m-K)
const T Pr = 0.713;		// Prandtl number
const T expansionCoeff = 3.66e-3;				//
const T charL = 0.1;

// Stores data from stl file in geometry in form of material numbers
void prepareGeometry( const ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR>& converter, IndicatorF3D<T>& indicator,
                      STLreader<T>& stlReader, SuperGeometry<T,3>& superGeometry )
{

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

  superGeometry.rename( 0,2,indicator );
  superGeometry.rename( 2,1,stlReader );
  superGeometry.clean();

  Vector<T,3> origin = superGeometry.getStatistics().getMinPhysR( 2 );
  origin[1] += converter.getPhysDeltaX()/2.;
  origin[2] += converter.getPhysDeltaX()/2.;

  Vector<T,3> extend = superGeometry.getStatistics().getMaxPhysR( 2 );
  extend[1] = extend[1]-origin[1]-converter.getPhysDeltaX()/2.;
  extend[2] = extend[2]-origin[2]-converter.getPhysDeltaX()/2.;

  // Set material number for inflow
  origin[0] = superGeometry.getStatistics().getMinPhysR( 2 )[0]-converter.getPhysDeltaX();
  extend[0] = 2*converter.getPhysDeltaX();
  IndicatorCuboid3D<T> inflow( extend,origin );
  superGeometry.rename( 2,3,inflow );

  // Set material number for outflow
  origin[0] = superGeometry.getStatistics().getMaxPhysR( 2 )[0]-converter.getPhysDeltaX();
  extend[0] = 2*converter.getPhysDeltaX();
  IndicatorCuboid3D<T> outflow( extend,origin );
  superGeometry.rename( 2,4,outflow );

  // Set material number for cylinder
  origin[0] = superGeometry.getStatistics().getMinPhysR( 2 )[0]+converter.getPhysDeltaX();
  extend[0] = ( superGeometry.getStatistics().getMaxPhysR( 2 )[0]-superGeometry.getStatistics().getMinPhysR( 2 )[0] )/2.;
  std::shared_ptr<IndicatorF3D<T>> cylinder = std::make_shared<IndicatorCuboid3D<T>>( extend, origin );
  superGeometry.rename( 2,5, cylinder );

  // Removes all not needed boundary voxels outside the surface
  superGeometry.clean();
  superGeometry.checkForErrors();

  superGeometry.print();

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

// Set up the geometry of the simulation
void prepareLattice( SuperLattice<T, NSDESCRIPTOR>& NSlattice,
                    SuperLattice<T, TDESCRIPTOR>& ADlattice,
                     ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
                     STLreader<T>& stlReader,
                     SuperGeometry<T,3>& superGeometry )
{
  OstreamManager clout( std::cout,"prepareLattice" );
  clout << "Prepare Lattice ..." << std::endl;

  const T omega = converter.getLatticeRelaxationFrequency();
  const T Tomega  = converter.getLatticeThermalRelaxationFrequency();

  // Material=1 -->bulk dynamics
  auto bulkIndicator = superGeometry.getMaterialIndicator({1});
  NSlattice.defineDynamics<BGKdynamics>(bulkIndicator);

  // Material=2 -->bounce back
  boundary::set<boundary::BounceBack>(NSlattice, superGeometry, 2);

  // Setting of the boundary conditions

  // if local boundary conditions are chosen
  //boundary::set<boundary::LocalVelocity>(sLattice, superGeometry, 3);
  //boundary::set<boundary::LocalPressure>(sLattice, superGeometry, 4);

  //if interpolated boundary conditions are chosen
  boundary::set<boundary::InterpolatedVelocity>(NSlattice, superGeometry, 3);
  boundary::set<boundary::InterpolatedPressure>(NSlattice, superGeometry, 4);

  // Material=5 -->bouzidi / bounce back
  //#ifdef BOUZIDI
  setBouzidiBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 5, stlReader);
  //#else
  //boundary::set<boundary::BounceBack>(NSlattice, superGeometry, 5);
  //#endif

  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 3);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 4);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 5);

  setBouzidiBoundary<T,TDESCRIPTOR,BouzidiAdeDirichletPostProcessor>(ADlattice, superGeometry, 5, stlReader);

  // Initial conditions
  AnalyticalConst3D<T,T> rhoF(converter.getLatticeDensity(density));
  AnalyticalConst3D<T,T> uF(0, 0);
  AnalyticalConst3D<T,T> T_wall(converter.getLatticeTemperature(Twall));
  AnalyticalConst3D<T,T> T_fluid(converter.getLatticeTemperature(Tfluid));

  // Initialize all values of distribution functions to their local equilibrium
  NSlattice.defineRhoU(superGeometry, 1, rhoF, uF);
  NSlattice.iniEquilibrium(superGeometry, 1, rhoF, uF);

  ADlattice.defineRho(superGeometry, 1, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 1, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 2, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 2, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 3, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 3, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 4, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 4, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 5, T_wall);
  ADlattice.iniEquilibrium(superGeometry, 5, T_wall, uF);
  setBouzidiAdeDirichlet(ADlattice, superGeometry, 5, T_wall);

  NSlattice.setParameter<descriptors::OMEGA>(omega);
  ADlattice.setParameter<descriptors::OMEGA>(Tomega);

  // Make the lattice ready for simulation
  ADlattice.initialize();
  ADlattice.initialize();

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

// Generates a slowly increasing inflow for the first iTMaxStart timesteps
void setBoundaryValues( SuperLattice<T, NSDESCRIPTOR>& NSlattice,
                        SuperLattice<T, TDESCRIPTOR>& ADlattice,
                        ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter, int iT,
                        SuperGeometry<T,3>& superGeometry )
{
  OstreamManager clout( std::cout,"setBoundaryValues" );

  // No of time steps for smooth start-up
  int iTmaxStart = converter.getLatticeTime( maxPhysT*0.4 );
  int iTupdate = 30;

  if ( iT%iTupdate == 0 && iT <= iTmaxStart ) {
    // Smooth start curve, sinus
    // SinusStartScale<T,int> StartScale(iTmaxStart, T(1));

    // Smooth start curve, polynomial
    PolynomialStartScale<T,int> StartScale( iTmaxStart, T( 1 ) );

    // Creates and sets the Poiseuille inflow profile using functors
    int iTvec[1] = {iT};
    T frac[1] = {};
    StartScale( frac,iTvec );
    std::vector<T> maxVelocity( 3,0 );
    maxVelocity[0] = 2.25*frac[0]*converter.getCharLatticeVelocity();

    T distance2Wall = converter.getPhysDeltaX()/2.;
    RectanglePoiseuille3D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall, distance2Wall, distance2Wall );
    NSlattice.defineU( superGeometry, 3, poiseuilleU );

    clout << "step=" << iT << "; maxVel=" << maxVelocity[0] << std::endl;

    NSlattice.setProcessingContext<Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
      ProcessingContext::Simulation);
  }
}

// Computes the pressure drop between the voxels before and after the cylinder
double getResults( SuperLattice<T, NSDESCRIPTOR>& NSlattice,
                   SuperLattice<T, TDESCRIPTOR>& ADlattice,
                   ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter, int iT,
                   SuperGeometry<T,3>& superGeometry, util::Timer<T>& timer,
                   STLreader<T>& stlReader,
                   Gnuplot<T>& gplot,
                   bool eoc )
{

  OstreamManager clout( std::cout,"getResults" );

  SuperLatticePhysVelocity3D<T, NSDESCRIPTOR> velocity( NSlattice, converter );
  SuperLatticePhysPressure3D<T, NSDESCRIPTOR> pressure( NSlattice, converter );
  //SuperLatticeYplus3D<T, DESCRIPTOR> yPlus( NSlattice, converter, superGeometry, stlReader, 5 );
  //SuperLatticeRefinementMetricKnudsen3D<T, NSDESCRIPTOR> quality( NSlattice, converter );
  //SuperRoundingF3D<T, T> roundedQuality ( quality, RoundingMode::NearestInteger );
  //SuperDiscretizationF3D<T> discretization ( roundedQuality, 0., 2. );
  SuperLatticePhysTemperature3D<T, NSDESCRIPTOR,TDESCRIPTOR> temperature(ADlattice, converter);

  const int vtkIter  = converter.getLatticeTime( .3 );
  const int statIter = converter.getLatticeTime( .1 );

  T dragCoefficient = 0.;

  if (!eoc) {
    SuperVTMwriter3D<T> vtmWriter( "cylinder3d" );
    //vtmWriter.addFunctor( quality );
    vtmWriter.addFunctor( velocity );
    vtmWriter.addFunctor( pressure );
   // vtmWriter.addFunctor( yPlus );
    vtmWriter.addFunctor( temperature );
    if ( iT==0 ) {
        // Writes the geometry, cuboid no. and rank no. as vti file for visualization
        SuperLatticeCuboid3D<T, NSDESCRIPTOR> cuboid( NSlattice );
        SuperLatticeRank3D<T, NSDESCRIPTOR> rank( NSlattice );
        vtmWriter.write( cuboid );
        vtmWriter.write( rank );

        vtmWriter.createMasterFile();
    }

    // Writes the vtk files
    if (iT%vtkIter == 0) {
        vtmWriter.write(iT);

        {
        SuperEuklidNorm3D<T> normVel( velocity );
        BlockReduction3D2D<T> planeReduction( normVel, Vector<T,3>({0, 0, 1}) );
        // write output as JPEG
        heatmap::write(planeReduction, iT);
        }

        //{
        //BlockReduction3D2D<T> planeReduction( discretization, Vector<T,3>({0, 0, 1}) );
        //heatmap::plotParam<T> jpeg_scale;
        //jpeg_scale.colour = "blackbody";
        //jpeg_scale.name = "quality";
        //heatmap::write( planeReduction, iT, jpeg_scale );
       // }
    }
  }

  // Writes output on the console
  if ( iT%statIter == 0 ) {
    NSlattice.setProcessingContext(ProcessingContext::Evaluation);

    // Timer console output
    timer.update( iT );
    timer.printStep();

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

    // Drag, lift, pressure drop
    AnalyticalFfromSuperF3D<T> intpolatePressure( pressure, true );
    SuperLatticePhysDrag3D<T,NSDESCRIPTOR> drag( NSlattice, superGeometry, 5, converter );

    olb::Vector<T, 3> point1V = superGeometry.getStatistics().getCenterPhysR( 5 );
    olb::Vector<T, 3> point2V = superGeometry.getStatistics().getCenterPhysR( 5 );
    T point1[3] = {};
    T point2[3] = {};
    for ( int i = 0; i<3; i++ ) {
      point1[i] = point1V[i];
      point2[i] = point2V[i];
    }
    point1[0] = superGeometry.getStatistics().getMinPhysR( 5 )[0] - converter.getPhysDeltaX();
    point2[0] = superGeometry.getStatistics().getMaxPhysR( 5 )[0] + converter.getPhysDeltaX();

    T p1, p2;
    intpolatePressure( &p1,point1 );
    intpolatePressure( &p2,point2 );

    clout << "pressure1=" << p1;
    clout << "; pressure2=" << p2;

    T pressureDrop = p1-p2;
    clout << "; pressureDrop=" << pressureDrop;

    T dragA[3];
    int input1[0];
    drag( dragA, input1 );
    clout << "; drag=" << dragA[0] << "; lift=" << dragA[1] << std::endl;

    dragCoefficient = dragA[0];

    //int input[4] = {};
    //SuperMax3D<T> yPlusMaxF( yPlus, superGeometry, 1 );
    //T yPlusMax[1];
    //yPlusMaxF( yPlusMax,input );
    //clout << "yPlusMax=" << yPlusMax[0] << std::endl;

    if (!eoc) {
      // set data for gnuplot: input={xValue, yValue(s), names (optional), position of key (optional)}
      gplot.setData( converter.getPhysTime( iT ), {dragA[0], 5.58}, {"drag(openLB)", "drag(schaeferTurek)"}, "bottom right", {'l','l'} );

      // every (iT%vtkIter) write an png of the plot
      if ( iT%( vtkIter ) == 0 ) {
        // writes pngs: input={name of the files (optional), x range for the plot (optional)}
        gplot.writePNG( iT, maxPhysT );
      }
    }
  }

  return dragCoefficient;

}

double simulateCylinder( int N, Gnuplot<T>& gplot, bool eoc )
{

  // === 1st Step: Initialization ===
  singleton::directories().setOutputDir( "./tmp/" );
  OstreamManager clout( std::cout,"main" );
  // display messages from every single mpi process
  //clout.setMultiOutput(true);

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

ThermalUnitConverter< T, NSDESCRIPTOR, TDESCRIPTOR > converter(	
	(T)	charL/N,
	(T)	3.46e-4, // CFL*L/0.2,        //(tau - 0.5) / descriptors::invCs2<T,NSDESCRIPTOR>() * util::pow((charL/N),2) / viscosity, //  //
	(T)	charL,
	(T)	Re*viscosity / charL,
	(T)	viscosity,
	(T)	density,
	(T)	conductivity,
	(T)	Pr * conductivity / (viscosity * density),
  (T)	expansionCoeff,
	(T)	Tfluid,
	(T)	Twall
  );

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

  // === 2nd Step: Prepare Geometry ===

  // Instantiation of the STLreader class
  // file name, voxel size in meter, stl unit in meter, outer voxel no., inner voxel no.
  STLreader<T> stlReader( "../../laminar/cylinder3d/cylinder3d.stl", converter.getPhysDeltaX(), 0.001);
  IndicatorLayer3D<T> extendedDomain( stlReader, converter.getPhysDeltaX() );

  // Instantiation of a cuboidDecomposition with weights
#ifdef PARALLEL_MODE_MPI
  const int noOfCuboids = singleton::mpi().getSize();
#else
  const int noOfCuboids = 7;
#endif
  CuboidDecomposition3D<T> cuboidDecomposition( extendedDomain, converter.getPhysDeltaX(), noOfCuboids );

  // Instantiation of a loadBalancer
  HeuristicLoadBalancer<T> loadBalancer( cuboidDecomposition );

  // Instantiation of a superGeometry
  SuperGeometry<T,3> superGeometry( cuboidDecomposition, loadBalancer );

  prepareGeometry( converter, extendedDomain, stlReader, superGeometry );

  // === 3rd Step: Prepare Lattice ===
  SuperLattice<T, NSDESCRIPTOR> NSlattice( superGeometry );
  SuperLattice<T, TDESCRIPTOR> ADlattice( superGeometry );

  //prepareLattice and set boundaryCondition
  prepareLattice( NSlattice, ADlattice, converter, stlReader, superGeometry );

  T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
                               converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();
  SuperLatticeCoupling coupling(
    NavierStokesAdvectionDiffusionCoupling{},
    names::NavierStokes{}, NSlattice,
    names::Temperature{}, ADlattice);
  coupling.setParameter<NavierStokesAdvectionDiffusionCoupling::T0>(
    converter.getLatticeTemperature(Tfluid));
  coupling.setParameter<NavierStokesAdvectionDiffusionCoupling::FORCE_PREFACTOR>(
    boussinesqForcePrefactor * Vector<T,3>{0.0,1.0,0.0});

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

  T drag = 0;

  for (std::size_t iT = 0; iT < converter.getLatticeTime( maxPhysT ); ++iT) {
    // === 5th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues( NSlattice, ADlattice, converter, iT, superGeometry );

    // === 6th Step: Collide and Stream Execution ===
    ADlattice.collideAndStream();
    NSlattice.collideAndStream();

    coupling.execute();

    // === 7th Step: Computation and Output of the Results ===
    if (iT % converter.getLatticeTime( .1 ) == 0)
    {
      drag = getResults( NSlattice, ADlattice, converter, iT, superGeometry, timer, stlReader, gplot, eoc );
    }
  }

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

  return drag;
}

#endif

[FBukreev]

Hi,

I would add zeroGradient boundary conditions for ADE lattice to the inlet and outlet of the channel and in the coupler you can comment out the Boussinesq force if you dont consider it.

[Anonymous]

Hi,

Thanks for replying.

I did try the following modification, but hasn’t led to any improvement.

Appreciate your help.

Best, P Shenai

void prepareLattice( SuperLattice<T, NSDESCRIPTOR>& NSlattice,
                    SuperLattice<T, TDESCRIPTOR>& ADlattice,
                     ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
                     STLreader<T>& stlReader,
                     SuperGeometry<T,3>& superGeometry )
{
  OstreamManager clout( std::cout,"prepareLattice" );
  clout << "Prepare Lattice ..." << std::endl;

  const T omega = converter.getLatticeRelaxationFrequency();
  const T Tomega  = converter.getLatticeThermalRelaxationFrequency();

  // Material=1 -->bulk dynamics
  auto bulkIndicator = superGeometry.getMaterialIndicator({1});
  NSlattice.defineDynamics<BGKdynamics>(bulkIndicator);

  // Material=2 -->bounce back
  boundary::set<boundary::BounceBack>(NSlattice, superGeometry, 2);

  // Setting of the boundary conditions

  // if local boundary conditions are chosen
  //boundary::set<boundary::LocalVelocity>(sLattice, superGeometry, 3);
  //boundary::set<boundary::LocalPressure>(sLattice, superGeometry, 4);

  //if interpolated boundary conditions are chosen
  boundary::set<boundary::InterpolatedVelocity>(NSlattice, superGeometry, 3);
  boundary::set<boundary::InterpolatedPressure>(NSlattice, superGeometry, 4);

  // Material=5 -->bouzidi / bounce back
  //#ifdef BOUZIDI
  setBouzidiBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 5, stlReader);
  //#else
  //boundary::set<boundary::BounceBack>(NSlattice, superGeometry, 5);
  //#endif

  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
  //ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 3);
  setZeroGradientBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry, 3);
  //ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 4);
  setZeroGradientBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry, 4);
  ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 5);

  setBouzidiBoundary<T,TDESCRIPTOR,BouzidiAdeDirichletPostProcessor>(ADlattice, superGeometry, 5, stlReader);

  // Initial conditions
  AnalyticalConst3D<T,T> rhoF(converter.getLatticeDensity(density));
  AnalyticalConst3D<T,T> uF(0, 0);
  AnalyticalConst3D<T,T> T_wall(converter.getLatticeTemperature(Twall));
  AnalyticalConst3D<T,T> T_fluid(converter.getLatticeTemperature(Tfluid));

  // Initialize all values of distribution functions to their local equilibrium
  NSlattice.defineRhoU(superGeometry, 1, rhoF, uF);
  NSlattice.iniEquilibrium(superGeometry, 1, rhoF, uF);

  ADlattice.defineRho(superGeometry, 1, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 1, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 2, T_fluid);
  ADlattice.iniEquilibrium(superGeometry, 2, T_fluid, uF);
  //ADlattice.defineRho(superGeometry, 3, T_fluid);
  //ADlattice.iniEquilibrium(superGeometry, 3, T_fluid, uF);
  //ADlattice.defineRho(superGeometry, 4, T_fluid);
  //ADlattice.iniEquilibrium(superGeometry, 4, T_fluid, uF);
  ADlattice.defineRho(superGeometry, 5, T_wall);
  ADlattice.iniEquilibrium(superGeometry, 5, T_wall, uF);
  setBouzidiAdeDirichlet(ADlattice, superGeometry, 5, T_wall);

  NSlattice.setParameter<descriptors::OMEGA>(omega);
  ADlattice.setParameter<descriptors::OMEGA>(Tomega);

  // Make the lattice ready for simulation
  ADlattice.initialize();
  ADlattice.initialize();

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

[FBukreev]

sorry I cannot understand the reason for the problem on that way. It could be also numerical paameters that you use. The divergence time step should be also investigated carefully.

[Anonymous]

Thank you for your help! Does the modified prepareLattice block of code look right? I made the changes you recommended, but I am not sure if its implemented properly.

Best,
P Shenai

[mathias]

We can not give such details support here. But you may consider https://www.openlb.net/consortium/ .

[Author]

Posts