Hello Dumux,
I tried to calculate the fluxes over a certain part of the boundary of the
model domain. For this purpose I used an ElementIterator to iterate over the
elements of the grid. In each element I iterate over the intersections and in
each intersection over the vertices. With information from the
IntersectionIterator and the vertex index I can get the boundary face index and
the boundary face position. I test whether I am at the desired location on the
grid and then calculate the local residual of the corresponding element. This
works just fine.
The idea was to use the computeFlux-function of the local residual afterwards
with the corresponding boundary face index and onBoundary=true to get the flux
over the boundary. However, this will not work because the fluidState is not
defined when the local fluxVariables are constructed and this will cause
Valgrind to complain.
To me it came as a surprise that it is not possible to use the
computeFlux-function outside the localresidual since this is quite an important
quantity to access when traversing the grid. I see that several people have
tried just the same approach but commented the code afterwards, probably for
the reasons mentioned above. The workaround which was used instead is to
locally construct the fluxVariables of the element and then do just the same
flux calculation as in the localresidual. This way, one is forced to copy code
and have code which does exactly the same in two different locations. I think
this should be improved and the fluxes over the interfaces should be readily
available when iterating over the grid.
The solution to this problem would be to add a second computeFlux-function to
the local residual. This function would get all the parameters needed for the
flux calculation without the need to touch the private functions of the local
residual. This would be a public version of computeFlux (the current
computeFlux function is private in nature since it can't be called outside the
local residual without causing problems and maybe should be set private for
consistency).
I believe the problem is the same for all (implicit) models. I attached my code
which is specific to my multidomain/fuel cell model. The first function is
called in the multidomain problem. The computeFlux-Function is part of the
localresidual.
These are just some ideas. Please let me know what you think.
Best regards
Georg Futter
--------------------------
German Aerospace Center (DLR)
Institute of Engineering Thermodynamics | Computational Electrochemistry |
Pfaffenwaldring 38-40 | 70569 Stuttgart
Dipl.-Ing. Georg Futter | Ph.D. student
Telefon 0711/6862-8135 | [email protected]<mailto:[email protected]>
www.DLR.de<http://www.dlr.de/>
/*!
* \brief Calculates the electrical current density per cell area
* from the fluxes at the anode current collector [A/m^2].
*/
void calculateAnodeCurrentDensityFromFlux()
{
anodeCurrentDensityFromFlux_ = 0.0;
FVElementGeometry1 fvGeometry1;
ElementVolumeVariables1 elemVolVars1;
PrimaryVariables1 flux(0.0);
ElementIterator1 eIt1 = this->sdGridView1().template begin<0>();
ElementIterator1 eEndIt1 = this->sdGridView1().template end<0>();
for (; eIt1 != eEndIt1; ++eIt1)
{
fvGeometry1.update(this->sdGridView1(), *eIt1);
elemVolVars1.update(this->sdProblem1(),
*eIt1,
fvGeometry1,
/*oldSol=*/false);
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
typedef Dune::ReferenceElement<Scalar, dim> ReferenceElement;
const ReferenceElement &refElement =
ReferenceElements::general(eIt1->geometry().type());
IntersectionIterator1 isIt1 = this->sdGridView1().ibegin(*eIt1);
const IntersectionIterator1 &endIt1 =
this->sdGridView1().iend(*eIt1);
for (; isIt1 != endIt1; ++isIt1)
{
// handle only intersections on the boundary
if (!isIt1->boundary())
continue;
// assemble the boundary for all vertices of the current face
const int fIdx = isIt1->indexInInside();
const int numFaceVertices = refElement.size(fIdx, 1, dim);
// loop over the single vertices on the current face
for (int faceVertexIdx = 0; faceVertexIdx < numFaceVertices;
++faceVertexIdx)
{
const int boundaryFaceIdx =
fvGeometry1.boundaryFaceIndex(fIdx, faceVertexIdx);
GlobalPosition globalPos =
fvGeometry1.boundaryFace[boundaryFaceIdx].ipGlobal;
if(fcGeometry_.atAnodeCurrentCollectorTest(globalPos))
{
this->localResidual1().computeFlux(flux,
this->sdProblem1(),
*eIt1,
fvGeometry1,
boundaryFaceIdx,
elemVolVars1,
/*onBoundary=*/
true);
anodeCurrentDensityFromFlux_ += flux[elecChargeEqIdx1];
// A
}
} // loop over vertices
} // loop over intersections
} // loop over elements
anodeCurrentDensityFromFlux_ /= fcGeometry_.cellArea(); // --> A/m^2
}
/*!
* \brief Evaluates the total flux of all conservation quantities
* over a face of a subcontrol volume.
*
* \param flux The flux over the SCV (sub-control-volume) face for each
component
* \param problem The problem
* \param element The element
* \param fvGeometry The geometry of the element
* \param fIdx The index of the SCV face
* \param elemVolVarsCur The current element volume variables
* \param onBoundary A boolean variable to specify whether the flux
variables
* are calculated for interior SCV faces or boundary faces,
default=false
*/
void computeFlux(PrimaryVariables &flux,
const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const unsigned int fIdx,
const ElementVolumeVariables &elemVolVarsCur,
const bool onBoundary=false) const
{
FluxVariables fluxVars(problem,
element,
fvGeometry,
fIdx,
elemVolVarsCur,
onBoundary);
flux = 0.0;
MassResid::computeFlux(flux, fluxVars, elemVolVarsCur);
Valgrind::CheckDefined(flux);
/*
* EnergyResid also called in the MassResid
* 1) Makes some sense because energy is also carried by mass
* 2) The component-wise mass flux in each phase is needed.
*/
const GlobalPosition globalPos =
fvGeometry.subContVolFace[fIdx].ipGlobal;
ChargeResid::computeFlux(flux, fluxVars, elemVolVarsCur, globalPos);
Valgrind::CheckDefined(flux);
}
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
