Calculate diffusive flux from solution in interactive mode
I've made a lot of progress learning sfepy this weekend but I'm stuck on post-processing. I'm in interactive mode (using Jupyter) using a toy model that solves the laplace equation on a rectangle with dirichlet conditions on the ends (image below). Now, I'd like to evaluate the flux at designated positions, say on a line or a set of points.
I'm able to evaluate the flux field using:
flux=pb.evaluate('ev_diffusion_velocity.3.Omega(s.D, c)', copy_materials=False, mode='qp')
where 'c' is the solution field variable and s is the material. This produces a set of flux vectors at quadrature points. What I can't seem to do is put this data in a proper field variable so that I can evaluate it with probes, say for example at the boundaries. I've looked at the examples but have found them too complex to reduce to a few lines of code. Can anyone suggest the easiest way to get this data into a field variable?
https://lh3.googleusercontent.com/-IGio8BtAJJE/WOF2SY3wmCI/AAAAAAAAAKI/q9S1I...
Hi Scott,
check the example [1] - it does pretty much all that you described, i.e., putting a quadrature point-based quantity into a FE field, and then probing that field. Let us know if you need more details, or help with reducing the code to a your problem.
r.
[1] http://sfepy.org/doc-devel/examples/linear_elasticity/its2D_interactive.html
On 04/03/2017 03:12 AM, Scott Calabrese Barton wrote:
I've made a lot of progress learning sfepy this weekend but I'm stuck on post-processing. I'm in interactive mode (using Jupyter) using a toy model that solves the laplace equation on a rectangle with dirichlet conditions on the ends (image below). Now, I'd like to evaluate the flux at designated positions, say on a line or a set of points.
I'm able to evaluate the flux field using:
flux=pb.evaluate('ev_diffusion_velocity.3.Omega(s.D, c)', copy_materials=False, mode='qp')
where 'c' is the solution field variable and s is the material. This produces a set of flux vectors at quadrature points. What I can't seem to do is put this data in a proper field variable so that I can evaluate it with probes, say for example at the boundaries. I've looked at the examples but have found them too complex to reduce to a few lines of code. Can anyone suggest the easiest way to get this data into a field variable?
https://lh3.googleusercontent.com/-IGio8BtAJJE/WOF2SY3wmCI/AAAAAAAAAKI/q9S1I...
Robert, thank you for your response. I've now written a function that calculates the gradient as a field variable from a scalar field variable:
def gradient(pb,region,scalar,optorder=1): from sfepy_scb.projections import project_by_component import numpy as nm
gfield = Field.from_args('gfield', nm.float64, 2, region,
approx_order=optorder - 1)
gradient = FieldVariable('gradient', 'parameter', dfield,
primary_var_name='(set-to-None)')
cfield = Field.from_args('component', nm.float64, 1, region,
approx_order=optorder - 1)
component = FieldVariable('component', 'parameter', cfield,
primary_var_name='(set-to-None)')
ev = pb.evaluate
order = 2 * (optorder - 1)
theterm= 'ev_grad.%d.%s(%s)' % (order,region.name,scalar.name)
gradient_qp=ev(theterm, mode='qp')
gradient_qp=nm.swapaxes(gradient_qp,2,3) #make the row vectors into columns
project_by_component(gradient, gradient_qp, component, order)
return gradientI can then use the output to evaluate line probes:
gradc= gradient(pb,omega,c) p=LineProbe([0,0],[5,0],0) gline=p(gradc)
Note that in the second line of gradient() I've imported my own version of project_by_component. That's because the original version expects a 3-component gradient and mine has only 2 components. I had to add a check for the vector dimension:
dim = nm.shape(tensor_qp)[2]
for ic in range(dim):My next task is to submit that as an issue.
I would also like to make a similar function to calculate flux, using something like:
theterm= 'ev_diffusion_velocity.%d.%s(s.D, %s)' % (order,region.name,scalar.name) flux_qp=ev(theterm, copy_materials=False, mode='qp')
The problem is that I don't know how to pass the material.diffusivity in as an argument. The best is could do is pass a string 's.D' because s.D does not evaluate to anything. Any suggestions?
Hi Scott,
On 9.4.2017 20:56, Scott Calabrese Barton wrote:
Robert, thank you for your response. I've now written a function that calculates the gradient as a field variable from a scalar field variable:
def gradient(pb,region,scalar,optorder=1): from sfepy_scb.projections import project_by_component import numpy as nm
gfield = Field.from_args('gfield', nm.float64, 2, region, approx_order=optorder - 1) gradient = FieldVariable('gradient', 'parameter', dfield, primary_var_name='(set-to-None)') cfield = Field.from_args('component', nm.float64, 1, region, approx_order=optorder - 1) component = FieldVariable('component', 'parameter', cfield, primary_var_name='(set-to-None)') ev = pb.evaluate order = 2 * (optorder - 1) theterm= 'ev_grad.%d.%s(%s)' % (order,region.name,scalar.name) gradient_qp=ev(theterm, mode='qp') gradient_qp=nm.swapaxes(gradient_qp,2,3) #make the row vectors intocolumns
project_by_component(gradient, gradient_qp, component, order) return gradientI can then use the output to evaluate line probes:
gradc= gradient(pb,omega,c) p=LineProbe([0,0],[5,0],0) gline=p(gradc)
Note that in the second line of gradient() I've imported my own version of project_by_component. That's because the original version expects a 3-component gradient and mine has only 2 components. I had to add a check for the vector dimension:
dim = nm.shape(tensor_qp)[2] for ic in range(dim):My next task is to submit that as an issue.
OK, thanks for submitting the issue - it seems like an easy-to-do improvement.
I would also like to make a similar function to calculate flux, using something like:
theterm= 'ev_diffusion_velocity.%d.%s(s.D, %s)' % (order,region.name,scalar.name) flux_qp=ev(theterm, copy_materials=False, mode='qp')
The problem is that I don't know how to pass the material.diffusivity in as an argument. The best is could do is pass a string 's.D' because s.D does not evaluate to anything. Any suggestions?
Either you can pass the material as a keyword argument to Problem.evaluate() (like "s=s", with s = Material(D=...)), or you can create the term instance directly, as in tests/test_high_level.py: test_term_evaluation() - there is unfortunately an example for this. Let us know if you need more help - in that case send also your script so that it can be tried.
r.
participants (2)
-
Robert Cimrman -
Scott Calabrese Barton