Calculating flux from simulation results - Improving Interpolation

[Jexyll_S]

Hi everyone,
I am trying to calculate the magnetic flux through a small, circular aperture using python scripts. My current process is to extract the values for the B-field from the simulation results, then add all the B-fields in a radius and multiply by the aperture area. The code for this is below.

import numpy
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit
	# Creating the analysis pipeline
	# Adding a new SimulationExtractor
	simulation = document.AllSimulations["Self Inductance TMS Device"]
	simulation_extractor = simulation.Results()

	# Adding a new EmSensorExtractor
	em_sensor_extractor = simulation_extractor["Overall Field"]
	em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All"
	document.AllAlgorithms.Add(em_sensor_extractor)

	B_field = em_sensor_extractor.Outputs["B(x,y,z,f0)"]
	B_field.Update()
	out = B_field.Data.Field(0)
	
	# Extract and correct dimensionality
	dims = np.array(B_field.Data.Grid.Dimensions)
	if np.prod(dims) != len(out):
		dims = dims - 1
	
	# get x component of B-field. Reshape to correspond to volume
	# Creates an array of [x,y,z] values
	out_Bx = out[:,0].reshape(dims, order = 'F')
	
	# Get cooordinates for cell nodes (not cell centres)
	xaxis_nodes = field.Data.Grid.XAxis
	yaxis_nodes = field.Data.Grid.YAxis
	zaxis_nodes = field.Data.Grid.ZAxis
	
	# Extrapolate coordinate points
	xaxis = (xaxis_nodes[:-1] + xaxis_nodes[1:])/2
	yaxis = (yaxis_nodes[:-1] + yaxis_nodes[1:])/2
	zaxis = (zaxis_nodes[:-1] + zaxis_nodes[1:])/2
	
	# Define circular apperture
	centre = np.array([0, 0, -15])*1e-3
	radius = 1*1e-3
	
	# Find index closest to x-coordinate of centre
	centre_idx_x = (np.abs(xaxis - centre[0])).argmin()
	
	# Calculate flux through apperture
	flux = 0.0
	for y_idx, ycoord in enumerate(yaxis):
		for z_idx, zcoord in enumerate(zaxis):
			if ((ycoord - centre[1])**2 + (zcoord - centre[2])**2) <= radius**2:
				# print("y and z coords are {} and {}".format(ycoord, zcoord))
				# print("Value being added to flux is %f" % str(out_Bx[centre_idx_x, ycoord, zcoord]))
				flux = flux + out_Bx[centre_idx_x, y_idx, z_idx]
	
	flux = flux * np.pi * radius**2
	print("Final flux is {}".format(flux))

In addition, I also attempted using the flux evaluator tool. However, when I create circle entities, I am unable to extract any surfaces. Instead, I had to use an infinitely thin cylinder and apply suitable transforms to position it appropriately.

import numpy as np
import s4l_v1.analysis as analysis
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.units as units
from s4l_v1 import ReleaseVersion
from s4l_v1 import Unit
from s4l_v1 import Vec3

	# Define the version to use for default values
	ReleaseVersion.set_active(ReleaseVersion.version5_2)
	
	cylinder2 = model.CreateSolidCylinder(Vec3(0,0,0),Vec3(0,0,0),0.5)
	rotation = model.Rotation(Vec3(0,1,0), np.pi/2)
	translation = model.Translation(Vec3(0,0,-15))
	
	cylinder2.ApplyTransform(rotation)
	cylinder2.ApplyTransform(translation)
	
	cylinder2.Name = "Cylinder 2"
	
	# Creating the analysis pipeline
	# Adding a new SimulationExtractor
	simulation = document.AllSimulations["Self Inductance TMS Device"]
	simulation_extractor = simulation.Results()

	# Adding a new ModelToGridFilter
	inputs = []
	model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs)
	model_to_grid_filter.Name = "Cylinder 2"
	model_to_grid_filter.Entity = model.AllEntities()["Cylinder 2"]
	model_to_grid_filter.MaximumEdgeLength = 7.07106769085e-05, units.Meters
	model_to_grid_filter.UpdateAttributes()
	document.AllAlgorithms.Add(model_to_grid_filter)

	# Adding a new EmSensorExtractor
	em_sensor_extractor = simulation_extractor["Overall Field"]
	em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All"
	document.AllAlgorithms.Add(em_sensor_extractor)

	# Adding a new FieldInterpolationFilter
	inputs = [em_sensor_extractor.Outputs["B(x,y,z,f0)"], model_to_grid_filter.Outputs["Surface"]]
	field_interpolation_filter = analysis.core.FieldInterpolationFilter(inputs=inputs)
	field_interpolation_filter.UpdateAttributes()
	document.AllAlgorithms.Add(field_interpolation_filter)

	# Adding a new SurfaceFieldFluxEvaluator
	inputs = [field_interpolation_filter.Outputs["B(x,y,z,f0)"]]
	surface_field_flux_evaluator = analysis.core.SurfaceFieldFluxEvaluator(inputs=inputs)
	surface_field_flux_evaluator.UpdateAttributes()
	surface_field_flux_evaluator.Update()
	document.AllAlgorithms.Add(surface_field_flux_evaluator)
	
	
	flux_data = surface_field_flux_evaluator.Outputs["Total Flux(B(x,y,z,f0))"].Data
	mag_flux = flux_data.GetComponent(0)[0]
	print("magnetic flux is {}".format(mag_flux))

I prefer the second method, since it interpolates values, providing a more accurate measure of the flux. I wanted to post this code to provide some more references which new users of sim4life can draw from, and also to get feedback on my coding techniques.

Cheers,
Jexyll

[Sylvain]

This looks quite good, thank you for providing your scripts.
Note that, for the second method, you could create your cylinder directly at the desired position:

cylinder2 = model.CreateSolidCylinder(Vec3(0,1,-15),Vec3(0,1,-15),0.5)

There is also the XCoredModeling.CoverWires() function that allows you to make a surface (i.e. face) out of a circle entity.

Last, but not least, note that the line model_to_grid_filter.MaximumEdgeLength is ultimately what determines the tradeoff between accuracy and computational cost of the interpolation (since it defines the resolution of the triangulated mesh on which the interpolation is done).