[sverdoold]

Hi,
it took me quite some time but I currently have an electrostaticFoam case running. Now I need the computed result at points which are to be specified. I figured out that I can/should use the sample utility to do exactly this.
Now comes my problem however: I need the values of the rhoFlux field. For some reason, I only succeed in extracting the phi and/or the rho values with the sample utility.
Does anyone know why this is and whether it will be possible to access the rhoFlux values in this case?

Thanks in advance,

Sjaak

[sverdoold]

A small update on my previous problem:
I figured out that I should change the standard electrostaticFoam slightly in order to have the solver write things for rhoFlux to disk.
By changing those values, I do see files bein created containing (I think) rhoFlux stuff.

However, when I try to sample at specific points, nothing gets outputted! Does anyone have any clue what's going on and what I might be doing wrong?

Thanks!

[msrinath80]

I don't know anything about electrostaticFoam. But I do know one thing. There are no magicians here. Unless you give specifics of the error message and explain what you were trying to do, not many people have the time to ponder over your issue. Hope that helps.

[sverdoold]

ok, sorry! I'll try to be more clear.....

I am trying to determine the electric field for a simple but variable geometry. This electric field has to be used in some other software so I would like to output this field at positions that I can specify myself. I just want to output the values to a file so the other software can read the values in and start using them.

By now I have the situation that my other software can output the geometry it is working on and I can process this geometry so it is usable by OpenFoam. I can set the boundary correct boundary conditions, and the calculations seem to be going fine. So far so good.

Then I need to get the rhoFlux field (which is the physical electric field) at the points I specify. I can use the sample utility to do this, provided I give it a correct sampleDict file. This is also no problem for the rho and the phi variables. For the rhoFlux field however, nothing is outputted, and this is exactly the field that I am interested in!

When I took a closer look at the createFields.H file of the solver, I noticed that for the rhoFlux field the IOobject was set to NO_WRITE whereas for rho and phi it was set to AUTO_WRITE. Therefore I changed this property for the rhoFlux field to AUTO_WRITE as well and recompiled the solver.
After running the solver again for the same case, files concerning rhoFlux showed up in the directories that are created for each time step. However, I still could not extract any rhoFlux values using the sample utility whereas the rho and phi values are no problem at all.

Everything seems to be fine, no error messages are generated, so most likely I am overseeing something. Is there maybe some setting somewhere which determines which fields can be sample'd and which cannot?

[msrinath80]

ok, trivial question. You did mention rhoFlux at the end of sampleDict right?

E.g:
fields
(
U
p
rhoFlux
);

[sverdoold]

yes I did.
Also the output from the sample utility shows all times it is processing, but no files are outputted.
When I have rho, phi and rhoFlux at the end of sampleDict the rho and the phi are outputted, but the rhoFlux is not...

A typical output with rho,phi and rhoFlux in the sample.Dict is:

Time = 0.01
Sampling volScalarField rho
Sampling volScalarField phi
Writing fields to "/home/sjaak/OpenFOAM/sjaak-1.4.1/run/electrostaticFoam/BoundaryTest1/samples/0. 01/Potential_rho_phi.xy"

so it seems to me that rhoFlux is not being sampled at all

[msrinath80]

Check out the note by Mattijs in this post:

http://www.cfd-online.com/OpenFOAM_D...es/1/3611.html

Is your point sitting at/near a boundary?

[sverdoold]

I read the post and just to be sure, I put a couple of points in my sampleDict of which I am sure that they are not located at a geometry boundary. Besides I made sure that I am using the cloud option.
The sampleDict actually is;

FoamFile
{
version 2.0;
format ascii;

root "/home/sjaak/OpenFOAM/sjaak-1.4.1/run/electrostaticFoam";
case "BoundaryTest1";
instance "system";
local "";

class dictionary;
object sampleDict;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

arguments "/home/sjaak/OpenFOAM/sjaak-1.4.1/run/electrostaticFoam/BoundaryTest1" off off off;

interpolationScheme cell;

writeFormat raw;

sampleSets
(
cloud
{
name Potential;
axis xyz;
points
(
(0.1 0.05 0.25)
(0.2 0.05 0.25)
(0.3 0.05 0.25)
);
}
);

fields
(
rhoFlux
rho
phi
);

where I also tried with and without commenting out the argument line, this does not make a difference though...

[msrinath80]

OK, I think this might help. If rhoFlux is a vector (I am assuming it represents ElectricField), then try using:

rhoFlux.component(0)

instead of just rhoFlux in the fields section. That should return the X-component of rhoFlux. If that works, then

rhoFlux.component(1)
rhoFlux.component(2)

will give you the other two components.

Nevertheless I am unsure why you will need to do this. Also please check whether the file rhoFlux contains something meaningful.

[sverdoold]

mmm I think we're getting close to the problem now..(although I'm not sure whether it is easy to solve...)
When I change the sampleDict to use rhoFlux.component(0) I get the following output:

Exec : sample . BoundaryTest1
Date : Oct 07 2007
Time : 02:06:19
Host : sjaak
PID : 1363
Root : /home/sjaak/OpenFOAM/sjaak-1.4.1/run/electrostaticFoam
Case : BoundaryTest1
Nprocs : 1
Create time

Create mesh for time = 0

Deleting samples/ directory

Time = 0
Sampling volScalarField rho
Sampling volScalarField phi
Writing fields to "/home/sjaak/OpenFOAM/sjaak-1.4.1/run/electrostaticFoam/BoundaryTest1/samples/0/ Potential_rho_phi.xy"

Time = 0.0005


--> FOAM FATAL ERROR : component function not supported for field rhoFlux.component(0) of type surfaceScalarField

From function sample
in file sample.C at line 569.

FOAM exiting

So what I understand from this is that basically rhoFlux is not a vector which is weird. I assumed that it was a vector because looking at the solver, it should really represent the electric field. (I am still using the "standard" solver with only the AUTO_WRITE values in for rhoFlux changed in createFields.H)

Looking at createFields.H rhoFlux was indeed defined as surfaceScalarField which I currently think is weird...(the gradient of a field will be a vector isn't it?). Changing the field type to surfaceVectorField let however to many error when I tried to recompile the solver:

createFields.H:70: error: no matching function for call to 'Foam::GeometricField<foam::vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>::GeometricField(Foam::IOobject, Foam::tmp<foam::geometricfield<double,> >)'
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:611: note: candidates are: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const Foam::GeometricField<type,>&, const Foam::wordList&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:576: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const Foam::GeometricField<type,>&, const Foam::word&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:545: note: Foam::GeometricField<type,>::GeometricField(const Foam::word&, const Foam::tmp<foam::geometricfield<type,> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:512: note: Foam::GeometricField<type,>::GeometricField(const Foam::word&, const Foam::GeometricField<type,>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:480: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const Foam::GeometricField<type,>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:448: note: Foam::GeometricField<type,>::GeometricField(const Foam::tmp<foam::geometricfield<type,> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:416: note: Foam::GeometricField<type,>::GeometricField(const Foam::GeometricField<type,>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:378: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, Foam::Istream&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:337: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:313: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, const Foam::dimensionSet&, const Foam::Field<type>&, const Foam::PtrList<patchfield<type> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:283: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, const Foam::dimensioned<type>&, const Foam::wordList&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:254: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, const Foam::dimensioned<type>&, const Foam::word&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:227: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, const Foam::dimensionSet&, const Foam::wordList&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:197: note: Foam::GeometricField<type,>::GeometricField(const Foam::IOobject&, const typename GeoMesh::Mesh&, const Foam::dimensionSet&, const Foam::word&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
electrostaticFoam.C:59: error: no match for 'operator=' in 'rhoFlux = Foam::operator*(const Foam::tmp<foam::geometricfield<double,> >&, const Foam::tmp<foam::geometricfield<double,> >&) [with PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh](((const Foam::tmp<foam::geometricfield<double,> >&)((const Foam::tmp<foam::geometricfield<double,> >*)(& Foam::fvc::snGrad(const Foam::GeometricField<type,>&) [with Type = double]()))))'
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:1029 : note: candidates are: void Foam::GeometricField<type,>::operator=(const Foam::GeometricField<type,>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:1054 : note: void Foam::GeometricField<type,>::operator=(const Foam::tmp<foam::geometricfield<type,> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
/home/sjaak/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude/GeometricField.C:1090 : note: void Foam::GeometricField<type,>::operator=(const Foam::dimensioned<type>&) [with Type = Foam::Vector<double>, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]
electrostaticFoam.C:65: error: no matching function for call to 'div(Foam::surfaceVectorField&, Foam::volScalarField&)'
make: *** [Make/linuxGccDPOpt/electrostaticFoam.o] Fout 1

At this moment those errors don't tell me much. I probably should take a closer look at many header files. However, I'll first get some sleep to have a clear mind when I do that (it's already 2:30 am over here...).
Please do not hesitate to give me some hints or tips if you see some obvious things....

[gschaider]

No. This doesn't help. rhoPhi is a surfaceScalarField (a scalar defined on the cellFaces). So OF is correct about not having any components.

I'm not 100% sure whether sample supports surfaceFields (but that the fact that it tries to do something with it says yes). I think that your choice of an interpolation-scheme ("cell") might be the problem: if the values are defined on the cell-faces it's kinda hard to use cell-values for interpolation

[hjasak]

Well, a surface field is not a good representation of a vector field - it exists only on faces. You can do something like

volVectorField rhoVelocity = fvc::reconstruct(rhoPhi):

and then sample that. Tis will produce a cell-centred equivalent of the information stored in the face fluxes.

Hrv

[mike_jaworski]

Hi, I'm a little late to the discussion, but maybe I can still add something...
I'm new to OF as well and figured electrostaticFoam would be a good place to start learning. it's also something I'm more familiar with (solving EM problems vs. cfd) so I'm a little more familiar with the physics.

Unfortunately, I'm not familiar with c++ as much, so I'm undergoing some confusion as well.

I don't think rhoFlux is the electric field but the divergence of the current density, but maybe one of the more versed coders can correct me. Here's why I think this:

The definition of rhoFlux is as follows:

rhoFlux = -k*mesh.magSf()*fvc::snGrad(phi)

If you look in the programmer's manual under divergence, it gives an integral form transforming the volume integral of the divergence into a [closed] surface integral.

Discretizing this you end up taking the summation of cell face areas multiplied by the surface normal value (see eq. 2.24 Programmers Manual) of the vector.

This summation of the cells is, I believe, accomplished by:
mesh.magSf()

k (I assume) is the conductivity so -k*grad(phi) = -k*E = J (the current density).
Or, div J = div( -k E) = div(-k grad(phi)) which, when you convert to a discretized object is the summation over the cell surfaces of the surface-normal gradient of the potential or

rhoFlux = -k*mesh.magSf()*fvc::snGrad(phi)

so rhoFlux is the divergence of current density, I think.

Now, hopefully someone can correct me if I'm wrong. I'm learning c++ as I go, but I think I pulled the physical meaning of the code back out.

In order to get an electric field, then, one can define the appropriate class/variable (just like rhoFlux) in createFields.H and also in electroStaticFoam.C for the following:

volVectorfield E = fvc::grad(phi)

This situation is a little odd, though, since the equation for the time derivative of rho involves the convective derivative (fvm::div(rhoFlux,rho)) and I don't recall that being present in Maxwell's equations, but I haven't a resource at the moment to double check this. More later I suppose.

Hopefully someone can double check what I think the code is saying since I'm a novice and I haven't had a chance to try it out.

Regards,
Mike

[mike_jaworski]

Hi again,
I forgot a negative sign in the calculation for electric field strength:

E = -1*fvc::grad(phi)

Also, I check with Maxwell's equations and rederived charge continuity.

d(rho)/dt - div J = 0 is what it should be.

If I understand the coding, this equation *should* be represented as:

fvm::ddt(rho) - rhoFlux == 0

Or if you define E as above

fvm::ddt(rho) + fvc::div(k * E) == 0

As long as k is a scalar. This would be an interesting way to do it (and I'm going to try it tomorrow) but if one defines k as a vector to represent an anisotropic medium, this would be:

fvm::ddt(rho) + fvc::div(k & E) == 0

One could add an anisotropic dielectric permitivity in the mix as well by altering the first equation in the solver:

fvm::laplacian(epsilonRelative, phi) + rho/epsilon0 == 0

where epsilonRelative is the dielectric tensor in relative units.

Since the "== 0" don't appear anywhere, I'm guessing the solver assumes the equation in the "solve(... );" term sums to zero.

I'll give a lot of these a shot tomorrow, but I do think there might be an error in the way the code is set up now, but I may not be c++ literate enough to fully understand it.

Regards,
Mike

P.S. Sorry for the double post.

[Daniele111]

hi
I would use electrostaticFoam. My domain contains insulated object and conductor, what kind of boundary conditions i must use? On conductor I impose the potential phi but on on insulated part?
Thanks
Daniele