[olwi]

Hello,

I solve an equation for electric potential, epot, in an electrically conducting fluid. Then I want to calculate the electric currents, which are (for now ignoring the issue of discretisation complications) computed like this:
J = sigma * ( -grad(epot) + (U ^ B0) )

J is a volVectorField, which is written with the other fields for postprocessing.

How can I give the J field the correct values in the boundary nodes, and not only in the internal nodes? U and B0 have nice values in the boundary nodes, because they are solved for (U), or set to a constant value on initialization (B0). grad(epot), however, is probably not defined for the boundary field (gradient calculation uses epot values on cell faces, which is ok.). Or? What I'm sure of is that J is zero on all boundaries in paraview. I have a "symmetry" plane (kind of...) on which J must obtain a non-zero normal value, so this is unfortunate. (I also want the current to have its correct value on the boundary, because it should be available for other reasons later on.)

Creators of the existing openFOAM MHD solver might wonder why I do it this way (solving the epot equation, etc.). The reason is I want to compare with existing implementations of these equations in two commercial codes. I will have to ask questions on other issues as well to make this work, but that's for later...

Thanks in advance for any advice.

/Ola

[hjasak]

If you want to specify J at the boundary, you will need to create it given a list of boundary patch types. If you give nothing (this is what happens now), you will get a calculated type for the boundary patch field, which is a default choice.

OpenFOAM also has boundaries whoe type is constrained by the mesh definition, e.g. empty, processor, symmetry, cyclic and wedge. This will be handled automatically, i.e. the code won't let you specify a fixed value patch type on the symmetry boundary.

Hrv

[olwi]

Thanks, Hrv. I now set boundary patches in a "J" file located in the "0" directory, and it works... The "symmetry" plane in question is symmetric (~zeroGradient) for all variables execect the electric potential epot which is fixed value 0 (reflection symmetry).

Some follow-up questions arise...
1. With boundary patch info defined as above, WHEN are the boundary values updated? At output, or already when the assignment is done ("=" operator?).
2. If I just define a scalar field inside the code, like
volScalarField a;
Can boundary patch info still be supplied (other than through an IOobject, that is) in order for the boundary values to be updated?

On a slightly diffrent theme: Assume I want to compute my current on cell faces instead, I would define a variable
faceVectorField Jf;
This would be calculated as:
Jf = sigma_f*( (-grad(epot))_f + U_f x B0_f );
a) I assume there is a nifty interpolation operator around for sigma_f, U_f and B0_f (values on faces). Would the facial boundary field values on the RHS properly update the boundary values on the LHS in the assignment?
b) Is there a (full vector) facial gradient operator, grad()_f above, whose result is consistent with the result of snGrad()?
c) Knowing Jf, is there a smart way of reconstructing a cell-centre J vector? (The whole point is to make sure div(J)=0...)

I realize I ask to many questions... Any partial input is welcome...

/Ola

[hjasak]

Well, you've hacked it, but never mind for the moment (you should have done this in the code instead).

1) Boundary conditions

Formally, the values are updated when you call correcBoundaryConditions() member function. However, the assignment operation will happen both on the internal field and the patch fields, but some of the fvPatchFields have special versions of the operator. For example, operator= on a fixedValue field does nothing.

Also, beware of correcting boundary conditions may need other stuff to evaluate. An example would be a totalPressure b.c. which needs the existence of velocity U to evaluate itself. If you read in p (with this b.c.), and call p.correctBoundaryConditions() when U is not present, the code will fall over.

2) Defining the field

Have a look at the GeometricField constructors, the relevant one looks like this:

//- Constructor given IOobject, mesh, dimensioned<type> and patch type.
GeometricField
(
const IOobject&,
const Mesh&,
const dimensioned<type>&,
const word& patchFieldType=PatchField<type>::calculatedType()
);

which (should) answer your question :-)

The default patch field type is "calculated", which means "I don't know what this should be, so I'll give you a generic one with regular algebra. If you try and ask the calculated type for matrix coefficients, the code falls over (e.g. when I mis-spell the b.c. type on pressure in icoFoam) :-) For a calculated b.c. evaluation does nothing: values are assembled through field calculus. In other words, you don't need the evaluation because the stuff is already there.

3) Other stuff
- that would be a surfaceScalarField
- you are missing some interpolation of gradients; more likely, you want to calculate the surfacenormal gradient, called snGrad. Watch out, it's dotted with the face area vector.
- boundary conditions on surface fields are a bit tricky: in reality you can only have fixedValue or calculated and the code will enforce this to you. If you don't believe me, try figuring out what zeroGradient means for a field defined on faces
- use fvc::reconstruct - built specially for the purpose (e.g. creating cell centre velocity from the face flux). Beware of face dot-product.

Hrv

[olwi]

Thanks Hrvoje,

I know I'm hacking... I'm trying to hack a proof-of-concept in late evenings, hoping to learn more about openFOAM and C++ in the process. Once there I'll convince someone it's worthwhile doing it properly, and hopefully I'll understand how by then.

The reconstruct method looks promising to me. I wasn't sure I could reconstruct a cell centre vector, only knowing the face-normal fluxes. (That's why I wanted to compute a full surfaceVectorField (sic) first.) If you say it works, I believe you...

Your info will keep me busy another evening or two.

/Ola

[olwi]

Just some feedback to Hrvoje: It works. Well!
There's always more questions, but I put them in separate threads, letting this one slide into oblivion.

/Ola