elm-pb-mri Version of elm-pb with diffusion terms added

Start of multi-rate model development, with slow rebuild of profiles.

examples/CMakeLists.txt

@@ -18,6 +18,7 @@ add_subdirectory(dalf3)
 add_subdirectory(eigen-box)
 add_subdirectory(elm-pb)
 add_subdirectory(elm-pb-outerloop)
+add_subdirectory(elm-pb-mri)
 add_subdirectory(fci-wave)
 add_subdirectory(finite-volume/diffusion)
 add_subdirectory(finite-volume/fluid)

examples/elm-pb-mri/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 3.13)
+
+project(elm_pb_mri LANGUAGES CXX)
+
+if (NOT TARGET bout++::bout++)
+  find_package(bout++ REQUIRED)
+endif()
+
+bout_add_example(elm_pb_mri
+  SOURCES elm_pb_mri.cxx
+  EXTRA_FILES data/BOUT.inp)

examples/elm-pb-mri/data/BOUT.inp

@@ -0,0 +1,257 @@
+# settings file for BOUT++
+# High-Beta reduced MHD case
+
+##################################################
+# Global settings used by the core code
+
+nout = 40          # number of time-steps
+timestep = 1       # time between outputs
+wall_limit = 1.55  # wall time limit (in hours)
+
+zperiod = 15        # Fraction of a torus to simulate
+MZ = 16             # Number of points in Z
+
+grid = "cbm18_dens8.grid_nx68ny64.nc"  # Grid file
+
+[mesh]
+
+staggergrids = false    # Use staggered grids
+
+[mesh:paralleltransform]
+type = shifted # Use shifted metric method
+
+##################################################
+# derivative methods
+
+[mesh:ddx]
+first = C4  # order of first x derivatives
+second = C4 # order of second x derivatives
+upwind = W3 # order of upwinding method W3 = Weno3
+
+[mesh:ddy]
+first = C2
+second = C2
+upwind = C2
+flux = C2
+
+[mesh:ddz]
+first = C4  # Z derivatives can be done using FFT
+second = C4
+upwind = W3
+
+##################################################
+# FFTs
+
+[fft]
+
+fft_measurement_flag = measure  # If using FFTW, perform tests to determine fastest method
+
+##################################################
+# Solver settings
+
+[solver]
+
+# mudq, mldq, mukeep, mlkeep preconditioner options
+atol = 1.0e-8 # absolute tolerance
+rtol = 1.0e-5  # relative tolerance
+
+use_precon = false    # Use preconditioner: User-supplied or BBD
+
+mxstep = 5000   # Number of internal steps between outputs
+
+##################################################
+# settings for high-beta reduced MHD
+
+[highbeta]
+
+density = 1.0e19       # number density of deuterium [m^-3]
+                       # used to produce output normalisations
+
+evolve_jpar = false     # If true, evolve J raher than Psi
+
+evolve_pressure = true # If false, switch off all pressure evolution
+
+phi_constraint = false # Solve phi as a constraint (DAE system, needs IDA)
+
+## Effects to include/exclude
+
+include_jpar0 = true     # determines whether to include jpar0 terms
+include_curvature = true # include curvature drive term?
+
+compress = false       # set compressible (evolve Vpar)
+nonlinear = false     # include non-linear terms?
+
+diamag = true         # Include diamagnetic effects?
+diamag_grad_t = false  # Include Grad_par(Te) term in Psi equation
+diamag_phi0 = true    # Balance ExB against Vd for stationary equilibrium
+
+##################################################
+# BRACKET_METHOD flags:
+# 0:BRACKET_STD; derivative methods will be determined
+#   by the choices C or W in this input file
+# 1:BRACKET_SIMPLE; 2:BRACKET_ARAKAWA; 3:BRACKET_CTU.
+
+bm_exb_flag = 0
+bm_mag_flag = 0
+##################################################################
+withflow = false     # With flow or not
+D_0 = 130000        # differential potential
+D_s = 20            # shear parameter
+K_H_term = false    # Contain K-H term
+sign = -1           # flow direction
+x0 = 0.855          # peak location
+D_min = 3000        # constant
+##################################################################
+
+eHall = false         # Include electron pressure effects in Ohm's law?
+AA = 2.0          # ion mass in units of proton mass
+
+noshear = false        # zero all shear
+
+relax_j_vac = false    # Relax to zero-current in the vacuum
+relax_j_tconst = 1e-2  # Time constant for vacuum relaxation
+
+## Toroidal filtering
+filter_z = true   # remove all except one mode
+filter_z_mode = 1  # Specify which harmonic to keep (1 = fundamental)
+low_pass_z = 16    # Keep up to and including this harmonic (-1 = keep all)
+zonal_flow = false    # keep this zonal harmonic?
+zonal_field = false     # keep this zonal harmonic?
+zonal_bkgd = true    # keep this zonal harmonic?
+
+## Jpar smoothing
+smooth_j_x = true   # Filter Jpar in the X direction
+
+## Magnetic perturbations
+include_rmp = false # Read RMP data from grid file
+
+simple_rmp = false  # Enable/disable a simple model of RMP
+rmp_n = 3           # Toroidal mode number
+rmp_m = 6           # Poloidal mode number
+rmp_factor = 1.e-4  # Amplitude of Apar [Tm]
+rmp_ramp = 1.e-4    # Timescale [s] of ramp
+rmp_polwid = -1.0   # Width of Gaussian factor (< 0 = No Gaussian)
+rmp_polpeak = 0.5   # Y location of maximum (fraction)
+
+## Vacuum region control
+
+vacuum_pressure = 0.02 # the pressure below which it is considered vacuum
+                       # fraction of peak pressure
+vacuum_trans = 0.01   # transition width (fraction of P)
+
+## Resistivity and Hyper-resistivity
+
+vac_lund = 1.0e8    # Lundquist number in vacuum  (negative -> infinity)
+core_lund = 1.0e8  # Lundquist number in core (negative -> infinity)
+hyperresist = 1.e-4 # Hyper-resistivity coefficient (like 1 / Lundquist number)
+
+## Inner boundary damping
+
+damp_width = -1       # Width of damping region (grid cells)
+damp_t_const = 1e-2  # Damping time constant
+
+## Perpendicular pressure diffusion
+terms_Gradperp_diffcoefs = false
+diffusion_perp = -1.0   # Perpendicular pressure diffusion (< 0 = none)
+
+## Parallel pressure diffusion
+
+diffusion_par = -1.0   # Parallel pressure diffusion (< 0 = none)
+diffusion_p4 = -1e-05   # parallel hyper-viscous diffusion for pressure (< 0 = none)
+diffusion_u4 = -1e-05    # parallel hyper-viscous diffusion for vorticity (< 0 = none)
+diffusion_a4 = -1e-05   # parallel hyper-viscous diffusion for vector potential (< 0 = none)
+
+## heat source in pressure in watts
+
+heating_P = -1   #   heat power in watts (< 0 = none)
+hp_width = 0.1     #   heat width, in percentage of nx (< 0 = none)
+hp_length = 0.3    #   heat length in percentage of nx (< 0 = none)
+
+## sink rate in pressure
+
+sink_P = -1   #   sink rate in pressure (< 0 = none)
+sp_width = 0.04     #   sink width, in percentage of nx (< 0 = none)
+sp_length = 0.15    #   sink length in percentage of nx (< 0 = none)
+
+
+## left edge sink rate in vorticity
+sink_Ul = -10.0        #   left edge sink rate in vorticity (< 0 = none)
+su_widthl = 0.06     #   left edge sink width, in percentage of nx (< 0 = none)
+su_lengthl = 0.1     #   left edge sink length in percentage of nx (< 0 = none)
+
+## right edge sink rate in vorticity
+sink_Ur = -10.0        #   right edge sink rate in vorticity (< 0 = none)
+su_widthr = 0.06     #   right edge sink width, in percentage of nx (< 0 = none)
+su_lengthr = 0.1     #   right edge sink length in percentage of nx (< 0 = none)
+
+## Viscosity and Hyper-viscosity
+
+viscos_par = -0.1   # Parallel viscosity (< 0 = none)
+viscos_perp = -1.0  # Perpendicular viscosity (< 0 = none)
+hyperviscos = -1.0  # Radial hyper viscosity
+
+## Compressional terms (only when compress = true)
+phi_curv = true    # Include curvature*Grad(phi) in P equation
+# gamma = 1.6666
+
+[phiSolver]
+inner_boundary_flags = 4 # INVERT_AC_LAP
+outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP
+
+[aparSolver]
+inner_boundary_flags = 4 # INVERT_AC_LAP
+outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP
+
+##################################################
+# settings for individual variables
+# The section "All" defines default settings for all variables
+# These can be overridden for individual variables in
+# a section of that name.
+
+[all]
+scale = 0.0 # default size of initial perturbations
+
+# boundary conditions
+# -------------------
+# dirichlet    - Zero value
+# neumann      - Zero gradient
+# zerolaplace  - Laplacian = 0, decaying solution
+# constlaplace - Laplacian = const, decaying solution
+#
+# relax( )   - Make boundary condition relaxing
+
+bndry_all = dirichlet_o2 # Default to zero-value
+
+[U]   # vorticity
+scale = 1e-05
+function = ballooning(gauss(x-0.5,0.1)*gauss(y-pi,0.6*pi)*sin(z),3)
+
+[P]  # pressure
+bndry_core = neumann
+#scale = 1.0e-5
+
+[Psi]  # Vector potential
+
+# zero laplacian
+bndry_xin = zerolaplace
+bndry_xout = zerolaplace
+
+[Psi_loc] # for staggering
+
+bndry_xin = zerolaplace
+bndry_xout = zerolaplace
+
+bndry_yup = free_o3
+bndry_ydown = free_o3
+
+[J]    # parallel current
+
+# Zero gradient in the core
+bndry_core = neumann
+
+[phi]
+
+bndry_xin = none
+bndry_xout = none
+#bndry_target = neumann
+