Hi Ranjeet,

to me the problem seems to be that you do not specify Dirichlet boundary conditions anywhere. You must have a Dirichlet boundary somewhere for the system to be defined properly and that is why it works when you use the original boundaryTypes() function.

Best wishes,
Dennis


On 18.09.2018 07:58, Ranjeet Singh wrote:
Hi, I am implementing the 2P PorousmediumFlowProblem with cornerpoint spe9 grid. Boundary conditions are no flow along all boundary and injection source at one corner (0th element) and production source at another opposite corner. Boundary and Initial condition are set as follows:

//problem.hh
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
    {
        BoundaryTypes bcTypes;
        bcTypes.setAllNeumann();
        return bcTypes;
    }

    NumEqVector source(const Element &element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolume &scv) const
    {
        NumEqVector values(0.0);

        int eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
        if (eIdx == injectionElement_)
            values[FluidSystem::phase1Idx] = injectionRate_/element.geometry().volume();

        if (eIdx == productionElement_)
            values[FluidSystem::phase0Idx] = -productionRate_/element.geometry().volume();

        return values;
    }

    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    {
        PrimaryVariables values;

        // hydrostatic pressure
        Scalar densityW = 1000;
        values[Indices::pressureIdx] = 5e5 + densityW*(this->gravity()*globalPos);  // inital pw
        values[Indices::saturationIdx] = 0.0; // initial sn

        return values;
    }
I guess I am doing something wrong with boundary condition. In the cornerpoint test problem, the boundary condition is implemented as,
    BoundaryTypes boundaryTypes(const Element &element,
                                const SubControlVolumeFace &scvf) const
    {
        BoundaryTypes bcTypes;
        // set no-flux on top and bottom, hydrostatic on the rest
        // use the scvf normal to decide
        const auto& normal = scvf.unitOuterNormal();
        using std::abs;
        if (abs(normal[dimWorld-1]) > 0.5)
            bcTypes.setAllNeumann();
        else
            bcTypes.setAllDirichlet();

        return bcTypes;
    }
when I implement the same BC, then solution gets converged. Could you please explain what I am doing wrong and also, in abs(normal[dimwold-1])> 0.5
, why it is compared with 0.5, not 0?

Thank you!

Regards,
Ranjeet


_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux

_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux

Reply via email to