How to compute energy density from electric field
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
# compute the electric field
electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg')
electric_field *= -1.0
out['electric_field'] = Struct(name='Electric field', mode='cell', data=electric_field, dofs=None)
# compute the energy density
energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef, ef)', ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
Thank you,
Geoff
On 08/07/2014 10:43 PM, Geoff Wright wrote:
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
Yes, this is easy.
# compute the electric field electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg') electric_field *= -1.0 out['electric_field'] = Struct(name='Electric field', mode='cell',data=electric_field, dofs=None)
# compute the energy density energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef, ef)',ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
You cannot use values in QP or element averages (energy_density) in place of term arguments - those has to be field variables. It is possible to get QP values back into a variable, but I would not recommend it - there would be interpolation involved and precision loss. Instead, compute all the derived quantities in QP:
from sfepy.linalg import dot_sequences
from sfepy.discrete.common.mappings import get_jacobian
# Electric field E in QP.
ef = pb.evaluate('ev_grad.i1.Omega(v)', mode='qp')
# E^T E in QP.
ed = dot_sequences(ef, ef, 'ATB')
# Scalar sigma in QP.
mat = pb.get_materials()['p']
sigma = mat.get_data(('Omega', 'i1'), 0, 'val')
# Energy density in QP.
ed = sigma * ed
# Jacobian in QP.
var = pb.get_variables()['v']
det = get_jacobian(var.field, pb.integrals['i1'])
# Element averages.
energy_density = (ed * det).sum(1) / det.sum(1)
energy_density = energy_density[..., None]
out['energy_density'] = Struct(name='Energy density', mode='cell',
data=energy_density, dofs=None)
# Total energy.
total_energy = (ed * det).sum()I tried the above in the sphere_disc_inside_sphere.py and it worked. Let me know it if is what you need.
r.
Yes - it works. Thanks!
On Friday, August 8, 2014 4:38:10 AM UTC-4, Robert Cimrman wrote:
On 08/07/2014 10:43 PM, Geoff Wright wrote:
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
Yes, this is easy.
# compute the electric field electric_field = pb.evaluate('ev_grad.i1.Omega( v )',mode='el_avg')
electric_field *= -1.0 out['electric_field'] = Struct(name='Electric field', mode='cell',data=electric_field, dofs=None)
# compute the energy density energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef,ef)',
ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
You cannot use values in QP or element averages (energy_density) in place of term arguments - those has to be field variables. It is possible to get QP values back into a variable, but I would not recommend it - there would be interpolation involved and precision loss. Instead, compute all the derived quantities in QP:
from sfepy.linalg import dot_sequences from sfepy.discrete.common.mappings import get_jacobian # Electric field E in QP. ef = pb.evaluate('ev_grad.i1.Omega(v)', mode='qp') # E^T E in QP. ed = dot_sequences(ef, ef, 'ATB') # Scalar sigma in QP. mat = pb.get_materials()['p'] sigma = mat.get_data(('Omega', 'i1'), 0, 'val') # Energy density in QP. ed = sigma * ed # Jacobian in QP. var = pb.get_variables()['v'] det = get_jacobian(var.field, pb.integrals['i1']) # Element averages. energy_density = (ed * det).sum(1) / det.sum(1) energy_density = energy_density[..., None] out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None) # Total energy. total_energy = (ed * det).sum()I tried the above in the sphere_disc_inside_sphere.py and it worked. Let me know it if is what you need.
r.
Hi Robert,
Quick question about the above -- when I visualize the electric field as returned from pb.evaluate('ev_grad.i.Omega(v)', mode='el_avg') it looks a little ugly because its drawing cell values i.e. each triangle has a constant value. For a nicer visualization I need to evaluate the gradient per vertex and save this into vtk in vertex mode. Then paraview will interpolate across the faces. So my question is which function to call to get the per vertex gradient? I tried calling ev_grad with mode='el' and it seems to not be implemented for ev_grad term because the shape I get back is the same as the el_alg.
Thanks,
Geoff
On Friday, August 8, 2014 10:25:43 AM UTC-4, Geoff Wright wrote:
Yes - it works. Thanks!
On Friday, August 8, 2014 4:38:10 AM UTC-4, Robert Cimrman wrote:
On 08/07/2014 10:43 PM, Geoff Wright wrote:
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
Yes, this is easy.
# compute the electric field electric_field = pb.evaluate('ev_grad.i1.Omega( v )',mode='el_avg')
electric_field *= -1.0 out['electric_field'] = Struct(name='Electric field', mode='cell',data=electric_field, dofs=None)
# compute the energy density energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef,ef)',
ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
You cannot use values in QP or element averages (energy_density) in place of term arguments - those has to be field variables. It is possible to get QP values back into a variable, but I would not recommend it - there would be interpolation involved and precision loss. Instead, compute all the derived quantities in QP:
from sfepy.linalg import dot_sequences from sfepy.discrete.common.mappings import get_jacobian # Electric field E in QP. ef = pb.evaluate('ev_grad.i1.Omega(v)', mode='qp') # E^T E in QP. ed = dot_sequences(ef, ef, 'ATB') # Scalar sigma in QP. mat = pb.get_materials()['p'] sigma = mat.get_data(('Omega', 'i1'), 0, 'val') # Energy density in QP. ed = sigma * ed # Jacobian in QP. var = pb.get_variables()['v'] det = get_jacobian(var.field, pb.integrals['i1']) # Element averages. energy_density = (ed * det).sum(1) / det.sum(1) energy_density = energy_density[..., None] out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None) # Total energy. total_energy = (ed * det).sum()I tried the above in the sphere_disc_inside_sphere.py and it worked. Let me know it if is what you need.
r.
On 08/14/2014 04:26 PM, Geoff Wright wrote:
Hi Robert,
Quick question about the above -- when I visualize the electric field as returned from pb.evaluate('ev_grad.i.Omega(v)', mode='el_avg') it looks a little ugly because its drawing cell values i.e. each triangle has a constant value. For a nicer visualization I need to evaluate the gradient per vertex and save this into vtk in vertex mode. Then paraview will interpolate across the faces. So my question is which function to call to get the per vertex gradient? I tried calling ev_grad with mode='el' and it seems to not be implemented for ev_grad term because the shape I get back is the same as the el_alg.
This works:
from sfepy.discrete import FieldVariable
from sfepy.discrete.fem import Field
ef = pb.evaluate('ev_grad.i1.Omega( v )', mode='qp')
gfield = Field.from_args('grad', nm.float64, 3,
pb.domain.regions['Omega'])
var = FieldVariable('ef', 'parameter', gfield,
primary_var_name='(set-to-None)')
var.set_data_from_qp(ef, pb.integrals['i1'])
out.update(var.create_output())but beware: the gradient is not defined in the vertices - the FE basis is only C^0 there. So the above approach averages the gradient from the elements sharing the vertex, see [1]. This may smooth out some interesting details.
r.
[1] http://sfepy.org/doc-devel/src/sfepy/discrete/fem/fields_base.html?highlight...
Thanks,
Geoff
On Friday, August 8, 2014 10:25:43 AM UTC-4, Geoff Wright wrote:
Yes - it works. Thanks!
On Friday, August 8, 2014 4:38:10 AM UTC-4, Robert Cimrman wrote:
On 08/07/2014 10:43 PM, Geoff Wright wrote:
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
Yes, this is easy.
# compute the electric field electric_field = pb.evaluate('ev_grad.i1.Omega( v )',mode='el_avg')
electric_field *= -1.0 out['electric_field'] = Struct(name='Electric field', mode='cell',data=electric_field, dofs=None)
# compute the energy density energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef,ef)',
ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
You cannot use values in QP or element averages (energy_density) in place of term arguments - those has to be field variables. It is possible to get QP values back into a variable, but I would not recommend it - there would be interpolation involved and precision loss. Instead, compute all the derived quantities in QP:
from sfepy.linalg import dot_sequences from sfepy.discrete.common.mappings import get_jacobian # Electric field E in QP. ef = pb.evaluate('ev_grad.i1.Omega(v)', mode='qp') # E^T E in QP. ed = dot_sequences(ef, ef, 'ATB') # Scalar sigma in QP. mat = pb.get_materials()['p'] sigma = mat.get_data(('Omega', 'i1'), 0, 'val') # Energy density in QP. ed = sigma * ed # Jacobian in QP. var = pb.get_variables()['v'] det = get_jacobian(var.field, pb.integrals['i1']) # Element averages. energy_density = (ed * det).sum(1) / det.sum(1) energy_density = energy_density[..., None] out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None) # Total energy. total_energy = (ed * det).sum()I tried the above in the sphere_disc_inside_sphere.py and it worked. Let me know it if is what you need.
r.
Okay yeah that works -- thanks!
Geoff
On Thursday, August 7, 2014 4:43:08 PM UTC-4, Geoff Wright wrote:
Hi Robert,
Me again :) but I think this one is easy... in the previous thread about flux I was computing electric field from voltage in the post processing step. I would like to compute energy density also, from electric field. I believe its something like this, in the post_process hook:
# compute the electric field electric_field = pb.evaluate('ev_grad.i1.Omega( v )', mode='el_avg') electric_field *= -1.0 out['electric_field'] = Struct(name='Electric field', mode='cell',data=electric_field, dofs=None)
# compute the energy density energy_density = pb.evaluate('dw_volume_dot.i1.Omega(p.sigma, ef,ef)', ef=electric_field) out['energy_density'] = Struct(name='Energy density', mode='cell', data=energy_density, dofs=None)
However, I get an error when its parsing the variables. How do I specify ef variable properly? I attached the relevant files in case its necessary.
Thank you,
Geoff
participants (2)
-
Geoff Wright -
Robert Cimrman