GAS_AdaptiveViscosity.h

Go to the documentation of this file.

#ifndef _GAS_AdaptiveViscosity_
#define _GAS_AdaptiveViscosity_

#include "GAS_SubSolver.h"
#include "GAS_Utils.h"

#include <SIM/SIM_RawField.h>
#include <SIM/SIM_RawIndexField.h>
#include <SIM/SIM_ScalarField.h>
#include <SIM/SIM_VectorField.h>

#include <UT/UT_FaceGradedOctreeLabels.h>
#include <UT/UT_SparseMatrix.h>
#include <UT/UT_VoxelArray.h>


class GAS_OctreeFaceVectorFieldInterpolator
{
    enum NodeLabels
    {
        INACTIVE_NODE,
        ACTIVE_NODE,
        DEPENDENT_NODE
    };

    using OctreeFaceLabels = UT_FaceGradedOctreeLabels::OctreeFaceLabels;
    using OctreeCellLabels = UT_FaceGradedOctreeLabels::OctreeCellLabels;

public:

    // The valid faces can in one of three states:
    // 1. active --> the face is present in the octree AND inside the simulation domain.
    //                  Represented with an OctreeFaceLabel of ACTIVE
    // 2. inactive -> the face is present in the octree but outside of the simulation domain
    //                  Represented with an OctreeFaceLabel of INACTIVE
    // 3. unassigned -> the face is not present in the octree.
    //                  Represented with an OctreeFaceLabel of UNASSIGNED

    GAS_OctreeFaceVectorFieldInterpolator(const UT_FaceGradedOctreeLabels &octree_labels,
                                            const UT_Array<UT_Array<SIM_RawField>> &staggered_velocity,
                                            const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices);

    fpreal getValue(const UT_Vector3 &sample_point, const int axis) const;

private:

    void findOccupiedNodeTiles(UT_Array<bool> &occupied_tile_list,
                                const UT_VoxelArray<int> &octree_labels,
                                const int level) const;

    void setActiveNodes(const int level);

    void sampleActiveNodes(UT_Array<SIM_RawField> &node_weights,
                            SIM_RawIndexField &node_flags,
                            const int level);

    void bubbleActiveNodeValues(UT_Array<UT_Array<SIM_RawField>> &node_weights,
                                UT_Array<SIM_RawIndexField> &node_flags,
                                const int level);

    void finishIncompleteNodes(UT_Array<UT_Array<SIM_RawField>> &node_weights,
                                UT_Array<SIM_RawIndexField> &node_flags,
                                const int level);

    void normalizeActiveNodes(UT_Array<SIM_RawField> &node_weights,
                                UT_Array<SIM_RawIndexField> &node_flags,
                                const int level);

    void distributeNodeValuesDown(const int level);

    UT_Array<SIM_RawField> m_node_labels;
    UT_Array<UT_Array<SIM_RawField>> m_node_values;

    const UT_FaceGradedOctreeLabels &m_octree_labels;
    const UT_Array<UT_Array<SIM_RawField>> &m_velocity;
    const UT_Array<UT_Array<SIM_RawIndexField>> &m_velocity_indices;
};

class GAS_API GAS_AdaptiveViscosity : public GAS_SubSolver
{
    using OctreeCellLabels = UT_FaceGradedOctreeLabels::OctreeCellLabels;

public:

    struct StressStencilFace
    {
        StressStencilFace(exint face_index, fpreal coefficient)
            : m_face_index(face_index)
            , m_coefficient(coefficient)
            {}

        exint m_face_index;
        fpreal m_coefficient;
    };

    static constexpr exint UNLABELED_FACE = -1;
    static constexpr exint OUTSIDE_FACE = -2;
    static constexpr exint COLLISION_FACE = -3;
    static constexpr exint FLUID_FACE = -4;
    static constexpr exint BOUNDARY_FACE = -5;

    static constexpr exint OUTSIDE_STRESS = -1;
    static constexpr exint ACTIVE_STRESS = -2;
    static constexpr exint INACTIVE_STRESS = -3;

    GET_DATA_FUNC_F(SIM_NAME_SCALE, ViscosityScale);
    GET_DATA_FUNC_F("minviscosity", MinViscosity);

    GET_DATA_FUNC_F(SIM_NAME_TOLERANCE, SolverTolerance);
    GET_DATA_FUNC_I("maxiterations", MaxIterations);
    GET_DATA_FUNC_I("numbersupersamples", NumberSuperSamples);

    GET_DATA_FUNC_I("floatprecision", FloatPrecision);

    GET_DATA_FUNC_I("octreelevels", OctreeLevels);
    GET_DATA_FUNC_I("finebandwidth", FineBandwidth);

    GET_DATA_FUNC_B("applysolidweights", DoApplySolidWeights);
    
    GET_DATA_FUNC_B("generateoctree", GenerateOctree);
    GET_DATA_FUNC_B("onlygenerateoctree", OnlyGenerateOctree);

    GET_DATA_FUNC_F("waterline",Waterline);
    GET_DATA_FUNC_B("usewaterline",UseWaterline);
    GET_DATA_FUNC_V3("waterlinedirection",WaterlineDirection);

    GET_DATA_FUNC_F("extrapolation", Extrapolation);

    GET_DATA_FUNC_B(GAS_NAME_USETIMESTEP, UseTimeStep);
    GET_DATA_FUNC_F(GAS_NAME_TIMESCALE, TimeScale);

protected:
    explicit GAS_AdaptiveViscosity(const SIM_DataFactory *factory);
            ~GAS_AdaptiveViscosity() override;

    // The overloaded callback that GAS_SubSolver will invoke to
    // perform our actual computation.  We are giving a single object
    // at a time to work on.
    
    bool
    solveGasSubclass(SIM_Engine &engine, SIM_Object *obj,
                     SIM_Time time, SIM_Time timestep) override;

    void initializeSubclass() override;

private:
    // We define this to be a DOP_Auto node which means we do not
    // need to implement a DOP_Node derivative for this data.  Instead,
    // this description is used to define the interface.
    static const SIM_DopDescription *getDopDescription();
    /// These macros are necessary to bind our node to the factory and
    /// ensure useful constants like BaseClass are defined.
    DECLARE_STANDARD_GETCASTTOTYPE();
    DECLARE_DATAFACTORY(GAS_AdaptiveViscosity, GAS_SubSolver,
                        "Gas Adaptive Viscosity", getDopDescription());

    void
    buildIntegrationWeights(UT_Array<SIM_RawField> &face_volume_weights,
                            SIM_RawField &center_volume_weights,
                            UT_Array<SIM_RawField> &edge_volume_weights,
                            const SIM_RawField &surface,
                            const SIM_RawField &collision,
                            const SIM_BoundaryLine &world_boundary_line,
                            const SIM_VectorField &velocity,
                            const fpreal extrapolation,
                            const bool do_apply_collision_weights) const;

    //
    // Helper methods to build the layer of fine cells along the liquid boundary
    //

    // TODO: consider just putting these as a function inside the C file?

    void
    setFineCellMaskInterior(UT_Array<UT_Array<UT_Vector3I>> &parallel_boundary_cell_list,
                            UT_VoxelArrayF &fine_cell_mask,
                            const SIM_RawField &surface,
                            const SIM_RawField &collision,
                            const SIM_RawField &center_volume_weights,
                            const fpreal extrapolation) const;

    void
    findOccupiedFineTiles(UT_Array<bool> &is_tile_occupied_list,
                            const UT_Array<UT_Vector3I> &fine_cell_list,
                            const UT_VoxelArrayF &fine_cell_mask) const;

    void
    uncompressFineTiles(UT_VoxelArrayF &fine_cell_mask,
                        const UT_Array<bool> &is_tile_occupied_list) const;

    void
    setFineCells(UT_VoxelArrayF &fine_cell_mask,
                    const UT_Array<UT_Vector3I> &fine_cell_list) const;

    void
    loadNextOutsideFineCellLayer(UT_Array<UT_Array<UT_Vector3I>> &parallel_next_cell_list,
                                    const UT_Array<UT_Vector3I> &old_cell_list,
                                    const UT_VoxelArrayF &fine_cell_mask,
                                    const SIM_RawField &center_volume_weights) const;

    void
    loadNextInsideFineCellLayer(UT_Array<UT_Array<UT_Vector3I>> &parallel_next_cell_list,
                                const UT_Array<UT_Vector3I> &old_cell_list,
                                const UT_VoxelArrayF &fine_cell_mask,
                                const SIM_RawField &center_volume_weights) const;

    void
    setCustomFineCells(UT_VoxelArrayF &fine_cell_mask,
                        const SIM_RawField &custom_fine_cells) const;

    void
    buildOctreeCellLabels(UT_FaceGradedOctreeLabels &octree_labels,
                            const SIM_RawField &surface,
                            const SIM_RawField &collision,
                            const SIM_RawField &center_volume_weights,
                            const SIM_RawField *custom_fine_cells,
                            const int octree_levels,
                            const fpreal extrapolation,
                            const int inside_fine_cell_width,
                            const int outside_fine_cell_width) const;

    //
    // Helper methods for building regular grid velocity labels
    //

    void
    findOccupiedRegularVelocityTiles(UT_Array<bool> &is_tile_occupied_list,
                                        const SIM_RawField &surface,
                                        const SIM_RawIndexField &regular_velocity_labels,
                                        const int axis) const;

    void
    classifyRegularVelocityFaces(SIM_RawIndexField &regular_velocity_labels,
                                    const SIM_RawField &surface,
                                    const SIM_RawField &collision,
                                    const SIM_RawField &center_volume_weights,
                                    const UT_Array<SIM_RawField> &edge_volume_weights,
                                    const SIM_RawField &face_volume_weights,
                                    const SIM_VectorField &velocity,
                                    const SIM_BoundaryLine &index_boundary_line,
                                    const int axis,
                                    const fpreal extrapolation) const;

    void
    buildRegularVelocityLabels(UT_Array<SIM_RawIndexField> &regular_velocity_labels,
                                const SIM_RawField &surface,
                                const SIM_RawField &collision,
                                const SIM_RawField &center_volume_weights,
                                const UT_Array<SIM_RawField> &edge_volume_weights,
                                const UT_Array<SIM_RawField> &face_volume_weights,
                                const SIM_VectorField &velocity,
                                const SIM_BoundaryLine &world_boundary_line,
                                const fpreal extrapolation) const;

    //
    // Helper methods for building octree velocity indices
    //

    void
    findOccupiedOctreeVelocityTiles(UT_Array<bool> &is_tile_occupied_list,
                                    const UT_VoxelArray<int> &octree_labels,
                                    const SIM_RawIndexField &octree_velocity_indices,
                                    const int axis) const;

    void
    classifyOctreeVelocityFaces(SIM_RawIndexField &octree_velocity_indices,
                                const UT_FaceGradedOctreeLabels &octree_labels,
                                const SIM_RawField &surface,
                                const SIM_RawField &collision,
                                const SIM_RawField &center_volume_weights,
                                const UT_Array<SIM_RawField> &edge_volume_weights,
                                const SIM_RawField &face_volume_weights,
                                const SIM_VectorField &velocity,
                                const SIM_BoundaryLine &world_boundary_line,
                                const int axis, const int level,
                                const fpreal extrapolation) const;

    exint
    buildOctreeVelocityIndices(UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                                const SIM_RawField &surface,
                                const SIM_RawField &collision,
                                const UT_FaceGradedOctreeLabels &octree_labels,
                                const SIM_RawField &center_volume_weights,
                                const UT_Array<SIM_RawField> &edge_volume_weights,
                                const UT_Array<SIM_RawField> &face_volume_weights,
                                const SIM_VectorField &velocity,
                                const SIM_BoundaryLine &world_boundary_line,
                                const fpreal extrapolation) const;

    //
    // Helper methods for building edge stress stencil indices
    //

    void
    findOccupiedEdgeStressTiles(UT_Array<bool> &is_tile_occupied_list,
                                const UT_VoxelArray<int> &octree_labels,
                                const SIM_RawIndexField &edge_stress_indices,
                                const int axis) const;

    void
    classifyEdgeStresses(SIM_RawIndexField &edge_stress_indices,
                            const UT_FaceGradedOctreeLabels &octree_labels,
                            const SIM_RawField &surface,
                            const SIM_RawField &edge_volume_weights,
                            const int axis, const int level) const;

    exint
    buildEdgeStressIndices(UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                            const SIM_RawField &surface,
                            const UT_FaceGradedOctreeLabels &octree_labels,
                            const UT_Array<SIM_RawField> &edge_volume_weights) const;

    //
    // Helper methods for building center stress stencils
    //
    
    void
    classifyCenterStresses(SIM_RawIndexField &center_stress_indices,
                            const UT_VoxelArray<int> &octree_labels,
                            const SIM_RawField &center_volume_weights,
                            const int level) const;

    exint
    buildCenterStressIndices(UT_Array<SIM_RawIndexField> &center_stress_indices,
                                const UT_FaceGradedOctreeLabels &octree_labels,
                                const SIM_RawField &center_volume_weights) const;

    //
    // Help methods to build edge stress stencils
    //

    void
    getEdgeStressFaces(UT_Array<GAS_AdaptiveViscosity::StressStencilFace> &stencil_faces,
                        UT_Array<fpreal> &boundary_faces,
                        const UT_Vector3I &edge,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                        const UT_FaceGradedOctreeLabels &octree_labels,
                        const SIM_VectorField *collision_velocity,
                        const int axis, const int level) const;

    fpreal
    edgeOctreeVolumes(const UT_Vector3I &edge,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                        const UT_FaceGradedOctreeLabels &octree_labels,
                        const int axis, const int level) const;

    struct EdgeStressStencilParms
    {
        EdgeStressStencilParms(const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                                const UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                                const UT_FaceGradedOctreeLabels &octree_labels,
                                const SIM_VectorField *collision_velocity,
                                const UT_Array<SIM_RawField> &edge_volume_weights,
                                const SIM_RawField &viscosity,
                                const fpreal viscosity_scale,
                                const fpreal min_viscosity,
                                const fpreal dt)
            : m_octree_velocity_indices(octree_velocity_indices),
              m_edge_stress_indices(edge_stress_indices),
              m_octree_labels(octree_labels),
              m_collision_velocity(collision_velocity),
              m_edge_volume_weights(edge_volume_weights),
              m_viscosity(viscosity),
              m_viscosity_scale(viscosity_scale),
              m_min_viscosity(min_viscosity),
              m_dt(dt)
        {
        }

        const UT_Array<UT_Array<SIM_RawIndexField>> &m_octree_velocity_indices;
        const UT_Array<UT_Array<SIM_RawIndexField>> &m_edge_stress_indices;
        const UT_FaceGradedOctreeLabels &m_octree_labels;
        const SIM_VectorField *m_collision_velocity;
        const UT_Array<SIM_RawField> &m_edge_volume_weights;
        const SIM_RawField &m_viscosity;
        const fpreal m_viscosity_scale;
        const fpreal m_min_viscosity;
        const fpreal m_dt;
    };

    void
    buildEdgeStressStencils(UT_Array<UT_Array<StressStencilFace>> &edge_stress_stencils,
                            UT_Array<UT_Array<fpreal>> &edge_stress_boundary_stencils,
                            UT_Array<fpreal> &edge_stress_weights,
                            const EdgeStressStencilParms &parms,
                            const int axis, const int level) const;

    //
    // Helper methods to build cell centered stress stencils
    //

    void
    getCenterStressFaces(UT_Array<GAS_AdaptiveViscosity::StressStencilFace> &stencil_faces,
                            UT_Array<fpreal> &boundary_faces,
                            const UT_Vector3I &cell,
                            const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                            const UT_FaceGradedOctreeLabels &octree_labels,
                            const SIM_VectorField *collision_velocity,
                            const int axis, const int level) const;

    struct CenterStressStencilParms
    {
        CenterStressStencilParms(const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                                    const UT_Array<SIM_RawIndexField> &center_stress_indices,
                                    const UT_FaceGradedOctreeLabels &octree_labels,
                                    const SIM_VectorField *collision_velocity)
            : m_octree_velocity_indices(octree_velocity_indices),
              m_center_stress_indices(center_stress_indices),
              m_octree_labels(octree_labels),
              m_collision_velocity(collision_velocity)
        {
        }

        const UT_Array<UT_Array<SIM_RawIndexField>> &m_octree_velocity_indices;
        const UT_Array<SIM_RawIndexField> &m_center_stress_indices;
        const UT_FaceGradedOctreeLabels &m_octree_labels;
        const SIM_VectorField *m_collision_velocity;
    };

    void
    buildCenterStressStencils(UT_Array<UT_Array<StressStencilFace>> &center_stress_stencils,
                                UT_Array<UT_Array<fpreal>> &center_stress_boundary_stencils,
                                const CenterStressStencilParms &parms,
                                const int axis, const int level,
                                const exint center_stress_count) const;

    void
    buildCenterStressWeights(UT_Array<fpreal> &center_stress_weights,
                                const SIM_RawIndexField &center_stress_indices,
                                const SIM_RawField &center_volume_weights,
                                const SIM_RawField &viscosity,
                                const fpreal viscosity_scale,
                                const fpreal min_viscosity,
                                const fpreal dt,
                                const int level) const;

    //
    // Helper method to update regular grid velocities with interpolated solution from the octree
    //

    void
    applyVelocitiesToRegularGrid(SIM_RawField &regular_velocity,
                                    const GAS_OctreeFaceVectorFieldInterpolator &interpolator,
                                    const SIM_RawIndexField &octree_velocity_indices,
                                    const SIM_RawIndexField &regular_velocity_labels,
                                    const SIM_RawField &collision,
                                    const SIM_VectorField *collision_velocity,
                                    const UT_Array<SIM_RawField> &octree_velocity,
                                    const int axis) const;

    //
    // Helper methods to unit test adaptivity
    //

    bool
    octreeVelocityGradingUnitTest(const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                                    const UT_FaceGradedOctreeLabels &octree_labels) const;

    bool
    edgeStressUnitTest(const UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                        const SIM_RawField &surface,
                        const UT_FaceGradedOctreeLabels &octree_labels) const;

    bool
    centerStresUnitTest(const UT_Array<SIM_RawIndexField> &center_stress_indices,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &edge_stress_indices,
                        const UT_Array<UT_Array<SIM_RawIndexField>> &octree_velocity_indices,
                        const UT_FaceGradedOctreeLabels &octree_labels) const;   
};

#endif