diff --git a/sirepo/package_data/static/html/silas-thermal-transport.html b/sirepo/package_data/static/html/silas-thermal-transport.html
index 28778ae096..ff5e7fdd8d 100644
--- a/sirepo/package_data/static/html/silas-thermal-transport.html
+++ b/sirepo/package_data/static/html/silas-thermal-transport.html
@@ -18,8 +18,5 @@
-
diff --git a/sirepo/package_data/static/js/silas.js b/sirepo/package_data/static/js/silas.js
index 0e0ef693a5..913db13ea1 100644
--- a/sirepo/package_data/static/js/silas.js
+++ b/sirepo/package_data/static/js/silas.js
@@ -506,6 +506,7 @@ SIREPO.viewLogic('thermalTransportCrystalView', function(appState, panelState, s
$scope.$on('crystal.changed', () => {
appState.models.thermalTransportCrystal.crystal = appState.models.crystal;
+ appState.saveQuietly('thermalTransportCrystal');
});
});
diff --git a/sirepo/package_data/static/json/silas-schema.json b/sirepo/package_data/static/json/silas-schema.json
index 85ca9954ae..d93d410c8d 100644
--- a/sirepo/package_data/static/json/silas-schema.json
+++ b/sirepo/package_data/static/json/silas-schema.json
@@ -220,8 +220,7 @@
"grid_points_r": ["Grid Points R", "Integer", 200],
"grid_points_w": ["Grid Points W", "Integer", 0],
"grid_points_z": ["Grid Points Z", "Integer", 200],
- "edge_fraction": ["Fraction of Radial Edge", "Float", 0.98],
- "boundary_tolerance": ["Boundary Tolerance [cm]", "Float", 0.1]
+ "edge_fraction": ["Fraction of Radial Edge", "Float", 0.98]
},
"watch": {
"_super": ["_", "model", "beamElement"],
@@ -295,7 +294,7 @@
]
},
"tempHeatMapAnimation": {
- "title": "Temperature Heatmap",
+ "title": "Radial/Axial Temperature Profile",
"advanced": []
},
"crystal3dAnimation": {
@@ -438,15 +437,10 @@
["R", [
"grid_points_r"
]],
- ["W", [
- "grid_points_w"
- ]],
["Z", [
"grid_points_z"
]]
- ],
- "edge_fraction",
- "boundary_tolerance"
+ ]
],
"advanced": []
},
diff --git a/sirepo/package_data/template/silas/crystal.py.jinja b/sirepo/package_data/template/silas/crystal.py.jinja
index f9f5de5e4a..bfd99ce870 100644
--- a/sirepo/package_data/template/silas/crystal.py.jinja
+++ b/sirepo/package_data/template/silas/crystal.py.jinja
@@ -1,92 +1,23 @@
-import dolfin
-import fenics
-import mshr
-import numpy
-import pandas
from numpy import array, meshgrid, unique
+from pykern.pkcollections import PKDict
from rslaser.optics import Crystal
from rslaser.thermal import ThermoOptic
-from pykern.pkcollections import PKDict
from sirepo.template import template_common
+import dolfin
+import numpy
-# ---
-
-## set up the problem and define (some of) the main parameters
-
-# # simulation parameters
-# # -- crystal properties
-# # cm^2/s diffusion constant of sapphire
-# a_saph = {{ crystalCylinder_diffusionConstant }}
-# # -- time step
-# T = {{ crystalSettings_time }}
-# n_steps = {{ crystalSettings_steps }}
-# dt = T / n_steps # size of time step
-# nip = {{ crystalSettings_plotInterval }} # number of intervals between records
-# # -- crystal dimensions
-# diam = {{ crystalCylinder_diameter }}
-# leng = {{ crystalLength }}
-# # -- mesh density within cylinder
-# md = {{ crystalSettings_meshDensity }}
crystal = Crystal(params=PKDict({{crystalParams}}))
-# # derived parameters
-# radius = diam / 2 # radius
-# lh = leng / 2 # half-length
-# radius2 = radius * radius
-
-# print("duration: ", T)
-# print("time-step: ", dt)
-# print("plot-step: ", dt * nip)
-
-
-# ---
-
-# domain = mshr.Cylinder(fenics.Point(0, 0, lh), fenics.Point(0., 0., -lh), radius, radius)
-# mesh = mshr.generate_mesh(domain, md)
-# # and define function space
-# V = fenics.FunctionSpace(mesh, 'P', 1)
-
-# ---
-
-# define boundary and initial conditions, and sources
-# T0 = {{ crystalCylinder_T0 }}
-# dT = {{ crystalCylinder_dT }}
-# wdT = {{ crystalCylinder_wdT }}
-# sg_px = {{ crystalCylinder_supergaussian }}
-
-# define Dirichlet boundary condition for sides
-# tol = 1e-13
-# def boundary(x, on_boundary):
-# return on_boundary and fenics.near(x[0]*x[0] + x[1]*x[1], radius2, tol)
-
-# bc = fenics.DirichletBC(V, fenics.Constant(0.), boundary)
-
-# # exp(decay) x SG
-# hse = fenics.Expression(
-# 'T0 + dT * (exp(-(x[2]+l/2)/dl) + exp((x[2]-l/2)/dl)) * exp(-0.5 * pow((x[0]*x[0] + x[1]*x[1])/(w*w), px))',
-# degree=1, T0=T0, dT=dT, w=wdT, l=leng, dl=0.8*leng, px=sg_px)
-
-# # define initial value
-# u_n = fenics.interpolate(hse, V)
-
-# #---
-
-# # define the variational problem: u' = D.∆u
-# u = fenics.TrialFunction(V)
-# v = fenics.TestFunction(V)
-# f = fenics.Constant(0)
-
-# F = u*v*fenics.dx + dt*a_saph*fenics.dot(fenics.grad(u), fenics.grad(v))*fenics.dx - (u_n + dt*f)*v*fenics.dx # w/ Dirichlet + initial condition
-# a, L = fenics.lhs(F), fenics.rhs(F)
-
-
def thermo_optic_sim():
# Initialize a simulator with a crystal object set
- thermo = ThermoOptic(crystal)
- # thermo.mesh # Uncomment to view 3D crystal mesh
+ thermo = ThermoOptic(crystal, {{thermalCrystal.mesh_density}})
# Set evaluation points & boundary conditions
- thermo.set_points((200, 0, 201))
+ thermo.set_points((
+ {{thermalTransportSettings_grid_points_r}},
+ {{thermalTransportSettings_grid_points_w}},
+ {{thermalTransportSettings_grid_points_z}},
+ ), {{thermalTransportSettings_edge_fraction}})
thermo.set_boundary(2*(crystal.radius*100.)**2./40.)
thermo.set_load("{{ pump_pulse_profile }}")
@@ -106,6 +37,7 @@ def thermo_optic_sim():
RZ, ZR = meshgrid(rrs, zs)
temp_profile = numpy.array([rs, Ts[ptzs==zface], zs, Ts[ptrs==rcenter]])
pTs = Ts.reshape((len(rs), len(zs)), order='F')
+ pTs = array((pTs[::-1]).tolist()+pTs[1:].tolist())
heat_map = pTs.T
return PKDict(
thermo=thermo,
@@ -119,8 +51,6 @@ def thermo_optic_sim():
res = thermo_optic_sim()
-# ---
-
# for 3D plots, get the facets, and build an array
# containing the indices of their coordinates
inds = []
@@ -129,45 +59,8 @@ for item in dolfin.cpp.mesh.facets(res.thermo.mesh):
# we will provide these indices to plotly so it can draw proper surfaces
inds = numpy.array(inds)
-#ii = inds[:, 0]
-#jj = inds[:, 1]
-#kk = inds[:, 2]
-
-# get node coördinate values and ranges
-# xvals = thermo.mesh.coordinates()[:, 0]
-# yvals = thermo.mesh.coordinates()[:, 1]
-# zvals = thermo.mesh.coordinates()[:, 2]
-# xmin, xmax = xvals.min(), xvals.max()
-# ymin, ymax = yvals.min(), yvals.max()
-# zmin, zmax = zvals.min(), zvals.max()
-
-# ---
-
-
-#def write_row(name, time, values, mode="a"):
-# with open("{{ crystalCSV }}", mode) as f:
-# f.write("{},{},{}\n".format(name, time, ",".join([str(x) for x in values])))
-
-
-# tol = {{crystalSettings_tolerance}}
-# TODO(pjm): add tolerance to thermalTransportSettings
-tol = 8e-3
-# avoid points outside domain (for a coarse mesh, increase this value)
-#xv = numpy.linspace(xmin * (1 - tol), xmax * (1 - tol), 201)
-#zv = numpy.linspace(zmin * (1 - tol), zmax * (1 - tol), 201)
-#radpts = [(x_, 0, 0) for x_ in xv]
-#axipts = [(0, 0, z_) for z_ in zv]
-
-#write_row("xv", 0, xv, "w")
-#write_row("zv", 0, zv)
-
-# for t, u in evolve():
-# ux = numpy.array([u(pt) for pt in radpts])
-
-# uvals = u.compute_vertex_values()
numpy.save("indices.npy", inds)
numpy.save("vertices.npy", res.thermo.mesh.coordinates())
-# numpy.save("intensity.npy", uvals)
numpy.save("tempProfile.npy", res.temp_profile)
template_common.write_dict_to_h5(
PKDict(
@@ -178,8 +71,8 @@ template_common.write_dict_to_h5(
crystal.params.pop_inversion_mesh_extent,
],
y=[
- -crystal.params.pop_inversion_mesh_extent,
- crystal.params.pop_inversion_mesh_extent,
+ -crystal.params.length / 2,
+ crystal.params.length / 2,
],
),
),
diff --git a/sirepo/template/silas.py b/sirepo/template/silas.py
index 84983e39cb..c9f8326425 100644
--- a/sirepo/template/silas.py
+++ b/sirepo/template/silas.py
@@ -115,24 +115,34 @@ def sim_frame_tempProfileAnimation(frame_args):
f = frame_args.run_dir.join(_TEMP_PROFILE_FILE)
d = numpy.load(str(f))
return template_common.parameter_plot(
- [n for n in .98*d[1 if frame_args.tempProfilePlot == "radialPlot" else 2]],
+ [
+ n * 1e-2
+ for n in 0.98 * d[0 if frame_args.tempProfilePlot == "radialPlot" else 2]
+ ],
[
PKDict(
- points=[n for n in .98*d[0 if frame_args.tempProfilePlot == "radialPlot" else 3]],
- label="(T-T_0), K",
+ points=[
+ n
+ for n in 0.98
+ * d[1 if frame_args.tempProfilePlot == "radialPlot" else 3]
+ ],
+ label="(T-T₀), K",
),
],
PKDict(),
PKDict(
- x_label=("Radial" if frame_args.tempProfilePlot == "radialPlot" else "Longitudinal") + " Position [cm]",
+ x_label=(
+ "Radial"
+ if frame_args.tempProfilePlot == "radialPlot"
+ else "Longitudinal"
+ )
+ + " Position [m]",
),
)
def sim_frame_tempHeatMapAnimation(frame_args):
- with h5py.File(
- frame_args.run_dir.join(_TEMP_HEATMAP_FILE), "r"
- ) as f:
+ with h5py.File(frame_args.run_dir.join(_TEMP_HEATMAP_FILE), "r") as f:
d = template_common.h5_to_dict(f)
r = d.ranges
z = d.intensity
@@ -140,9 +150,9 @@ def sim_frame_tempHeatMapAnimation(frame_args):
title="",
x_range=[r.x[0], r.x[1], len(z)],
y_range=[r.y[0], r.y[1], len(z[0])],
- x_label="Horizontal Position [m]",
- y_label="Vertical Position [m]",
- z_label="Temperature",
+ x_label="Radial Position [m]",
+ y_label="Longitudinal Position [m]",
+ z_label="Temperature (T-T₀), K",
z_matrix=z,
)
@@ -447,6 +457,7 @@ def _generate_parameters_file(data):
v.pump_pulse_profile = _get_crystal(data).pump_pulse_profile
v.crystalLength = _get_crystal(data).length
v.crystalCSV = _CRYSTAL_CSV_FILE
+ v.thermalCrystal = _get_crystal(data)
return res + template_common.render_jinja(SIM_TYPE, v, "crystal.py")
if data.report in _SIM_DATA.SOURCE_REPORTS:
data.models.beamline = []
@@ -462,13 +473,7 @@ def _generate_parameters_file(data):
def _get_crystal(data):
- crystals = [
- x for x in data.models.beamline if x.type == "crystal" and x.origin == "new"
- ]
- for e in crystals:
- if e.id == data.models.thermalTransportCrystal.crystal_id:
- return e
- return crystals[0]
+ return data.models.thermalTransportCrystal.crystal
def _initial_intensity_percent_complete(run_dir, res, data, model_names):