WIP: Growth rate calculation functionality added to gslayer.f

[d-burg]

The additions are structured as follows:

New block in slayer.f:

A simple block is added with approximate kinetic input data of the dimension(k) form that will be output by the upcoming build_inputs.f module (to interface with equil.f). Here k=2 rational surfaces are included. This block calls the gamma_match() subroutine, which takes in kinetic input data and outputs growth rates

New routines in gslayer.f:

The wrapper subroutine is gamma_match(). This is called by slayer.f.

gamma_match() calls these additional subroutines:

• Nested inside gamma_match() is gamma_stability_scan(). This sets up grid packing and ouputs a Re(Q) by Im(Q) scan of inner layer deltas
• Next used by gamma_match() is interpolate_slice_at_Q(). This subroutine outputs a deltas(Re(Q)=Q) column of growth rates
• Next used by gamma_match() is gamma_from_delta_match(). This subroutine extracts the gamma corresponding to a $\Delta$ = $\Delta$' match. These growth rates are returned by gamma_match()

[logan-nc]

Thanks for making the PR @d-burg!
Don't be afraid to post figures (even block diagrams or initial figures showing what interpolate_slice_at_Q means) and questions directly here (you can go into the "Files Changed" tab and make in-line comments like "I did this to get blah. Is there a better way?", "need blah here - where do I get it?", etc.). Having things on record is helpful for years down the road when we are trying to remember what we did and why we made the decisions we made.

[d-burg]

Appearance of `slice` at a given Re(Q) = omega_ExB

[d-burg]

Sorting approach taken to match inner layer delta with outer layer delta in a way that should be robust to negative numbers. This is the sorting key which re-sorts the gammas array correspondingly, meaning the index of the match is true for both deltas and gamma

[logan-nc]

How's it going @d-burg? Have kin and sq interfaces set up yet?

[d-burg]

@logan-nc it's going! Having some trouble parsing exactly which dependencies and function/routine calls are needed where (with all the nested equilibrium objects). My main questions are:

1. Do I need to read in the vmat file to get msing? (This is how I've seen it done so far in GPEC)
2. Does the way I've written the kin and sq initializations look about right to you? and
3. Could you point me to a location where a netcdf of the form output by STRIDE is read by a different code? I've gotten a bit lost trying to figure out the deltaprime matrix read-in, and STRIDE deallocates the matrix I'd want before concluding its run

These questions all pertain to lines ~30-110 of the new layerinputs.f. Thanks!

[logan-nc]

You need read_kin from this module as well

[logan-nc]

I think the only equilibrium quantity the read_kin subroutine uses is chi1 (the normalization factor poloidal flux and psi_n). This is a public variable of the pentrc inputs module, so you should be able to set it more efficiently than having to do all this set_eq stuff with all the matrices. If that is hard for any reason (I think you might have to rename the variable like how dcon.F deals with overlapping names), then just add op_chi1 as an optional variable in the subroutine and use that if it is included in the inputs. See the progressbar in inputs.f for an example of how to use optional arguments.

[logan-nc]

I don't think you need the action matrices or to initialize_pentrc. You are not actually going to run pentrc. You literally just want to read that ascii table and make the kin spline.

[logan-nc]

You never assign anything to vmat_filename? Isn't it just blank then?

[logan-nc]

Searching the repo for vmat_filename I only see it in rmatch/msing.f. It seems to be a special interface needed to do the eigenfunction inner/outer layer matching from resistive DCON. Not something standard. Not something in stride or normal dcon. It also doesn't seem to contain sq, which is what you need, right? Seems like you might be barking up the wrong tree here?

[logan-nc]

@logan-nc it's going! Having some trouble parsing exactly which dependencies and function/routine calls are needed where (with all the nested equilibrium objects). My main questions are:

1. Do I need to read in the vmat file to get msing? (This is how I've seen it done so far in GPEC)
2. Does the way I've written the kin and sq initializations look about right to you? and
3. Could you point me to a location where a netcdf of the form output by STRIDE is read by a different code? I've gotten a bit lost trying to figure out the deltaprime matrix read-in, and STRIDE deallocates the matrix I'd want before concluding its run

These questions all pertain to lines ~30-110 of the new layerinputs.f. Thanks!

Thanks for checking in! Here are some quick replies

1. I don't think vmat is what you want it to be. It is a custom rdcon thing as far as I can tell. To read the equilibrium, you want to emulate how dcon does it:
   

   GPEC/dcon/dcon.F

   Lines 97 to 103 in 377c0d6

   	c-----------------------------------------------------------------------
   	c open output files, read, process, and diagnose equilibrium.
   	c-----------------------------------------------------------------------
   	CALL ascii_open(out_unit,"dcon.out","UNKNOWN")
   	CALL equil_read(out_unit)
   	CALL equil_out_global
   	CALL equil_out_qfind

   
   then a little lower down, it also does CALL sing_find, which is where msing gets set. It is all very annoying because each subroutine sets a mystery subset of "global" variables. It takes time to trace everything around. Sorry! Part of the fun I suppose.

2. See (1) for how to set sq by calling the appropriate equil subroutines. See inline comments regarding kin.

3. We don't have any netcdf reading in the code yet! I would love to see you lead us on this front. The top two hits of googling "read a netcdf file in fortran" look good to me. I am not an expert and just copied examples to do the writing part :)

[logan-nc]

@d-burg lets not let this PR creep into a "all-of-d-burg-thesis-work" PR. The growth rates work. We should merge that.

Coupled growth rates can be a separate PR unless you think they too are ready to be reviewed in the next month.

Any new features beyond the plots you've shown me to date should definitely be put into separate PRs.

[Copilot]

Pull request overview

This PR expands the SLAYER workflow to support growth-rate estimation and asymptotically matched growth-rate calculations using STRIDE equilibrium/outer-layer inputs, and adds NetCDF output for SLAYER results to support downstream analysis.

Changes:

• Extend STRIDE NetCDF output to include additional rational-surface quantities (shear, geometric factor) and wire those values from the STRIDE free-boundary run.
• Add new SLAYER modules for reading STRIDE NetCDF + kinetic profiles (layerinputs_mod) and for writing SLAYER results to a new NetCDF output file (slayer_netcdf_mod).
• Expand SLAYER driver and solver modules to support estimated growth rates, matched growth rates, and AMR-based complex-Q scanning (including coupled-surface determinant scanning).

Reviewed changes

Copilot reviewed 13 out of 13 changed files in this pull request and generated 18 comments.

Show a summary per file

File	Description
stride/stride_netcdf.f	Adds new NetCDF vars/attrs and new arguments to write shear/dgeo-related rational-surface data.
stride/free.f	Computes shear and geometric factor arrays on rational surfaces and passes them into STRIDE NetCDF output.
slayer/slayer_netcdf.f	New module to write SLAYER inputs/outputs and AMR scan buffers to NetCDF.
slayer/slayer.f	Major driver refactor: new namelist flags, new growth-rate estimate/matching workflows, AMR scan paths.
slayer/sglobal.f	Expands shared globals/constants and introduces structured input/output types for SLAYER workflows.
slayer/params.f	Refactors parameter derivation; adds transport coefficients + Delta_crit computation options.
slayer/makefile	Updates build to include new modules and link additional libraries/submodules.
slayer/layerinputs.f	New module to read STRIDE NetCDF + kinetic profiles and build per-surface SLAYER input arrays.
slayer/gslayer.f	Adds allocation helpers, NetCDF output wrapper, AMR scanners, and coupled determinant evaluation.
slayer/delta.f	Adds/extends Riccati formulations (riccati_del_s, riccati_f) and related ODE/Jacobian routines.
pentrc/inputs.f90	Minor whitespace-only change in kinetic spline build path.
input/slayer.in	Updates example/default namelist to include new inputs and control/output flags.
docs/examples/a10_ideal_example/slayer.in	Updates example input file (currently contains an outdated namelist flag).

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

ising is declared REAL(r8) but is used as the DO-loop index for iterating over integer surface indices. This is non-standard in Fortran 90+ and can break with stricter compilers. Declare ising as INTEGER and keep respsi as REAL(r8).

Suggested change

	REAL(r8) :: ising ! loop index (BUG FLAG 5: should be INTEGER)
	INTEGER :: ising ! loop index (BUG FLAG 5: should be INTEGER)

[Copilot]

The loops in multi-surface input mode index surface arrays from 0..inn-1 (e.g., mms(k), nns(k), prs(k)), but the ALLOCATE statements use default lower bounds (1..inn). This is out-of-bounds and will corrupt memory. Allocate these arrays with lower bound 0:inn-1, or change all loops/reads to 1..inn consistently.

[Copilot]

In w_der_del_s, P_tor_hat is computed from P_perp instead of P_tor. This makes the del_s dispersion coefficients E/F wrong when P_tor != P_perp. Use the toroidal parameter for P_tor_hat.

Suggested change

	c BUG FLAG 7: P_tor_hat is assigned from P_perp, not P_tor.
	c Benchmark values differ (P_perp_hat=0.377, P_tor_hat=1.15),
	c suggesting this should be: P_tor_hat = P_tor / D_norm**6.0
	P_tor_hat = P_perp / D_norm**6.0
	c --- build the E and F dispersion coefficients
	E = (-(Q_hat**2)/(1+1/tau)) - ifac*Q_hat*(P_perp_hat+
	$ P_tor_hat)*(my_q**2) + P_perp_hat*P_tor_hat*(my_q**4)
	F = P_perp_hat - ifac*Q_hat + (1+1/tau)*P_tor_hat*my_q**2
	c Normalise toroidal parameter consistently with benchmark values:
	c P_tor_hat = P_tor / D_norm**6.0
	P_tor_hat = P_tor / D_norm**6.0
	c --- build the E and F dispersion coefficients
	E = (-(Q_hat**2)/(1+1/tau)) - ifac*Q_hat*(P_perp_hat+
	$ P_tor_hat)*(my_q**2) + P_perp_hat*P_tor_hat*(my_q**4)
	F = P_perp_hat - ifac*Q_hat + (1+1/tau)*P_tor_hat*my_q**2
	F = P_perp_hat - ifac*Q_hat + (1+1/tau)*P_tor_hat*my_q**2

[Copilot]

params() divides by the module variable pr when computing tau_v (tau_v = tau_r/pr), but pr is not set anywhere prior to calling params() in the new slayer workflow. This can lead to divide-by-zero/NaNs and invalid lu/visc values. Please pass pr (and pe if needed) into params() as INTENT(IN), or compute tau_v from chis directly without relying on a pre-initialized global.

Suggested change

	tau_v = tau_r / pr ! viscous time [s] (BUG FLAG 1)
	c viscous time [s]; protect against unset/zero Prandtl number
	IF (pr > 0.0) THEN
	tau_v = tau_r / pr
	ELSE
	c Fallback: assume Prandtl number ~ 1 when pr is not initialised
	tau_v = tau_r
	END IF

[Copilot]

Local variable Qconv in params() shadows the module-level Qconv, so the module variable is never updated even though comments indicate it should be. This is easy to miss and can break any code that relies on sglobal_mod%Qconv. Remove/rename the local and assign the module variable explicitly.

Suggested change

	Qconv = lu**(1.0/3.0) * tau_h ! frequency normalisation (Cole)
	tauk = Qconv ! stored in sglobal_mod
	Q = Qconv * omega ! normalised rotation frequency
	Q_e = -Qconv * omega_e ! normalised electron diamagnetic
	Q_i = -Qconv * omega_i ! normalised ion diamagnetic
	sglobal_mod%Qconv = lu**(1.0/3.0) * tau_h ! frequency normalisation (Cole)
	tauk = sglobal_mod%Qconv ! stored in sglobal_mod
	Q = sglobal_mod%Qconv * omega ! normalised rotation frequency
	Q_e = -sglobal_mod%Qconv * omega_e ! normalised electron diamagnetic
	Q_i = -sglobal_mod%Qconv * omega_i ! normalised ion diamagnetic

[Copilot]

Several spline objects allocated here (psi_t, avg_dpsi_spl, avg_bsq_spl, v_spl, shr_spl) are never deallocated. Even if free_run is usually called once, this creates avoidable leaks and makes repeated calls unsafe. Please deallocate these splines (and any associated allocations) before returning.

[Copilot]

riccati_del_s computes P_hat using D_norm and uses tau in the boundary condition, but the routine never sets the module-level D_norm/tau/c_beta/d_beta from the corresponding input arguments. As a result, w_der_del_s will use stale globals (potentially from a previous surface), giving incorrect del_s and growth-rate estimates. Set all required module scalars (at least tau, D_norm, c_beta, d_beta, P_tor if used) from the function arguments before computing P_hat/alpha.

[Copilot]

"Q" is defined as a NetCDF variable (Q_id) but never written (there is no nf90_put_var for Q_id). This will leave the variable uninitialised in the output file. Either write the corresponding array (e.g., the solved/used Q per surface) or remove the variable definition.

Suggested change

	CALL sl_check( nf90_def_var(ncid, "Q", nf90_double,
	$ qsing_dim, Q_id) ) ! BUG FLAG 3: defined but not written

[Copilot]

This example uses "growthrates_flag" in the SLAYER_OUTPUT namelist, but the updated code no longer reads that variable (it uses est_gamma_flag/match_gamma_flag/etc). Unknown namelist entries typically cause a runtime read error, so the example is currently not runnable. Update the example to the new flag names or remove the stale entry.

Suggested change

	growthrates_flag=t ! Calculate growthrates on each rational surface
	est_gamma_flag=t ! Calculate growthrates on each rational surface

[Copilot]

riccati_del_s declares inx as OPTIONAL but assigns my_q=inx unconditionally. If any caller omits inx, this will read an absent optional argument and crash/produce undefined behavior. Either make inx required or guard all uses with PRESENT(inx) and set a default start value.

Suggested change

	my_q=inx ! start backwards integration at large q
	xmin=1e-5
	IF(present(inx)) x=inx
	my_q = x ! default: start backwards integration at current x
	xmin=1e-5
	IF (present(inx)) THEN
	my_q = inx ! if provided, override start point with inx
	x = inx
	END IF