Electrostatics with volume charge

[drmike]

I am learning how to use Elmer and it is clear there is a lot I don't understand. As a starting point I'd like to put a charge distribution in a vacuum volume and just compute potential and E field. How do I feed the charge density to the StatElecSolver?

[kevinarden]

https://www.nic.funet.fi/pub/sci/physic ... Manual.pdf

Section 16 says you can apply a charge density as a body force

[drmike]

I see - I want to follow the example "Body Force as a Function of Space" in the solver manual. Thanks!

[drmike]

It's clear I don't know what I'm doing. I created a simple function and then introduced a body force with

Code: Select all

Body 1
  Equation = 1
  Material = 1
  Body Force = 1
End
...
Body Force 1
     Charge Density
     Variable Potential
     Real Proceedure "initrho" "initrho"
End
Solver 1 :: Reference Norm = 5.84522573E-01

With the variable being a dummy input. But I get an error:

Code: Select all

ERROR:: LoadInputFile:  Unknown specifier: [reference norm]
ERROR:: LoadInputFile:  In section: [body force 1]
ERROR:: LoadInputFile:  For property name:[solver 1]
STOP 1

Clearly the body force section is totally wrong. How do I make it correct?

[kevinarden]

Body Force 1
Charge Density = Variable Potential
Real Procedure "initrho" "initrho"
End

[drmike]

I changed it to

Code: Select all

Body Force 1
     Charge Density Potential
     Real Procedure "initrho" "initrho"
End

and it runs. I think the main problem was that I spelled "Procedure" wrong!!! Hard to see two e's where there is only supposed to be one.
Now I just have to make the charge distribution charged!

[drmike]

I have a lot more to learn...

The output reports

Code: Select all

Loading user function library: [initrho]...[initrho]
LoadInputFile: Number of BCs: 2
LoadInputFile: Number of Body Forces: 1

and I put the line

Code: Select all

     Call INFO('initrho', 'Processing INIITRHO');
     firsttime = .false.

in the "first time" loop. But it never shows up in the output. So the function is loaded, but never actually called.
I have the .sif file attached, hopefully there is something obviously stupid that will get a charge density to show up.

[kevinarden]

Your sif worked fine on my system

StatElecSolver: Starting Assembly...
GetNOFActive: Number of active elements: 25819
initrho: Processing INIITRHO
DefUtils::DefaultDirichletBCs: Setting Dirichlet boundary conditions

perhaps the problem is in initrho can you post that?

[drmike]

Sure:

Code: Select all

Function initrho(model, n, dummyin) result(rho)

  USE DefUtils
  Implicit None

! function variables

  Type(Model_t) :: model
  Integer :: n
  Real(Kind=dp) :: dummyin, rho

! internal variables

  Integer :: dim
  Real(Kind=dp) :: x, y, z, r1, r2, r, rp
  Logical :: firsttime = .true.

! keep these variables

  Save dim, firsttime, r1, r2

  if(firsttime) then
     dim = CoordinateSystemDimension()
     if(dim /= 3) then
        Call FATAL('initrho', 'Must be 3D')
     end if
     Call INFO('initrho', 'Processing INIITRHO');
     firsttime = .false.
     r1 = 1.0
     r2 = 0.35
  end if

! get local coordinates

  x = model % Nodes % x(n)
  y = model % Nodes % y(n)
  z = model % Nodes % z(n)
  rp = sqrt(x**2 + y**2)
  r = sqrt(x**2 + 2*y**2 + z**2 - 2*rp*r1 + r1**2)

  rho = 8.6e-5*exp(-12.5*(r/r2)**2)
end function initrho

I think I found one problem - I was playing with caf and cafrun, and I may have gotten elmerf90 to use caf instead of f95 during my fooling around. This morning all that is gone and I now see

Code: Select all

ElmerSolver par0.sif | grep rho
Loading user function library: [initrho]...[initrho]
initrho: Processing INIITRHO

I will play with the amplitude of the charge and see if that's my problem.
Edit: Changing the amplitude from e-5 to e+5 made no difference on the potential. Removing the factor 12.5 made no difference either. So the charge density is being computed, but not as "charge". So I'm still missing something.
As I said, I have a lot to learn! Thanks for helping.

[kevinarden]

it works here, I add a print line just to make sure

initrho.f90

(907 Bytes) Downloaded 8 times

StatElecSolver: Starting Assembly...
GetNOFActive: Number of active elements: 25819
initrho: Processing INIITRHO
I am here
DefUtils::DefaultDirichletBCs: Setting Dirichlet boundary conditions
SetNodalLoads: Checking for nodal loads for variable: potential