Welcome! Log In Create A New Profile

Advanced

need help with external Force

Posted by ChoelzlChoelzl  
need help with external Force
November 26, 2009 10:08AM
Hello guys,

I´m new to palabos and finally got stuck at some point:

My predecessor worked with openLB and used the following method to apply an external Force:

lattice.get(iX,iY).defineExternalField (
                                                 DESCRIPTOR<T>::ExternalField::forceBeginsAt,
                                                 DESCRIPTOR<T>::ExternalField::sizeOfForce,
                                                 force );


I tried to do something similar with palabos, and modified the code from the example rayleighTaylor2D.cpp - were a gravity field is applied.
So, after a few changes my code looks like:

template<typename T, template<typename U> class Descriptor>
class ExternalForceInitializer : public OneCellIndexedFunctional2D<T,Descriptor> {
public:
    ExternalForceInitializer(T force_)
        : force(force_)
    { }
	
    ExternalForceInitializer<T,Descriptor>* clone() const {
        return new ExternalForceInitializer<T,Descriptor>(*this);
    }
    virtual void execute(plint iX, plint iY, Cell<T,Descriptor>& cell) const {
        Array<T,2> force_x      (force, 0.);
        Array<T,2> zeroVelocity (0.,0.);
        T rho = (T)1.;
		
        iniCellAtEquilibrium(cell, rho, zeroVelocity);
		// Initialize the value of the force in the external scalars.
        cell.setExternalField (
               Descriptor<T>::ExternalField::forceBeginsAt,
               Descriptor<T>::ExternalField::sizeOfForce, &force_x[0] );

    }
private:
    T force;
};

// and in the main-function:

T force        = 0.15/(T)ny;

MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), 
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );
lattice.periodicity().toggle(0, true);
defineCylinderGeometry(lattice, parameters);
applyIndexed(lattice, Box2D(0, nx-1, 0, ny-1), new ExternalForceInitializer<T,DESCRIPTOR>(force) );

So far so good, the only problem is that there´s no movement at all, regardless of the value of <force>.

Any help is appreciated :)

chris

edit:and I´m using the
#define DESCRIPTOR plb::descriptors::ForcedD2Q9Descriptor

edit2: well it seems it´s only possible to compile with
#define DESCRIPTOR plb::descriptors::ForcedShanChenD2Q9Descriptor
and not the descriptor from my first edit, although I´m sure it did compile without errors with the above one on my first try.

edit3:
It seems i still need to apply a data-processing functional, but which one exactly ? So far, everything i tried (from the example i mentioned above or [www.lbmethod.org]) didn´t work.



Edited 3 time(s). Last edit at 11/26/2009 02:58PM by Choelzl.
Re: need help with external Force
November 27, 2009 01:28PM
I tried the following

#define DESCRIPTOR plb::descriptors::ForcedD2Q9Descriptor

//and

MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), 
                                                                                   new ForcedBGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

but it seems thats from the old openLB project, so far i can´t find the (renamed?) versions for palabos.....
any hints in the right direction are very welcome
Re: need help with external Force
November 27, 2009 02:56PM
Chris,

Your code looks perfectly OK, you almost got it. And, chosing the Descriptor

#define DESCRIPTOR plb::descriptors::ForcedD2Q9Descriptor

was the right thing to do. The last point you are missing is the dynamics. Instead of BGKdynamics, you should use GuoExternalForceBGKdynamics and the code will run smoothly.

This dynamics implements the Guo external force terms. More forcing terms will be available in the next Palabos release.
Re: need help with external Force
November 30, 2009 10:44AM
Thanks a lot jonas, everything is working fine now.

is there a table about all available dynamics somewhere - and which descriptor can be used with them (or the other way around) ?
Re: need help with external Force
November 30, 2009 10:54AM
Hi,

The most complete table is provided in the section "Implemented fluid models" in the user's guide:

[www.lbmethod.org]

Incidentally, I just noticed that the name of the descriptor to use with GuoExternalForceBGKdynamics is missing. It's been added for the next release.
Re: need help with external Force
March 19, 2010 06:00PM
Hi Choelzl,

Could you put here you source code which works?

Olivier.
Re: need help with external Force
April 08, 2010 07:57AM
Hi Olivier,

yes, I can put the code up here, but at the moment I have problems with my results. Once they are fixed I can post the code here.

Cheers,

Chris
Re: need help with external Force
May 12, 2010 09:22PM
Hi Chris,

Any success with your code?

Olivier.
Re: need help with external Force
May 17, 2010 11:41AM
Hi Olivier,

sorry I completely forgot to post it.

This is a simple example, I just compiled it and it should work.
For the real simulation take care to get usefull values for omega/tau !!

#include "palabos2D.h"
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
  #include "palabos2D.hh"   // include full template code
#endif

#include "poiseuille.h"
#include "poiseuille.hh"

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



using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR descriptors::ForcedD2Q9Descriptor

template<typename T, template<typename U> class Descriptor>
class ExternalForceInitializer : public OneCellIndexedFunctional2D<T,Descriptor> {
public:
    ExternalForceInitializer(T force_)
        : force(force_)
    { }
	
    ExternalForceInitializer<T,Descriptor>* clone() const {
        return new ExternalForceInitializer<T,Descriptor>(*this);
    }
    virtual void execute(plint iX, plint iY, Cell<T,Descriptor>& cell) const {
        T force_x = force;
	T force_y = 0;
        Array<T,2> forcefield (force_x, force_y);
	Array<T,2> zeroVelocity (0.,0.);
        T rho = (T)1.;
		
        //pcout << "Force: Fx = " << force_x << ", Fy = " << force_y << endl;
	iniCellAtEquilibrium(cell, rho, zeroVelocity);
		// Initialize the value of the force in the external scalars.
        cell.setExternalField (
               Descriptor<T>::ExternalField::forceBeginsAt,
               Descriptor<T>::ExternalField::sizeOfForce, &forcefield[0] );

    }
private:
    T force;
};

template<typename T>
class CylinderShapeDomain2D : public plb::DomainFunctional2D {
public:
    CylinderShapeDomain2D(plb::plint cx_, plb::plint cy_, plb::plint radius)
        : cx(cx_),
          cy(cy_),
          radiusSqr(plb::util::sqr(radius))
    { }
    virtual bool operator() (plb::plint iX, plb::plint iY) const {
        return plb::util::sqr(iX-cx) + plb::util::sqr(iY-cy) <= radiusSqr;
    }
    virtual CylinderShapeDomain2D<T>* clone() const {
        return new CylinderShapeDomain2D<T>(*this);
    }
private:
    plb::plint cx;
    plb::plint cy;
    plb::plint radiusSqr;
};

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

	global::directories().setOutputDir("/your/directory/path/");
	ostringstream value;

	
    // Define numeric parameters.
    IncomprFlowParam<T> parameters (
        (T) 1e-2,  // uMax
        (T) 10.,  // Re
        300,        // N
        1,        // lx
        1         // ly 
    );

	plint N = parameters.getResolution();
	plint nx = parameters.getNx();
	plint ny = parameters.getNy();
	pcout << "N: " << N << '\n';
	pcout << "Nx: " << nx << '\n';
	pcout << "Ny: " << ny << '\n';
	T force        = 0.015/(T)nx;	
	value << force;
	string title = "External Force applied on each lattice node: " + value.str();

	writeLogFile(parameters, title);
	
	
    MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), 
		new GuoExternalForceBGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

    
	
	lattice.periodicity().toggleAll(true);
	
	
    
	applyIndexed(lattice, Box2D(0, nx-1, 0, ny-1), new ExternalForceInitializer<T,DESCRIPTOR>(force) );
        
	
	defineDynamics (lattice, lattice.getBoundingBox(), new CylinderShapeDomain2D<T>(150, 150, 30), new BounceBack<T,DESCRIPTOR>);
	
	
	
    lattice.initialize();

    // Main loop over time iterations.
    for (plint iT=0; iT<1500; ++iT) {
        if (iT%100==0) {
            pcout << "Writing image at dimensionless time " << iT*parameters.getDeltaT() << endl;
            ImageWriter<T> imageWriter("leeloo");
            imageWriter.writeScaledGif (
                    createFileName("velocity", iT, 6),
                    *computeVelocityNorm(lattice));
		    VtkImageOutput2D<T> vtkOut(createFileName("vtk", iT, 6), parameters.getDeltaX());
			vtkOut.writeData<2,double>(*computeVelocity(lattice), "velocity", parameters.getDeltaX()/parameters.getDeltaT());

            pcout << computeAverageEnergy(lattice) << endl;
	    
        }
        // Lattice Boltzmann iteration step.
        lattice.collideAndStream();
    }

    //delete boundaryCondition;


}
Re: need help with external Force
January 16, 2018 09:02AM
Hi!

I am trying to implement a time dependent force model in damBreak3D case. I tried using setExternalVector as is defined in womersley example. The code is compiling successfully but after the force is modified, the simulation crashes. Here is my code:



Language: C++
#include "palabos3D.h" #include "palabos3D.hh"   using namespace plb; typedef double T; #define DESCRIPTOR descriptors::ForcedD3Q19Descriptor #define DYNAMICS GuoExternalForceBGKdynamics<T, DESCRIPTOR>(omega) const T pi = (T)4.0*std::atan((T)1.0); T timeDepForce(T t, T A, T omega) { return A * std::cos(omega * t); }   // Smagorinsky constant for LES model. const T cSmago = 0.14;   // Physical dimensions of the system (in meters). const T lx = 3.22; const T ly = 1.0; const T lz = 1.0;   const T rhoEmpty = T(1);   plint writeImagesIter = 100; plint getStatisticsIter = 20;   plint maxIter; plint N; plint nx, ny, nz; T delta_t, delta_x; Array<T,DESCRIPTOR<T>::d> force(1,0,0); T nuPhys, nuLB, tau, omega, Bo, surfaceTensionLB, contactAngle;   std::string outDir; plint obstacleCenterXYplane, obstacleLength=0, obstacleWidth=0, obstacleHeight=0, beginWaterReservoir, waterReservoirHeight;   plint waterLevelOne, waterLevelTwo, waterLevelThree, waterLevelFour;   T calcForce(plint iT, T A, T Tp) { T omega=2*pi/Tp; return -1*A*omega*omega*std::sin(omega*iT)*delta_t*delta_t/delta_x; }   void setupParameters() { delta_x = lz / N; nx = util::roundToInt(lx / delta_x); ny = util::roundToInt(ly / delta_x); nz = util::roundToInt(lz / delta_x);   tau = (nuPhys*DESCRIPTOR<T>::invCs2*delta_t)/(delta_x*delta_x) + 0.5; omega = 1./tau; nuLB = (tau-0.5)*DESCRIPTOR<T>::cs2; // Viscosity in lattice units. surfaceTensionLB = 0;//rhoEmpty * gLB * N * N / Bo;   obstacleCenterXYplane = util::roundToInt(0.744*N);   beginWaterReservoir = 0; waterReservoirHeight = util::roundToInt(0.3*N);   waterLevelOne = util::roundToInt(0.496*N); waterLevelTwo = util::roundToInt(2.*0.496*N); waterLevelThree = util::roundToInt(3.*0.496*N); waterLevelFour = util::roundToInt((3.*0.496 + 1.150)*N); }   void prepareCase(MultiBlockLattice3D<T,DESCRIPTOR>& lattice) { force[0] = calcForce(0,5,1.74); force[1] = 0.; force[2] = -9.81*delta_t*delta_t/delta_x; // externalForce = Array<T,3>(0,0,0); setExternalVector(lattice,lattice.getBoundingBox(),DESCRIPTOR<T>::ExternalField::forceBeginsAt,force); }   // Specifies the initial condition for the fluid (each cell is assigned the // flag "fluid", "empty", or "wall"). int initialFluidFlags(plint iX, plint iY, plint iZ) { // Place an obstacle on the left end, which is hit by the fluid. bool insideObstacle = iX >= obstacleCenterXYplane-obstacleWidth/2 && iX <= obstacleCenterXYplane+obstacleWidth/2 && iY >= ny/2-obstacleLength/2 && iY <= ny/2+obstacleLength/2 && iZ <= obstacleHeight+1;   if (insideObstacle) { return freeSurfaceFlag::wall; } else if (iX >= beginWaterReservoir && iZ <= waterReservoirHeight) { return freeSurfaceFlag::fluid; } else { return freeSurfaceFlag::empty; } }   void writeResults(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, MultiScalarField3D<T>& volumeFraction, plint iT) { static const plint nx = lattice.getNx(); static const plint ny = lattice.getNy(); static const plint nz = lattice.getNz(); Box3D slice(0, nx-1, ny/2, ny/2, 0, nz-1); ImageWriter<T> imageWriter("leeloo"); imageWriter.writeScaledPpm(createFileName("u", iT, 6), *computeVelocityNorm(lattice, slice));   imageWriter.writeScaledPpm(createFileName("rho", iT, 6), *computeDensity(lattice, slice));   imageWriter.writeScaledPpm(createFileName("volumeFraction", iT, 6), *extractSubDomain(volumeFraction, slice));   // Use a marching-cube algorithm to reconstruct the free surface and write an STL file. std::vector<T> isoLevels; isoLevels.push_back((T) 0.5); typedef TriangleSet<T>::Triangle Triangle; std::vector<Triangle> triangles; isoSurfaceMarchingCube(triangles, volumeFraction, isoLevels, volumeFraction.getBoundingBox()); TriangleSet<T>(triangles).writeBinarySTL(createFileName(outDir+"/interface", iT, 6)+".stl");   VtkImageOutput3D<T> vtkOut(createFileName("volumeFraction", iT, 6), 1.); vtkOut.writeData<float>(volumeFraction, "vf", 1.); }   void writeStatistics(FreeSurfaceFields3D<T,DESCRIPTOR>& fields) { pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << std::endl; T averageMass = freeSurfaceAverageMass<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox()); pcout << "Average Mass: " << averageMass << std::endl; T averageDensity = freeSurfaceAverageDensity<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox()); pcout << "Average Density: " << std::setprecision(12) << averageDensity << std::endl;   T averageVolumeFraction = freeSurfaceAverageVolumeFraction<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox()); pcout << "Average Volume-Fraction: " << std::setprecision(12) << averageVolumeFraction << std::endl;   pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << std::endl; }     int main(int argc, char **argv) { plbInit(&argc, &argv); global::directories().setInputDir("./");   if (global::argc() != 8) { pcout << "Error missing some input parameter\n"; }   try { global::argv(1).read(outDir); global::directories().setOutputDir(outDir+"/");   global::argv(2).read(nuPhys); global::argv(3).read(Bo); global::argv(4).read(contactAngle); global::argv(5).read(N); global::argv(6).read(delta_t); global::argv(7).read(maxIter); } catch(PlbIOException& except) { pcout << except.what() << std::endl; pcout << "The parameters for this program are :\n"; pcout << "1. Output directory name.\n"; pcout << "2. kinematic viscosity in physical Units (m^2/s) .\n"; pcout << "3. Bond number (Bo = rho * g * L^2 / gamma).\n"; pcout << "4. Contact angle (in degrees).\n"; pcout << "5. number of lattice nodes for lz .\n"; pcout << "6. delta_t .\n"; pcout << "7. maxIter .\n"; pcout << "Reasonable parameters on a desktop computer are: " << (std::string)global::argv(0) << " tmp 1.e-5 100 80.0 40 1.e-3 80000\n"; pcout << "Reasonable parameters on a parallel machine are: " << (std::string)global::argv(0) << " tmp 1.e-6 100 80.0 100 1.e-4 80000\n"; exit (EXIT_FAILURE); }     setupParameters(); pcout << "delta_t= " << delta_t << std::endl; pcout << "delta_x= " << delta_x << std::endl; pcout << "delta_t*delta_t/delta_x= " << delta_t*delta_t/delta_x << std::endl; pcout << "externalForce= " << force[2] << std::endl; pcout << "relaxation time= " << tau << std::endl; pcout << "omega= " << omega << std::endl; pcout << "kinematic viscosity physical units = " << nuPhys << std::endl; pcout << "kinematic viscosity lattice units= " << nuLB << std::endl;   global::timer("initialization").start();     SparseBlockStructure3D blockStructure(createRegularDistribution3D(nx, ny, nz));   Dynamics<T,DESCRIPTOR>* dynamics = new SmagorinskyBGKdynamics<T,DESCRIPTOR>(omega, cSmago);   // If surfaceTensionLB is 0, then the surface tension algorithm is deactivated. // If contactAngle is less than 0, then the contact angle algorithm is deactivated. FreeSurfaceFields3D<T,DESCRIPTOR> fields( blockStructure, dynamics->clone(), rhoEmpty, surfaceTensionLB, contactAngle, force ); //integrateProcessingFunctional(new ShortenBounceBack3D<T,DESCRIPTOR>, fields.lattice.getBoundingBox(), fields.freeSurfaceArgs, 0);   // Set all outer-wall cells to "wall" (here, bulk-cells are also set to "wall", but it // doesn';t matter, because they are overwritten on the next line). setToConstant(fields.flag, fields.flag.getBoundingBox(), (int)freeSurfaceFlag::wall); // In the bulk (all except outer wall layer), initialize the flags as specified by // the function "initialFluidFlags". setToFunction(fields.flag, fields.flag.getBoundingBox().enlarge(-1), initialFluidFlags); prepareCase(fields.lattice); fields.lattice.periodicity().toggle(0,true); fields.defaultInitialize();   pcout << "Time spent for setting up lattices: " << global::timer("initialization").stop() << std::endl; T lastIterationTime = T();   for (plint iT = 0; iT <= maxIter; ++iT) {   if (iT % writeImagesIter==0) { Array<T,DESCRIPTOR<T>::d> force(calcForce(iT,61,1.74),0,0); setExternalVector(fields.lattice,fields.lattice.getBoundingBox(), DESCRIPTOR<T>::ExternalField::forceBeginsAt,force); } global::timer("iteration").restart();   T sum_of_mass_matrix = T(); T lost_mass = T(); if (iT % getStatisticsIter==0) { pcout << std::endl; pcout << "ITERATION = " << iT << std::endl; pcout << "Time of last iteration is " << lastIterationTime << " seconds" << std::endl; writeStatistics(fields); sum_of_mass_matrix = fields.lattice.getInternalStatistics().getSum(0); pcout << "Sum of mass matrix: " << sum_of_mass_matrix << std::endl; lost_mass = fields.lattice.getInternalStatistics().getSum(1); pcout << "Lost mass: " << lost_mass << std::endl; pcout << "Total mass: " << sum_of_mass_matrix + lost_mass << std::endl; pcout << "Interface cells: " << fields.lattice.getInternalStatistics().getIntSum(0) << std::endl; }   if (iT % writeImagesIter == 0) { global::timer("images").start(); writeResults(fields.lattice, fields.volumeFraction, iT); pcout << "Total time spent for writing images: " << global::timer("images").stop() << std::endl; }   // This includes the collision-streaming cycle, plus all free-surface operations. fields.lattice.executeInternalProcessors(); //fields.lattice.collideAndStream(); fields.lattice.evaluateStatistics(); fields.lattice.incrementTime(); //fields.lattice.collideAndStream();   lastIterationTime = global::timer("iteration").stop(); } }

Any help will be appreciated.
Thanks
Sorry, you do not have permission to post/reply in this forum.