Reply To: Gallium Melting from 2D to 3D

[Jijo]

Hey Ramiro Freile,

Thank you for your help. I have done the conversion but the output text only shows 1 material number rather than 4. Can you please check out my code? I would appreciate that.

The code:
#include “olb3D.h”
#include “olb3D.hh”

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

using T = double;

using NSDESCRIPTOR = D3Q19<POROSITY,VELOCITY_SOLID,FORCE>;
using TDESCRIPTOR = D3Q7<VELOCITY,TEMPERATURE>;

using TotalEnthalpyAdvectionDiffusionDynamics = TotalEnthalpyAdvectionDiffusionTRTdynamics<T,TDESCRIPTOR>;

// Parameters for the simulation setup
const T lx = 88.9e-3; // length of the channel
const T ly = 63.5e-3; // height of the channel
const T lz = 38.1e-3;
int N = 128; // resolution of the model
T tau = 0.51; // relaxation time
const T Ra = 1.656e6; // Rayleigh number
const T Pr = 0.0216; // Prandtl number
const T Ste = 0.039; // Stephan number
const T maxPhysT = 1140.; // simulation time

const T Tcold = 0;
const T Tmelt = (302.8 – 301.3)/(311.0 – 301.3);
const T Thot = 1;

const T lambda_s = 40.6; // W / m K
const T lambda_l = 25.3; // W / m K
const T R_lambda = lambda_s/lambda_l;

const T cp_s = 1.0; // J / kg K
const T cp_l = 0.96; // J / kg K
const T R_cp = cp_s/cp_l;

//for this case, the harmonic mean (cp_ref) is applicable
const T cp_ref = 2.0 * cp_s * cp_l / (cp_s + cp_l); // J / kg K

const T R_alpha = lambda_s / lambda_l * cp_l / cp_s;

const T L = cp_l * (Thot – Tmelt) / Ste; // J / kg

T lattice_Hcold, lattice_Hhot;
T physDeltaX, physDeltaT;

/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry<T,3>& superGeometry,
ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter)
{

OstreamManager clout(std::cout,”prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0,4);

std::vector<T> extend(3,T());
extend[0] = lx;
extend[1] = ly;
extend[2] = lz;

std::vector<T> origin(3,T());
origin[0] = converter.getPhysLength(1);
origin[1] = 0.5*converter.getPhysLength(1);
origin[2] = 0.0;
IndicatorCuboid2D<T> cuboid2(extend, origin);

std::vector<T> extendwallleft(3,T(0));
extendwallleft[0] = converter.getPhysLength(1);
extendwallleft[1] = ly;
extendwallleft[2] = 0.1;

std::vector<T> originwallleft(3,T(0));
originwallleft[0] = 0.0;
originwallleft[1] = 0.0;
originwallleft[2] = 0.0;
IndicatorCuboid3D<T> wallleft(extendwallleft, originwallleft);

std::vector<T> extendwallright(3,T(0));
extendwallright[0] = converter.getPhysLength(1);
extendwallright[1] = ly;
extendwallright[2] = 0.1;

std::vector<T> originwallright(3,T(0));
originwallright[0] = lx+converter.getPhysLength(1);
originwallright[1] = 0.0;
originwallright[2] = 0.0;

IndicatorCuboid3D<T> wallright(extendwallright, originwallright);

superGeometry.rename(4,2,1,wallleft);
superGeometry.rename(4,3,1,wallright);

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

superGeometry.print();

clout << “Prepare Geometry … OK” << std::endl;

}

void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,3>& superGeometry )
{

OstreamManager clout(std::cout,”prepareLattice”);
clout << “Prepare Lattice …” << std::endl;

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

NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);
NSlattice.defineDynamics<ForcedPSMBGKdynamics>(superGeometry.getMaterialIndicator({1, 2, 3, 4}));

ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
ADlattice.defineDynamics<TotalEnthalpyAdvectionDiffusionDynamics>(superGeometry.getMaterialIndicator({1, 2, 3}));
ADlattice.defineDynamics<BounceBack>(superGeometry, 4);

/// sets boundary
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry.getMaterialIndicator({2, 3}));
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3}));

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

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
ADlattice.setParameter<collision::TRT::MAGIC>(0.25);
ADlattice.setParameter<TotalEnthalpy::T_S>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::T_L>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::CP_S>(cp_s);
ADlattice.setParameter<TotalEnthalpy::CP_L>(cp_l);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_S>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5) * R_lambda);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_L>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5));
ADlattice.setParameter<TotalEnthalpy::L>(L);

/// define initial conditions
AnalyticalConst3D<T,T> rho(1.);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(lattice_Hcold);
AnalyticalConst3D<T,T> T_hot(lattice_Hhot);

/// for each material set Rho, U and the Equilibrium
NSlattice.defineRhoU(superGeometry.getMaterialIndicator({1, 2, 3, 4}), rho, u0);
NSlattice.iniEquilibrium(superGeometry.getMaterialIndicator({1, 2, 3, 4}), rho, u0);

ADlattice.defineField<descriptors::VELOCITY>(superGeometry.getMaterialIndicator({1, 2, 3}), u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({1, 3}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({1, 3}), T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);

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

clout << “Prepare Lattice … OK” << std::endl;
}

void setBoundaryValues( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,3>& superGeometry)
{

// nothing to do here

}

void getResults( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,3>& superGeometry,
util::Timer<T>& timer,
bool converged)
{

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

SuperVTMwriter3D<T> vtkWriter(“thermalNaturalConvection3D”);
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeField3D<T, TDESCRIPTOR, VELOCITY> velocity(ADlattice);
SuperLatticePhysPressure3D<T, NSDESCRIPTOR> pressure(NSlattice, converter);

SuperLatticeDensity3D<T, TDESCRIPTOR> enthalpy(ADlattice);
enthalpy.getName() = “enthalpy”;
SuperLatticeField3D<T, NSDESCRIPTOR, POROSITY> liquid_frac(NSlattice);
liquid_frac.getName() = “liquid fraction”;
SuperLatticeField3D<T, TDESCRIPTOR, TEMPERATURE> temperature(ADlattice);
temperature.getName() = “temperature”;
SuperLatticeField3D<T, NSDESCRIPTOR, FORCE> force(NSlattice);
force.getName() = “force”;
vtkWriter.addFunctor( geometry );
vtkWriter.addFunctor( pressure );
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( enthalpy );
vtkWriter.addFunctor( liquid_frac );
vtkWriter.addFunctor( temperature );
vtkWriter.addFunctor( force );

const int vtkIter = converter.getLatticeTime(0.5);

if (iT == 0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid3D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank3D<T, NSDESCRIPTOR> rank(NSlattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);

vtkWriter.createMasterFile();
}

/// Writes the VTK files
if (iT % vtkIter == 0 || converged) {

timer.update(iT);
timer.printStep();

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

/// ADLattice statistics console output
ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

if ( NSlattice.getStatistics().getAverageRho() != NSlattice.getStatistics().getAverageRho() or ADlattice.getStatistics().getAverageRho() != ADlattice.getStatistics().getAverageRho() ) {
clout << “simulation diverged! stopping now.” << std::endl;
exit(1);
}
vtkWriter.write(iT);
}
}

int main(int argc, char *argv[])
{

/// === 1st Step: Initialization ===
OstreamManager clout(std::cout,”main”);
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);

T char_lattice_u = 0.2;

if (argc >= 2) {
N = atoi(argv[1]);
}
if (argc >= 3) {
tau = atof(argv[2]);
}
if (argc >= 4) {
char_lattice_u = atof(argv[3]);
}

const T char_u = util::sqrt( 9.81 * 1.2e-4 * (311. – 302.8) * 6093. );
const T conversion_u = char_u / char_lattice_u;

physDeltaX = lx / N;
physDeltaT = physDeltaX / conversion_u;
physDeltaT = 6093. / 1.81e-3 / descriptors::invCs2<T,NSDESCRIPTOR>() * (tau – 0.5) * physDeltaX * physDeltaX;

lattice_Hcold = cp_s * Tcold;
lattice_Hhot = cp_l * Thot;

clout << “H_cold ” << lattice_Hcold << ” H_hot ” << lattice_Hhot << std::endl;

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const converter(
(T) physDeltaX, // physDeltaX
(T) physDeltaT, // physDeltaT
(T) lx, // charPhysLength
(T) char_u, // charPhysVelocity
(T) 1.81e-3 / 6093., // physViscosity
(T) 6093., // physDensity
(T) 32., // physThermalConductivity
(T) 381., // physSpecificHeatCapacity
(T) 1.2e-4, // physThermalExpansionCoefficient
(T) 302.8, // charPhysLowTemperature
(T) 311 // charPhysHighTemperature
);
converter.print();
clout << “lattice cp ” << converter.getLatticeSpecificHeatCapacity(cp_l) << std::endl;

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

std::vector<T> extend(3,T());
extend[0] = lx + 2*converter.getPhysLength(1);
extend[1] = ly + converter.getPhysLength(1);
extend[2] = lz + 2*converter.getPhysLength(1);
std::vector<T> origin(3,T());
IndicatorCuboid3D<T> cuboid(extend, origin);

/// Instantiation of a cuboidGeometry with weights
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());
cuboidGeometry.setPeriodicity(false, false, true); // x, y, z

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

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

prepareGeometry(superGeometry, converter);

/// === 3rd Step: Prepare Lattice ===

SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

std::vector<T> dir{0.0,1.0,0.0};

T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

NavierStokesAdvectionDiffusionCouplingGenerator3D<T,NSDESCRIPTOR>
coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(lx),
boussinesqForcePrefactor, Tcold, 1., dir);

NSlattice.addLatticeCoupling(superGeometry, 1, coupling, ADlattice);

//prepareLattice and setBoundaryConditions
prepareLattice(converter, NSlattice, ADlattice, superGeometry);

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

for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT)+1; ++iT) {

/// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

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

/// === 7th Step: Computation and Output of the Results ===
getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, false);
}

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

}