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