https://linen.dev logo
Join Discord
Powered by
# programming
  • j
    made possible by schizophrenia lol
  • f
    Currently working on verifying my CHT with a richard bénard simulation (the solver is not as in a boussinesq approximation, but fully compressible) although the pressure variations are 0. This is because my solver later will need large temperature variations. When verifying I am going for a Gr = 10^3, Ra = 7.1e-2 Using Pr = 0.71, beta = 0.01 (as I simulate ideal gas, from 100 to 99 K), Cp = 1000 (stabilization), P = 1e5, M = 8.314 g/mol, rho = 1 (solution of the paramters). dynamic viscosity = 0.001. For the wanted Gr I find g = 0.1 This leads to alpha = 1.40845e-3 kappa_f = 1.40845 I use Kr = kappa_s/kappa_F = 10 kappa_s = 14.0845 I initialize the temperature field already with T(t=0) = 100 -5/6 x Then I run my simulation, I calculate the heatFlux on the CHT wall with. Finally finding a Nu = 1.5, although it should be 1.04
  • f
    The heat flux is calculated as:
    Copy code
    void Foam::chtPimpleApp::Manager::write() const
{
    Region_FLUID::Storage& storage = regions().region_FLUID().storage();


    const volScalarField& T     = storage.thermo().T();
    const volScalarField& kappa = storage.Kappa();

    forAll(T.boundaryField(), patchI)
    {
        const fvPatchScalarField& Tpatch = T.boundaryField()[patchI];

        if (isType<solidWallMixedTemperatureCoupledFvPatchScalarField>(Tpatch))
        {
            const fvPatchScalarField& Tpatch = T.boundaryField()[patchI];

            // Thermal conductivity at this patch
            const scalarField& kappaPatch = kappa.boundaryField()[patchI];

            // Compute per-face heat flux: q = kappa * (gradT · n) * faceArea
            scalarField heatFluxPatch = kappaPatch * Tpatch.snGrad() * Tpatch.patch().magSf();

            surfaceScalarField heatFlux
            (
                IOobject
                (
                    "heatFlux",       
                    this->time().timeName(),  
                    T.mesh(),       
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                T.mesh(),            
                dimensionedScalar("zero", dimensionSet(1,0,-3,0,0,0,0), 0.0) // initialize all faces to 0
            );

            heatFlux.boundaryField()[patchI] = heatFluxPatch;

            heatFlux.write();
        }
    }
}
  • t
    I learned today that the "normal" optimized version of OF (
    WM_COMPILE_OPTION=Opt
    ) allows access beyond end of 
    UList
    , but the debug version (
    WM_COMPILE_OPTION=Debug
    ) does not and raises error like the one below. Does anyone know why the optimized UList behaves differently than the debug version? This means that one should always try to compile and and run on both Opt and Debug versions to reveal hidden bugs.
    Copy code
    --> FOAM FATAL ERROR: 
index 78000 out of range 0 ... 76799

    From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = Foam::Vector<double>; Foam::label = int]
  • s
    Bounds checking is only performed for const-qualified objects
  • t
    So, why is the behavior different in Opt vs. Debug builds?
  • s
    I would assume there is a bounds checking when the Debug flag is active
  • s
    youll have to deep dive in the code my friendo
  • s
    otherwise use the classic List<>
  • s
    if you want bounds checkking always
  • t
    Ok thanks. Unfortunately, the choice is not mine, this happens for 
    Foam::fvMesh::Cf()
    and many more of base fvMesh functions when I try to access the values on patch faces..
  • t
    The weird thing is the code works just fine with Opt build. It's just that Debug build is picky about it for some reason. I wonder if this could be somehow "working as intended"..
  • o
    just to comment about this, i had the time to test out this approach, and it is not good enought, and the path given is a little bit problematic, as i wanted to have a file on the root of the folder, and with the division of the files on the processors or not in the processors (i mean if run in parallel or in serial) the relative path will be different... and the absolute one will be case dependant so I will do something similar, but instead that will read the values and sed it from the file during the allRun, this looks more robust, simply that is not general and need to change the sed command in different simulations..
  • t
    Does anyone know how to get list of faces of a given cell from 
    fvMesh
    ? Feels like this should be a trivial thing, but I haven't found a way.. CFD-online says 
    mesh.cellFaces()
    but there is no such function. Looks like OF does not need this information, so it is not available. Have to build it yourself.
  • t
  • s
    it is cellFaces()
  • s
    that will give you a list of lists
  • t
    Recent .com or .org does not seem recognize that.. But no biggie, I'll just do it myself. 🤸
  • s
    Copy code
    const labelList& faceOwner = mesh.faceOwner();
  • t
    Indeed you are correct. checkIndex() is called ONLY when FULLDEBUG is true. So, it is a good idea to try compiling and running stuff with OpenFOAM built with 
    WM_COMPILE_OPTION=Debug
    . And sometimes run through valgrind, also.
    Copy code
    template<class T>
inline T& Foam::UList<T>::operator[](const label i)
{
    #ifdef FULLDEBUG
    checkIndex(i);
    #endif
    return v_[i];
}
  • s
    You thought I would lie?
  • t
    NO. What makes you say that?
  • f
    Love it when you thought your biconjugate gradient solver was actually converging compared to deflated conjugate gradient but the whole point of it was that previous people coded the following:
    Copy code
    // Read solver for V
            word solverV
            (
                mesh.solutionDict().solver(V.name()).lookup("solver")
            );

            if(solverV != "deflatedPCG")
            {
// TODO: Read reference values
//                 VEqn.setReference(VRefCell, VRefValue);
                VEqn.setReference(0, 0.0);
            }`
    Which doesnt make physical sense, I do get it when you have 
    p_rgh
    but there is not real reference for it. So I commented it out ran the biconjugate gradient again and no convergence... But some interesting insights
  • f
  • f
  • f
    Some diverging point on the surface
  • t
    mesh there has some weird cells?
  • f
    Not specifically if i recall but it was not converging anyways if i dont set ref value. But indeed must be mesh, even when I made it coarse/fine no difference. When changed coil thickness not yet a solution. When adding a buffer around it to snap on (currently working in that).
  • f
    Do you have any clue why setReference would lead to convergence although to what it converges is ultimately not physical. I think it just takes some reference cell to compute residual but that would generally be very weird. And therefore although said it is converged it technically is not, by means of a reference that doesnt make sense, although what the reference truely means I am not familiar with.
  • t
    IO question: I need to add a feature to smoothMesh which would write and read a bunch of `boolList`s which contain classification of mesh points. The list length is mesh.nPoints(). Reading will be done at first run iteration, writing at writeInterval and last iteration. IOList? Instead of boolList could be some kind of Field. IOField? How to do this? At least not IOList, it does not seem to exist in org. Maybe try to use labelIOList with values 0 and 1, that exists in both org and com 🤦