Calculating flux from simulation results - Improving Interpolation
-
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 -
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).