Hi Sam,

I know it seems a bit distracting, but it's sometimes nice to write up your elog entries in a more readable form. I've re-edited your comments to make them more easily read.

I used http://www.sciweavers.org/free-online-latex-equation-editor to convert your latex to png files for easy reading.

At this point, my are goals are to:

This complex conversion already seems to have been done in the paper, since the stationary solution of the form:

 

 

is assumed. k_n and ℓ_m are obtained from the adiabatic boundary conditions.

 Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12.

( although I'm not quite sure how to take the second derivative of a zeroth-order bessel function).

 

I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules.  

At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

Here is a set of slides by Yoichi Aso on how to handle gravity in Comsol.

This tutorial will go through how to show frequency convergence for a cylindrical cantilever using Ansys 14.5 and Mathematica.

The aLIGO model available on IFOSim uses the first surface of the mirror as Port 1 and the other surface as Port 2. This requires us to use negative radii of curvatures among others to get the correct modeling. However, this makes adding surface maps more difficult and was causing issues with getting the simulations to follow what is expected. The surface_map attribute for mirrors doesn't allow us to specify the mirror port in Finesse.

Following a discussion with Kevin, it seems that this approach is being deprecated. His suggestion is to use node 1 as the front of a mirror and node 2 as the back of the mirror. And likewise nodes 1 and 2 as the front of a beamsplitter and nodes 3 and 4 as the back. 

I am rewriting the IFOSim aLIGO model with this convention to allow for later simulations not to be plagued by issues and allow for easier simulations of the beam.

Attached are plots of how HOMs are introduced by defects in the mirror. This is just an example for a basic FP Cavity.

Initially, ideal conditions were considered, without any deformations or thermal heating effects in the mirrors. This offered a baseline model from which subsequent tests could draw comparisons.

Following this, more complex scenarios were introduced into the model: mirror-shape defects caused by suspension and thermal lensing due to heat absorption. These heat absorption rates were set at 0.5 ppm in the ETM coating and 0.3 ppm in the ITM coating.

The lensing calculations for these scenarios were performed using the power observed in the non-deformed mirror case. This lensing needs to be made dynamic in future codes to get the modeling correct. Nonetheless, introducing these conditions resulted in the emergence of higher-order modes (HOMs). Unfortunately, there seems to be a loss of power in the cavity that I haven't been able to resolve yet. Notably, this power loss persisted across the range of mirror position scans showing that this is not due to just a change in resonance length.

Also, the Hello Vinet Function in Finesse 3 seems to yield a higher degree of lensing than typically reported in the literature. I need to check with someone about what assumption I am taking incorrectly here. It's plausible that the heightened lensing effect is contributing to the observed power loss in the cavity, possibly due to high curvatures.

Looking forward, my focus will be on determining the source of this lensing discrepancy first. This analysis will involve a more detailed examination of the variables contributing to thermal lensing and potential causes for the higher-than-expected lensing observed in the simulation. Ultimately, the aim is to create a dynamic deformation model that can more accurately represent these phenomena within the aLIGO setup before introducing ring heaters and controls.

I have coded a rudimentary triangular model using parameters provided by Aaron. In this model, I have kept a simple PDH signal which is the REFL power demodulated with the EO frequency. The PDH signal is then fed to the laser to modulate the frequency.

There are some points of confusion I have in this model such as the high loop gain required to achieve locking. I think this may be due to the fact that the reflective power in SI units is incredibly small while the frequency shift is considerably larger in SI units.

Then I took this model through basic operations such as determining the optimal phase angle for the PDH signal and assessing if lock can be achieved. They went as required for a sanity check.

Then I introduced a mismatch for light incident on the first mirror. This resulted in the creation of higher-order modes in the cavity as expected.

I am still analyzing the effect of this mismatch for the length control offsets. Also, I am not sure why we need a sensing matrix in this situation, because we're currently dealing with a single signal source and one degree of freedom. Or maybe a more complex readout and control scheme may be more helpful. This remains to be seen,

Nonetheless, I've attached plots that illustrate the cavity modes that have emerged due to the mismatch.

At this point, my are goals are to 1) convert the time-dependent heat equation into stationary, complex form, 2) use the Levin approach to calculate the TR noise given this stationary PDE, and 3) verify the results in COMSOL

 I have looked at the Heinert paper and converted the time-dependent partial differential Heat equation to a stationary, complex one.  I did this by assuming a temperature distribution of the form T(\vec{r},t) = T_0 ( \vec{r} ) e^{i \omega t }, where  \omega is the frequency of oscillation of the input heat field: q( \vec{r} , t) = q( \vec{r} ) e^{i \omega t}.  Hence, I am assuming the steady-state solution of the temperature field.  From this ansatz, the PDE becomes

C_p T(\vec{r}) + \frac{ i \kappa}{ \omega }  \nabla^2 T ( \vec{r} ) = q ( \vec{ r} ).  

This complex conversion already seems to have been done in the paper, since the stationary solution of the form

T(r, z) = \sum_{n = 0}^{\infty} \sum_{m = 0}^{ \infty} T_{n m} J_0 (k_n r) \cos ( \ell_m z) is assumed.

k_n and \ell_m are obtained from the adiabatic boundary conditions.  Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12. ( although I'm not quite sure how to take the second derivative of a zeroth-order bessel function).

I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules.  At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

 I carried out the calculations for Heinert's simple model in greater detail and was able to obtain equation 12 (with the exception of a plus sign instead of a minus sign in the denominator; I guess the Fourier transform method gives a minus sign.  In the end, that doesn't matter, since the modulus square yields the same result in the denominator.)  I obtained the same spectral density S_z (\omega), with the exception of a missing factor of R^2.  I don't know where this R^2 went, or why it is that the spectral density is associated with the extrinsic variable z (instead of r, say).  It's supposed to be a "homogeneous [noise] readout along the z direction."

 It does appear that the simplified model is only relevant for the simulations.  To quote Heinert: "An efficient computation is only possible for the simple model, as the advanced model would require an element of size more than 10⁶ ."  I have run Koji's code that replicates Heinert's figure 3.  I have attached the resulting temperature distribution and noise amplitude curve.  In the noise amplitude curve, the red line is the analytical result, while the dots are from COMSOL.

The next step is to convert this code to an efficient complex time-independent solution.  As stated before, my main concern here is whether COMSOL actually solves the right equation in the stationary case.

 Here, I describe some test cases to see if COMSOL's solutions are agreeing with some simple analytical solutions.  Right now, I have two plots showing COMSOL's solution and my analytical solution on separate plots.  I will be plotting there difference to see if they really match up.

 The following document shows the relative difference between these two plots.

Agreement with Heinert's paper for cylindrical TR noise has now been achieved.  Using the stationary state assumption to calculate the temperature profile, the computation time was reduced compared to the previous time-dependent approach. Here are the plots showing the agreement.  I have shown the plots for a 1D axisymmetric model, in addition to a full 3D model in COMSOL.  Both give the same result.

(See Plots in attached document)

My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V).  I am trying to match Duan's analytical expression with Heinert's analytical expression.  This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. )   Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".

 I have plotted Heinert's analytical solution for TR noise using Duan's parameters.  Since TO and TE noise can be found by simply rescaling TR noise, these have been included in the plot as well.  The solid curve represents the analytical solution, while the tick marks represent COMSOL's solution.  I have used COMSOL for both a 1D axisymmetric and a 3D model.  Since Duan's cylinder has a radius of 125 microns, but a length of 1 m, the meshing was difficult for the 3D model.  I ended up shortening the length of the cylinder, converting to the actual length when finally calculating the thermal noise.

 I have now plotted the Duan-Heinert comparison for the case of an infinite upper bound.  It turns out that the curves differ by a maximum of 10 percent for low frequencies.  Such a discrepancy has been attributed to lack of experimental investigation into this regime (according to Duan). For our purposes, such a discrepancy is acceptable. We will therefore use the Heinert curve for subsequent calculations due to its faster computation time.

 In this document, I try to identify I good mesh by comparing the numerical solution from that mesh with my analytical model.  Since there are problems with carrying out the analytical calculation, it is still not entirely clear which mesh should be used.

 I have refined the analytical calculation procedure, as outlined in this new document.  The procedure indicates that the discrepancy between the analytical and numerical solutions are more likely attributed to meshing inaccuracies.

Meshing Surface Layers

Defining New Selections

I don't know why I wasn't seeing this problem with previous models (perhaps because I wasn't importing any CAD or STEP files), but my latest attempts at meshing and selecting specific domains of my model were being thwarted by inconsistent domain definitions. I was previously always manually selecting domains, which is confusing because all domains just get assigned a number when they are created. Worse, sometimes the numbers assigned to domains change when the model rebuilds, especially if there has been a significant change in the model geometry. This results in later steps selecting the wrong set of domains (or boundaries, etc).

To fix this problem, I created new selections (sets of domains, boundaries, or etc that receive their own label and can be selected as a group later in the model). The new selections include:

This is probably also a necessary step for getting reliable results when interfacing with MATLAB, and might explain some weird problems I was running in to a while ago and just made haphazard fixes for.

Mesh Steps

I use the following mesh steps to get what seems like a pretty reliable meshing:

1. Mesh the upper tapers (bulk and surface separately) with a free tetrahedral mesh
   1. I mostly use the defaults for an extremely coarse mesh, but the only parameter that seems to make a large difference is the minimum mesh size. I set this to the skin_depth for the surface layer, and (fiber_main_radius-skin_depth) for the bulk. fiber_main_radius-skin_depth should be the radius of the bulk domain, and the skin_depth characterizes the smallest length scale in the surface layers. I have some limited ability to tweak the minimum mesh sizes when the surface layer is comparable in size to the fiber_main_radius (so the surface layer is comparable to the entire fiber radius), but it seems to be best to keep the mesh this fine when the surface layer becomes small.
2. Mesh the main fiber (bulk and surface separately) with a swept mesh
   1. I create a distribution with a fixed number of elements, at ceil(fiber_main_length/fiber_stock_radius/2). This is somewhat arbitrary--the stock has the largest radial length scale, so I figured I'd divide up the main fiber into units that tall. A better thing would be to know how high a mode we are interested in studying, and break up the fiber into enough pieces to observe that mode. Seems fine for now, might want this distribution to be a bit coarser though.
3. Mesh the lower tapers (bulk and surface separately) with a free tetrahedral mesh
   1. Use the same size settings as on the upper tapers
4. Mesh the thick sections (bulk and surface separately) with a swept mesh
   1. The distribution uses a fixed number of ceil(fiber_thick_length/fiber_stock_radius/2) elements. Again perhaps this can be coarser; it also doesn't attempt to make a finer mesh at the thermoelastic cancellation region.
5. Mesh the necks (bulk and surface separately) with a free tetrahedral mesh.
   1. I use an extremely coarse mesh with the minimum mesh element size set to fiber_thick_radius-skin_depth for the bulk; for the surface it is an extra coarse mesh instead, because the extremely coarse mesh gave low quality mesh elements. I'm not sure why this is, I can't see much difference and the problem only seems to happen to one of the 8 neck sections (why not in all sections?). The problem only arises when the skin depth is less than 20um, for the other parameters at the nominal aLIGO values.
6. Mesh the stocks (bulk and surface separately) with a swept mesh
   1. Distribution is again ceil(fiber_stock_length/fiber_stock_radius/2), resulting in 2 vertical divisions of the stock.
7. Mesh the horns with a free tetrahedral mesh
   1. Size is an extremely coarse mesh, where the minimum element size is set to the skin_depth (because the horn sees the boundary of the stock surface, so its smallest adjacent element can have a scale down to the skin_depth), and maximum growth rate is increased to 9 (any higher results in low quality mesh elements; I set it high because most of the horn does not need to be meshed as finely as the part directly in contact with the fiber stock).

Note on the skin depth: I defined the thickness of the surface layer by the parameter skin_depth. Since we are interested in how the energy in the surface layer changes with the radius of the fiber, and expect that the surface depth doesn't change much with fiber radius (a very thick fiber will have a few micron surface layer; so will a fiber half that diameter). I managed to get the mesh working with a 15um surface layer, but was having trouble getting a 10um layer. I wasn't completely sure how small a surface layer would be necessary, but for a test mass 1/100 the size of the aLIGO masses, the main part of the fiber would have a radius of 20um, so if we want to study the surface layer for masses down to that size I figure the skin depth should be at most around 10-15um. It would be better to be motivated by the actual scale of the physics going on at the surface--how deep to micro cracks and other lossy imperfections go?

Study

I'm running the frequency study on a small number (2) of frequencies in each decade from 60Hz to 10kHz (so there should be 5-8 total frequencies in the study). I started it at 4:25pm on my local machine, and it hasn't gotten very far in 30 minutes, so I may abort the study and try to make the mesh coarser--especially the distribution settings, which are easy to change.

I had some trouble running this on optimus.

Optimus has COMSOL 5.1 installed, but I made these files in 5.3. I downloaded the comsol 5.3 dvd.iso file last night, but on install I'm now getting the error "No locks available." I wasn't sure if this is a file permissions issue (sometimes the file has been 'locked' when it is saved and I open it in the GUI), an svn issue (most of the information I find online is from people running into some problem with lockd or their NFS), or a licensing issue (I don't think this is likely, COMSOL on my laptop shows more licenses available).

Has anyone seen this? Do I need to do something else with /cvs after install?

Gautam advised me against trying to install version 5.3, lest it break version 5.1--I had already gone through the install, but looking at the install manual it says it shouldn't affect previous installs except that the default behavior when double clicking .mph files will always choose the latest version of COMSOL. Since we mostly run on the command line we should be fine. That said I haven't tested files with COMSOL 5.1.

Zach pointed me to sandbox1.ligo.caltech.edu at the lunch meeting, as well as these notes from Rana about running comsol on a remote server. I couldn't run the file myself on sandbox1 because I don't yet have a home directory there, but I asked Larry if he could set one up for me so I can use the clusters. Gautam helped me run them on sandbox1 by having me scp the file onto the shared 40m drive, then he moved it in to his home directory on sandbox1 and ran it. A version of the model asking for an analysis at 60Hz and 600Hz ran in 10-20 min, and looked good (I was able to scp the output to my laptop and open it in the gui, though I have a script not quite to the point where Matlab can do this for me, it's not ready yet and anyway I wanted to directly see that the output was normal). I modified the study to take 100 steps between 60Hz and 5kHz, and Gautam has now started that run, which if it scales linearly (hopefully better than linear, since it won't have to remesh and etc) will take about 8-16 hours.

Thank you Gautam for the extended mattermost session helping me run these!

Simulation results

First run

Relevant properties of the file that might need to be tweaked:

• F_load
  • The applied force may have been resulting in large displacements of the test mass (and thus large curvature in the fibers). Since the model contains geometric nonlinearities (I think it has to, no?), large displacements may not be able to be linearly extrapolated back to the small displacements we expect. Since I am running a frequency domain study, COMSOL expects a load force, which means I am specifying a force rather than a displacement. In future runs, I wanted to reduce the displacement by scaling down the force... I'll discuss my first attempt at this later.
• skin_depth
  • Gabriele mentioned that when he models surface losses for wafers he does not explicitly mesh a boundary layer, but rather asks COMSOL for a boundary integral (so just the 2D integral on the outer surface of the wafer), and then presumably makes some assumption about the thickness of this layer to go from J/m to J if need be. In contrast, my suspension thermal noise model defines the 'surface' of the fiber as the 15um outer layer, and I explicitly mesh that surface volume and do a 3D (volume) integral over the surface domain.
  • Using a 2D integral would significantly simplify the model and probably run faster and result in fewer bugs when the scale of the geometry changes the necessary fineness of the mesh. However, I am a bit uncertain about whether this is what we want. At least in my understanding, we are interested in asking about what happens when the surface layer depth becomes comparable to the radius of the fiber. Since the 2D integral would always represent an infinitesimal volume, and the geometry of an extended (deep) surface layer would differ significantly from that of an infinitesimal shell, so I am concerned that the strain energy density would vary significantly with the radial postition in the fiber. Nonetheless, even the mesh I am using now has only a few (2-4) mesh elements along the radial direction, so I'm not sure I'm doing much better even with the additional computation time. Therefore, in future runs I'd like to try Gabriele's suggestion.
• Frequencies
  • I was using 50Hz spacing on the frequencies from 60Hz to 5kHz. I haven't been very systematic about this, but I'm getting some convergence issues going to lower frequencies--maybe by taking the suggestion for eliminating the skin_depth I can make a finer mesh and go to lower frequencies. This frequency sweep is too coarse to find well-defined resonances, but they are suggested in the plot below.

More recent runs

I wanted to automate this loop over mass parameters in a matlab script, so I set up the Livelink handshake so Matlab would send the model to sandbox1 for solving, then MATLAB on my machine would work with the solved model to extract results. I realized later during running that this might not be optimal, since it will require me to keep the connection to the remote COMSOL server during the entire run, which is A. annoying and B. risky because these might take many hours to run and the connection can easily be severed, wrecking the job. I don't really have a workaround for this, so my current plan is to continue logging on to the remote server and running individual COMSOL files. As I'll discuss, this is probably necessary for now anyway since I ran in to some problems running these models, and might want to tweak the models on my local GUI before running them over a large (frequency) parameter space on sandbox1. 

aLIGO parameters, but decrease F_load

Based on the plots generated from the first run, I estimated that the TM displacement was comparable to the TM thickness; just as a first pass I figured 'small displacement' could be easily defined as 'displacements smaller than the thickness of the fiber' or maybe 'smaller than the thickness of the surface layer'. If the displacement is directly proportional to the force, this meant I wanted to scale the previous F_load by 1e-5, which I did. I wanted to just see what kinds of displacements I would get, so I asked to run the model only at 300Hz. After more than 45 min of running, the server threw an error that it couldn't find my STEP file for the horns, which were in a local directory but not in the remote (server's) directory.

I thought this was an odd error for two reasons. First, I had originally not even had the STEP in my local directory (it only needed to load once and then could be used many times for a given model in the GUI), and I was getting an error within a couple minutes of starting the job that it could not find the file. Adding the file to my local directory seemed to solve this problem, as the model was running for much longer. The error I got after it finally crashed was not the same error as before, but was still an error in loading the file, which makes me a bit confused about where it is actually looking for this file or when exactly the loading process happens. To solve this, I have copied the STEP to the remote directory and will run again with the import pointing to that remote file, which I suspect will solve the problem.

aLIGO with the mass scaled by 0.1, also decrease F_load by 1e-5*0.1

Since the first run seemed to mostly work up to the outstanding questions mentioned above, I decided to also run the model in batch mode directly from the sandbox1 command line, just as Gautam had. That is, I ssh on to sandbox1 and run the command

$ /usr/local/comsol53/multiphysics/bin/glnxa64/comsol batch -inputfile [inputname].mph -outputfile [outputname].mph

For this run, I scaled the mass of the model by 0.1, so we now have a 4kg TM. Most lengths of the model then scale by 0.1^(1/3), but those related to the radius of the fibers (and the bend point of the fibers) scale by 0.1^(1/2) because we want to maintain a constaint stress in the fiber. The main length of the fiber remained the same so the modes would be in about the same place. I document these more thoroughly in the matlab script, which I should upload to the git. I scale the force by 1e-5 as above (same back of the envelope calculation), with an additional factor of 0.1 to account for the lower mass.

The model ran in to some convergence issues after a few frequencies. I could solve this by changing the relative tolerance, but I think I will most likely instead pursue Gabriele's suggestion above and try to refine the mesh to improve the convergence. The solver nodes of the COMSOL models are still a bit mysterious to me--I don't really know much about when different methods or measures of convergence are appropriate. Probably playing with these could both improve the accuracy and efficiency of the model... but it seems like a hefty undertaking. Might be worht the long term payoff though.

The first few frequencies did converge, but as I was extracting the resuls ancha, the COMSOL license server, went down so I'll have to extract them later.

NOTE TO SELF: WHAT'S WITH THE INLINE IMAGE QUALITY?                DON't USE JPG !!!!

I made a copy of Gautam's 40m model to add the unstable filter cavity for the ponderomotive squeezing project. I wanted to make a more explanatory record of the changes I've made because I think some of them might be necessary for other scripts using gautam's original model, but I have not implemented them in that file (also just for my own paper trail).

Changes:

• swapped the role of nXBS and nYBS; before I think it was sending the reflected beam at the main BS to the X arm, and the transmitted beam to the Y arm
• I kept nPOY, but I am a bit confused by it--this is the beam returning from the x-arm that is transmitted in to the BS, but it isn't followed out of the BS; rather it seems this pick off beam is detected inside the BS substrate, which is odd. Anyway we aren't using it for now.
• Changed some labels to distinguish between the beamsplitter already in the IFO (IBS) and the BS used in the quantum phase compensator (QBS)
• Added a quantum phase compensator cavity, consisting of QP1 and QP2, as well as a mirror (QP0) to direct light transmitted from the QPC through QBS back to QBS

I talked with Shoaib about some changes I could make to the FEA model to improve convergence and reduce memory usage. Summary:

• use a hex mesh rather than tetrahedral
• Use more structured meshes. In particular, I can make an angled swept the mesh in the tapered portions rather than using a free mesh in these regions, defining the mesh only on one boundary
• Use a nonconformal mesh, so adjacent domains do not need to have matching meshes. This could allow transitions to coarser meshes as the fiber becomes thicker or contacts the horn
• try using curved elements so the curvature isn't just approximated by the number of regular mesh elements
• Might try scaling the model uniformly by some factor (100, 1000), which could avoid machine precision issues

I'll get these implemented this week and see if the computation goes through.

I started implemented some of these changes:

• Started the mesh with a boundary free quad mesh on the interface between the upper tapers and the main part of the fiber. I used the following size setting
  • Maximum element size is fiber_taper_length, which I felt was a good characteristic maximum because it wouldn't make sense to have radial elements larger than the length of the taper itself. This setting does not limit the mesh (elements are much smaller than the maximum)
  • Minimum element size: skin_depth/10, where I use the thickness of the surface layer as the characteristic length, and allow a mesh with elements substantially smaller than this characteristic length.
  • Maximum element growth rate I set at 8, which from some trial and error isn't strongly affecting the mesh as long as it isn't brought much lower than this (so the mesh is allowed to get coarse relatively quickly, though other settings limit this)
  • Curvature factor is 0.7, and actually does affect the mesh strongly; for now value seems to give a 'pretty' (symmetric to the eye) mesh.
  • The resolution of narrow regions is set to 5, where this value again strongly affects the mesh and this value was chosen because it looks 'pretty'
  • Note that in contrast to the previous way I had defined the mesh, I now am defining the surface and bulk of the mesh in the same step, which seems to help these two regions match up. Shoaib had mentioned using a nonconformal mesh, which I may end up implementing in a later update.
• Instead of doing a free mesh in each non-cylindrical region, I continued with a swept mesh along the length of the fiber (still sweeping in sections, though I'm not sure this is necessary)
• I also slightly increased the number of mesh elements along the length of the fibers (ie, used more elements in each section of the mesh sweeping).
• In this model, I still have a free tet mesh on the horn and test mass, since those seemed a bit more involved to move over to a hex mesh (there is no free hex mesh in COMSOL that I could see right away, and I'm not sure that the horn can be meshed with a swept mesh due to its irregularity. I'll look in to this further, but I had no immediate solution for making a hex mesh in the horn) and anyway they are much coarser than the fibers so I don't expect them to cause much trouble.

These changes were fine with a 10um mesh, but ended up with much too fine a mesh (~80000 boundary mesh elements in the initial mesh between the taper and main fiber) due to meshing the surface and boundary in a single stage in the first step. I separate them out to get a coarser mesh in the bulk, and also make the resolution a bit coarser.

When running comsol with matlab interface on sandbox1, it is usually most convenient to ssh with screen forwarding (eg '-CY') and launch COMSOL with matlab by following the instructions in the livelink manual.

Sometimes, it is necessary to run COMSOL without any display available. In that case, the Instructions in the manual are a little unclear. Here are the detailed steps that let me run my script '/home/amarkowi/metamaterials/run_spiral.m' with no screen forwarding.

To run COMSOL on sandbox1 with no graphical Interface, here are the steps that worked for me (Tue Aug 11 16:35:51 2020) from a mac on the Caltech VPN.

We could run the simulations on the 32 core machine in the 40m lab (optimus)? I think Mariia was running some of her studies on optimus, and even though we had some problems with the licensing initially, I think she resolved these and has detailed the procedure in her elogs...

It's a good idea, I'll check out her elogs and get it started tonight.

I found that I had the relative tolerance set too low (0.001, while the iterations were not converging further than 0.02 or so); I changed it to 0.1, which let me run over a few modes relatively quickly once I reduced the number of sections on the main part of the fiber to 10. This is not sufficient--at 600Hz the fiber looks completely jagged because the mesh is not fine enough, so using optimus will be necessary. 

 Here are the plots with their fractional residuals: abs(S_COMSOL - S_analytical)/S_analytical

 When you want to find out who's using up all the licenses, you can run lmstat -a do find out.

Specifically, with Mac COMSOL, I do:

> cd /Applications/COMSOL43a/license

> maci64/lmstat -a -c license.dat

Users of COMSOL:  (Total of 5 licenses issued;  Total of 1 license in use)

  "COMSOL" v4.3, vendor: LMCOMSOL

  floating license

    dmassey c22042.local (v4.3) (ancha/1718 1074), start Thu 12/13 18:26

etc.

 PDF please instead of EPS or BMP or JFIF or TARGA or GIF or ascii art.

I moved the only entry from the 'ENG_FEA' log into the COMSOL log and then renamed that logbook as 'FEA' since we don't need two FEA logs.

Also renamed 'AdhikariLab' log as ATF.

This summarizes how to get the remote comsol server to run. COMSOL 5.1.0.234 is now on tegmeni thanks to Larry.

On the server:

rana@tegmeni|~> /usr/local/comsol51/multiphysics/bin/comsol server -login force

This starts up a comsol server instance, listening on port 2036. '-login force' will ask you to supply a username and password which you make up. You will have to enter these later from the client side.

On my laptop, from the MATLAB prompt:

>> addpath('/Applications/COMSOL51/Multiphysics/mli/')

>> mphstart('tegmeni.ligo.caltech.edu', 2036,'uname','pword')

That's it! Now you just run the matlab script which runs the COMSOL file.

I'm attaching a tarball of the .mph file (written by Dmitry Kopstov from MSU) and the associated matlab scripts which change the parameters, as well as looping over test mass thickness to produce a plot of Brownian noise PSD v. 

I took Aidan's COMSOL model for the ITM from a couple years ago and updated it with some more details:

1. Through radiative cooling only, the ITM is cooled to 103 K. Taking it to 123 K will be accomplished by adding a ring heater to the ITM.
2. Assume 3 W of heating from main laser beam onto ITM HR face.
3. Emissivity of ITM barrel is 0.95. Emissivity of HR* and AR faces is 0.5.
4. The CP and the Inner Shield are kept fixed at 80 K. This is to simulated the effect of having conductive cooling with cold straps. This needs to be checked in more detail by actually modeling thermal straps.
5. Emissivity of the CP is 0.9.
6. The total length of the inner shield is 5 m. The CP is at z = 0 m and the ITM is at z = 2.25 m. We should check what the result would be if the shield is ~1m shorter or longer.

In the attached image, I have made one quadrant of the tubular cryo shield transparent just for clarity - the actual modelled tube is 3 cm thick, made of aluminum, has an emissivity of 0.95 on the inside and 0.03 on the outside (to simulate what we would get from polished aluminum or gold coating).

This files is in our GitLab: https://git.ligo.org/rana-adhikari/CryogenicLIGO/blob/master/FEA/ITM-ColdShield-CP.mph

*I am suspicious of just using a single emissivity number for the AR and HR coatings. Since we are concerned with wavelengths which are long w.r.t. the coating thickness, it may be that the HR and AR coatings have a complicated wavelength dependence in the 5-50 micron band.

I've just tried this out on my desktop machine using COMSOL 5.1 and its still working. Which COMSOL is installed on optimus at the 40m ?

Since 2014, the limit of 12 workers using the matlab parallel computing toolbox has been lifted. Today, I was able to get this to work. There's a trick.

Usually, when you start up matlab and run a parallel thing like 'parfor', it just uses a default profile 'local' which limits you to 12 workers. You can try to ask for more by doing 'parpool(40)' for 40 workers, but it will tell you that NumWorkers = 12 and you're out of luck. So instead:

myCluster = parcluster('local')

myCluster.NumWorkers = 40;

saveProfile(myCluster);

parpool('myCluster', 40)

It seems that it needs the max # of workers and the requested number of workers to be 40 to use 40, otherwise you'll just get 12 (as of matlab 2016a).

When saving your COMSOL files do these two things to make the files much smaller (good for saving in version control and sharing):

1. File -> Compact History
2. Preferences -> Files -> Optimize for File Size (not speed)

Also,

1. click 'Clear Mesh' under the mesh menu
2. 'Clear Solutions' under the Study menu

In this way the file sizes will be ~100 kB instead of 10's of MB.

For the Voyager test masses, we have been considering a barrel coating to increase the IR emissivity to increase the radiative cooling power. We also seek to estimate the added Brownian thermal noise that arises from this.

Dmitry Kopstov (from MSU) made a baseline model for this which we have been modifying. The latest is in the CryogenicLIGO git repo in the FEA directory as ETM_NASAbarrel.mph.

Comparison with Analytic Estimates:

Convergence Tests:

this is a time dependent model of the previous steady-state one

• Cold Shield and CP held at a constant 60 K
• 3 W heat input to the ITM from the main laser beam
• radiative cooling to the shield
• ITM barrel emissivity = 0.9
• ITM HR/AR emissivity = 0.5/0.5

So the cooldown time w/o a heat switch is ~4 days. Since this is less than the usual pumpdown time required to open the gate valves on the beamtubes, perhaps no heatswitch or in-vac cryogens are required.