diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5d280a88..5e5b28af 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -74,6 +74,10 @@ jobs: cd "$GITHUB_WORKSPACE" sh ./tests/custom_reference/run_test.sh + # coupled bc test + cd "$GITHUB_WORKSPACE" + sh ./tests/coupled_bcs/run_test.sh + git diff > changes.diff git diff --exit-code --name-only diff --git a/CHANGELOG.md b/CHANGELOG.md index 700b079a..b61ae319 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] ### Added +- Allow the boundary conditions of scalar fields (both the thermal temperature/entropy field and active scalar fields) to be linearly coupled together. This is documented in the "Under Development" section of the User Guide and has an example test. \[Cian Wilson; 6-21-2024; [#460](https://github.com/geodynamics/Rayleigh/pull/460)\] + - The documentation for solving for active and passive scalar fields has been expanded in the "Under Development" section of the User Guide. \[Cian Wilson; 6-20-2024; [#541](https://github.com/geodynamics/Rayleigh/pull/541)\] - A docker image for Stampede3 and Frontera (TACC) has been added to the repository. \[Rene Gassmoeller; 6-20-2024; [#402](https://github.com/geodynamics/Rayleigh/pull/402)\] diff --git a/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt b/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt index 5dd45a39..cee8c4e7 100644 --- a/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt +++ b/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt @@ -1,7 +1,7 @@ **fix_tvar_top** - Logical flag indicating whether thermal variable (T,S) should be fixed on the upper boundary. Default = .true. + Logical flag indicating whether thermal variable (T,S) should be fixed on the upper boundary. Default = .false. **fix_tvar_bottom** - Logical flag indicating whether thermal variable (T,S) should be fixed on the lower boundary. Default = .true. + Logical flag indicating whether thermal variable (T,S) should be fixed on the lower boundary. Default = .false. **fix_dtdr_top** Logical flag indicating whether the radial derivative of thermal variable (T,S) should be fixed on the upper boundary. Default = .false. **fix_dtdr_bottom** @@ -9,7 +9,7 @@ **T_top** Value of thermal variable (T,S) at the upper boundary. Default = 0. **T_bottom** - Value of thermal variable (T,S) at the lower boundary. Default = 1. + Value of thermal variable (T,S) at the lower boundary. Default = 0. **dTdr_top** Value of radial derivative of thermal variable (T,S) at the upper boundary. Default = 0. **dTdr_bottom** diff --git a/doc/source/User_Guide/coupled_bcs.txt b/doc/source/User_Guide/coupled_bcs.txt new file mode 100644 index 00000000..ef9cb6ac --- /dev/null +++ b/doc/source/User_Guide/coupled_bcs.txt @@ -0,0 +1,181 @@ +Rayleigh can couple the values and gradients of scalar fields (including the temperature/entropy, :math:`\Theta`, and active, +:math:`\chi_{a_i}`, scalar fields) at the +boundaries through this initial implementation of coupled boundary conditions. + +These boundary conditions take the form: + +.. math:: + :label: coupled_bcs_T + + \Theta = b_{\Theta,\Theta} + b_{\Theta,\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j} b_{\Theta,\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\Theta,\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling temperature/entropy to its gradient and/or to other scalar fields, + +.. math:: + :label: coupled_bcs_dTdr + + \frac{\partial \Theta}{\partial r} = b_{\frac{\partial \Theta}{\partial r},\frac{\partial \Theta}{\partial r}} + b_{\frac{\partial \Theta}{\partial r},\Theta}\Theta \ + + \sum_{j} b_{\frac{\partial \Theta}{\partial r},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\frac{\partial \Theta}{\partial r},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling the derivative of temperature/entropy to its value and/or to other scalar fields, + +.. math:: + :label: coupled_bcs_chia + + \chi_{a_i} = b_{\chi_{a_i},\chi_{a_i}} \ + + b_{\chi_{a_i},\Theta}\Theta \ + + b_{\chi_{a_i},\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j (j \ne i)} b_{\chi_{a_i},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\chi_{a_i},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling active scalar i to its gradient and/or to other scalar fields including temperature/entropy, and finally + +.. math:: + :label: coupled_bcs_dchiadr + + \frac{\partial \chi_{a_i}}{\partial r} = b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}} \ + + b_{\frac{\partial\chi_{a_i}}{\partial r},\Theta}\Theta \ + + b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j} b_{\frac{\partial\chi_{a_i}}{\partial r},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j (j \ne i)} b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling the gradient of active scalar i to its value and/or to other scalar fields including temperature/entropy. + +The values of the :math:`b_{i,i}` terms can be constant or spatially varying. The values of the :math:`b_{i,j}` coefficients are constant. All values can be set using options in the boundary conditions namespace. + +Setup +^^^^^ + +Boundary conditions +******************* + +Model parameters for the scalar fields follow the same convention as temperature but using the prefix `chi_a` or `chi_p` for active and passive +scalars respectively. + +**couple_tvar_top** + Logical flag indicating whether thermal variable (T,S) should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_tvar_bottom** + Logical flag indicating whether thermal variable (T,S) should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_dtdr_top** + Logical flag indicating whether radial gradient of thermal variable (T,S) should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_dtdr_bottom** + Logical flag indicating whether radial gradient of thermal variable (T,S) should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_chivar_a_top** + Logical flag indicating whether active scalar i should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_chivar_a_bottom** + Logical flag indicating whether active scalar i should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_dchidr_a_top** + Logical flag indicating whether radial gradient of active scalar i should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_dchidr_a_bottom** + Logical flag indicating whether radial gradient of active scalar i should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**T_top** + Set the thermal variable, :math:`b_{\Theta,\Theta}`, at the top of the domain (overloaded) +**T_top_file** + Set a spatially varying thermal variable, :math:`b_{\Theta,\Theta}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**T_dTdr_coeff_top** + Set the coupling coefficient between the thermal variable and the radial derivative of the thermal variable, :math:`b_{\Theta,\frac{\partial \Theta}{\partial r}}`, at the top of the domain +**T_chi_a_coeff_top(i)** + Set the coupling coefficient between the thermal variable and active scalar field i, :math:`b_{\Theta,\chi_{a_i}}`, at the top of the domain +**T_dchidr_a_coeff_top(i)** + Set the coupling coefficient between the thermal variable and the radial derivative of active scalar field i, :math:`b_{\Theta,\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain + +**T_bottom** + Set the thermal variable, :math:`b_{\Theta,\Theta}`, at the base of the domain (overloaded) +**T_bottom_file** + Set a spatially varying thermal variable, :math:`b_{\Theta,\Theta}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**T_dTdr_coeff_bottom** + Set the coupling coefficient between the thermal variable and the radial derivative of the thermal variable, :math:`b_{\Theta,\frac{\partial \Theta}{\partial r}}`, at the base of the domain +**T_chi_a_coeff_bottom(i)** + Set the coupling coefficient between the thermal variable and active scalar field i, :math:`b_{\Theta,\chi_{a_i}}`, at the base of the domain +**T_dchidr_a_coeff_bottom(i)** + Set the coupling coefficient between the thermal variable and the radial derivative of active scalar field i, :math:`b_{\Theta,\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain + +**dTdr_top** + Set the radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain (overloaded) +**dTdr_top_file** + Set a spatially varying radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**dTdr_T_coeff_top** + Set the coupling coefficient between the radial derivative of the thermal variable and the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\Theta}`, at the top of the domain +**dTdr_chi_a_coeff_top(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain +**dTdr_dchidr_a_coeff_top(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain + +**dTdr_bottom** + Set the radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain (overloaded) +**dTdr_bottom_file** + Set a spatially varying radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**dTdr_T_coeff_bottom** + Set the coupling coefficient between the radial derivative of the thermal variable and the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\Theta}`, at the base of the domain +**dTdr_chi_a_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain +**dTdr_dchidr_a_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain + + +**chi_a_top(i)** + Set active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the top of the domain (overloaded) +**chi_a_top_file(i)** + Set a spatially varying active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**chi_a_T_coeff_top(i)** + Set the coupling coefficient between active scalar field i and the thermal variable, :math:`b_{\chi_{a_i},\Theta}`, at the top of the domain +**chi_a_dTdr_coeff_top(i)** + Set the coupling coefficient between active scalar field i and the radial derivative of the thermal variable, :math:`b_{\chi_{a_i},\frac{\partial\Theta}{\partial r}}`, at the top of the domain +**chi_a_chi_a_coeff_top(i,j)** + Set the coupling coefficient between active scalar field i and active scalar field j, :math:`b_{\chi_{a_i},\chi_{a_j}}`, at the top of the domain. :math:`b_{\chi_{a_i},\chi_{a_i}}` is ignored +**chi_a_dchidr_a_coeff_top(i,j)** + Set the coupling coefficient between active scalar field i and the radial derivative of active scalar field j, :math:`b_{\chi_{a_i},\frac{\partial \chi_{a_j}}{\partial r}}`, at the top of the domain + +**chi_a_bottom(i)** + Set active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the base of the domain (overloaded) +**chi_a_bottom_file(i)** + Set a spatially varying active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**chi_a_T_coeff_bottom(i)** + Set the coupling coefficient between active scalar field i and the thermal variable, :math:`b_{\chi_{a_i},\Theta}`, at the base of the domain +**chi_a_dTdr_coeff_bottom(i)** + Set the coupling coefficient between active scalar field i and the radial derivative of the thermal variable, :math:`b_{\chi_{a_i},\frac{\partial\Theta}{\partial r}}`, at the base of the domain +**chi_a_chi_a_coeff_bottom(i,j)** + Set the coupling coefficient between active scalar field i and active scalar field j, :math:`b_{\chi_{a_i},\chi_{a_j}}`, at the base of the domain. :math:`b_{\chi_{a_i},\chi_{a_i}}` is ignored +**chi_a_dchidr_a_coeff_bottom(i,j)** + Set the coupling coefficient between active scalar field i and the radial derivative of active scalar field j, :math:`b_{\chi_{a_i},\frac{\partial \chi_{a_j}}{\partial r}}`, at the base of the domain + +**dchidr_a_top(i)** + Set the radial derivative of active scalar field i, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain (overloaded) +**dchidr_a_top_file(i)** + Set a spatially varying radial derivative of active scalar field i, :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**dchidr_a_T_coeff_top(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\Theta}`, at the top of the domain +**dchidr_a_dTdr_coeff_top(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain +**dchidr_a_chi_a_coeff_top(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\chi_{a_j}}`, at the top of the domain +**dchidr_a_dchidr_a_coeff_top(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_j}}{\partial r}}`, at the top of the domain. :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}` is ignored + +**dchidr_a_bottom(i)** + Set the radial derivative of active scalar field i, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain (overloaded) +**dchidr_a_bottom_file(i)** + Set a spatially varying radial derivative of active scalar field i, :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**dchidr_a_T_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\Theta}`, at the base of the domain +**dchidr_a_dTdr_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain +**dchidr_a_chi_a_coeff_bottom(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\chi_{a_j}}`, at the base of the domain +**dchidr_a_dchidr_a_coeff_bottom(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_j}}{\partial r}}`, at the base of the domain. :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}` is ignored + +Further Information +^^^^^^^^^^^^^^^^^^^ + +See `tests/coupled_bcs` for example input files. + + + diff --git a/doc/source/User_Guide/under_development.rst b/doc/source/User_Guide/under_development.rst index 8c4e31d4..2ae76b40 100644 --- a/doc/source/User_Guide/under_development.rst +++ b/doc/source/User_Guide/under_development.rst @@ -17,6 +17,14 @@ Arbitrary Scalar Fields .. include:: arbitrary_scalar_fields.txt +.. _coupled_bcs: + +Coupled Boundary Conditions +--------------------------- + +.. include:: coupled_bcs.txt + + diff --git a/src/IO/Spherical_IO.F90 b/src/IO/Spherical_IO.F90 index 51bb7d88..bb52c65b 100755 --- a/src/IO/Spherical_IO.F90 +++ b/src/IO/Spherical_IO.F90 @@ -663,12 +663,12 @@ Subroutine Write_Full_3D(qty) Implicit None Real*8, Intent(In) :: qty(:,my_rmin:,my_theta_min:) Integer :: i, funit - Character*4 :: qstring + Character*5 :: qstring Character*120 :: iterstring, data_file, grid_file ! Write the data file (Parallel I/O) Write(iterstring,i_ofmt) current_iteration - Write(qstring,'(i4.4)') current_qval + Write(qstring,'(i5.5)') current_qval data_file = trim(local_file_path)//'Spherical_3D/'//trim(iterstring)//'_'//qstring Call full_3d_buffer%cache_data(qty) Call full_3d_buffer%write_data(filename=data_file) diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index d242ac9c..7a19ed83 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -33,18 +33,28 @@ Module BoundaryConditions Implicit None - Logical :: Fix_Tvar_Top = .True. - Logical :: Fix_Tvar_Bottom = .True. + Logical :: Fix_Tvar_Top = .False. + Logical :: Fix_Tvar_Bottom = .False. Logical :: Fix_dTdr_Top = .False. Logical :: Fix_dTdr_Bottom = .False. - Logical :: Fix_chivar_a_Top(1:n_scalar_max) = .True. - Logical :: Fix_chivar_a_Bottom(1:n_scalar_max) = .True. + Logical :: Couple_Tvar_Top = .False. + Logical :: Couple_Tvar_Bottom = .False. + Logical :: Couple_dTdr_Top = .False. + Logical :: Couple_dTdr_Bottom = .False. + + Logical :: Fix_chivar_a_Top(1:n_scalar_max) = .False. + Logical :: Fix_chivar_a_Bottom(1:n_scalar_max) = .False. Logical :: Fix_dchidr_a_Top(1:n_scalar_max) = .False. Logical :: Fix_dchidr_a_Bottom(1:n_scalar_max) = .False. - Logical :: Fix_chivar_p_Top(1:n_scalar_max) = .True. - Logical :: Fix_chivar_p_Bottom(1:n_scalar_max) = .True. + Logical :: Couple_chivar_a_Top(1:n_scalar_max) = .False. + Logical :: Couple_chivar_a_Bottom(1:n_scalar_max) = .False. + Logical :: Couple_dchidr_a_Top(1:n_scalar_max) = .False. + Logical :: Couple_dchidr_a_Bottom(1:n_scalar_max) = .False. + + Logical :: Fix_chivar_p_Top(1:n_scalar_max) = .False. + Logical :: Fix_chivar_p_Bottom(1:n_scalar_max) = .False. Logical :: Fix_dchidr_p_Top(1:n_scalar_max) = .False. Logical :: Fix_dchidr_p_Bottom(1:n_scalar_max) = .False. @@ -58,16 +68,46 @@ Module BoundaryConditions Logical :: Dipole_Field_Bottom = .False. Logical :: adjust_dTdr_Top = .True. - Real*8 :: T_Bottom = 1.0d0 + Real*8 :: T_Bottom = 0.0d0 Real*8 :: T_Top = 0.0d0 Real*8 :: dTdr_Top = 0.0d0 Real*8 :: dTdr_Bottom = 0.0d0 - Real*8 :: chi_a_Bottom(1:n_scalar_max) = 1.0d0 + Real*8 :: T_chi_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dTdr_coeff_top = 0.0d0 + Real*8 :: T_chi_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dTdr_coeff_bottom = 0.0d0 + Real*8 :: dTdr_chi_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_T_coeff_top = 0.0d0 + Real*8 :: dTdr_chi_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_T_coeff_bottom = 0.0d0 + + Real*8 :: chi_a_Bottom(1:n_scalar_max) = 0.0d0 Real*8 :: chi_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_chi_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dchidr_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_chi_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dchidr_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_chi_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dchidr_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_chi_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dchidr_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_p_Bottom(1:n_scalar_max) = 1.0d0 Real*8 :: chi_p_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_p_Top(1:n_scalar_max) = 0.0d0 @@ -120,7 +160,17 @@ Module BoundaryConditions fix_dchidr_a_bottom, fix_dchidr_a_top, dchidr_a_top, dchidr_a_bottom, & chi_p_top_file, chi_p_bottom_file, dchidr_p_top_file, dchidr_p_bottom_file, & fix_chivar_p_top, fix_chivar_p_bottom, chi_p_bottom, chi_p_top, & - fix_dchidr_p_bottom, fix_dchidr_p_top, dchidr_p_top, dchidr_p_bottom + fix_dchidr_p_bottom, fix_dchidr_p_top, dchidr_p_top, dchidr_p_bottom, & + couple_Tvar_top, couple_Tvar_bottom, couple_dtdr_top, couple_dtdr_bottom, & + T_chi_a_coeff_top, T_dchidr_a_coeff_top, T_dTdr_coeff_top, & + T_chi_a_coeff_bottom, T_dchidr_a_coeff_bottom, T_dTdr_coeff_bottom, & + dTdr_chi_a_coeff_top, dTdr_dchidr_a_coeff_top, dTdr_T_coeff_top, & + dTdr_chi_a_coeff_bottom, dTdr_dchidr_a_coeff_bottom, dTdr_T_coeff_bottom, & + couple_chivar_a_top, couple_chivar_a_bottom, couple_dchidr_a_top, couple_dchidr_a_bottom, & + chi_a_chi_a_coeff_top, chi_a_dchidr_a_coeff_top, chi_a_T_coeff_top, chi_a_dTdr_coeff_top, & + chi_a_chi_a_coeff_bottom, chi_a_dchidr_a_coeff_bottom, chi_a_T_coeff_bottom, chi_a_dTdr_coeff_bottom, & + dchidr_a_chi_a_coeff_top, dchidr_a_dchidr_a_coeff_top, dchidr_a_T_coeff_top, dchidr_a_dTdr_coeff_top, & + dchidr_a_chi_a_coeff_bottom, dchidr_a_dchidr_a_coeff_bottom, dchidr_a_T_coeff_bottom, dchidr_a_dTdr_coeff_bottom Contains @@ -129,18 +179,85 @@ Subroutine Initialize_Boundary_Conditions() Real*8 :: tilt_angle_radians,a,b Real*8 :: Ftop, rhotk_top, Fbottom Integer :: i - - fix_tvar_top = .not. fix_dtdr_top - fix_tvar_bottom = .not. fix_dtdr_bottom + Integer :: n_active_bcs + + n_active_bcs = count( (/ fix_tvar_top, fix_dtdr_top, couple_tvar_top, couple_dtdr_top /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for tvar on top boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for tvar on top boundary.") + call stdout%print(" Defaulting to fix_tvar_top with T_top = 0.0.") + call stdout%print(" ") + fix_tvar_top = .true. + T_top = 0.0 + endif + + n_active_bcs = count( (/ fix_tvar_bottom, fix_dtdr_bottom, couple_tvar_bottom, couple_dtdr_bottom /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for tvar on bottom boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for tvar on bottom boundary.") + call stdout%print(" Defaulting to fix_tvar_bottom with T_bottom = 1.0.") + call stdout%print(" ") + fix_tvar_bottom = .true. + T_bottom = 1.0 + endif do i = 1, n_active_scalars - fix_chivar_a_top(i) = .not. fix_dchidr_a_top(i) - fix_chivar_a_bottom(i) = .not. fix_dchidr_a_bottom(i) + n_active_bcs = count( (/ fix_chivar_a_top(i), fix_dchidr_a_top(i), couple_chivar_a_top(i), couple_dchidr_a_top(i) /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on top boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on top boundary.") + call stdout%print(" Defaulting to fix_chivar_a_top.") + call stdout%print(" ") + fix_chivar_a_top(i) = .true. + endif + + n_active_bcs = count( (/ fix_chivar_a_bottom(i), fix_dchidr_a_bottom(i), & + couple_chivar_a_bottom(i), couple_dchidr_a_bottom(i) /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on bottom boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on bottom boundary.") + call stdout%print(" Defaulting to fix_chivar_a_bottom.") + call stdout%print(" ") + fix_chivar_a_bottom(i) = .true. + endif end do do i = 1, n_passive_scalars - fix_chivar_p_top(i) = .not. fix_dchidr_p_top(i) - fix_chivar_p_bottom(i) = .not. fix_dchidr_p_bottom(i) + n_active_bcs = count( (/ fix_chivar_p_top(i), fix_dchidr_p_top(i) /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on top boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on top boundary.") + call stdout%print(" Defaulting to fix_chivar_p_top.") + call stdout%print(" ") + fix_chivar_p_top(i) = .true. + endif + + n_active_bcs = count( (/ fix_chivar_p_bottom(i), fix_dchidr_p_bottom(i) /) ) + if (n_active_bcs > 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on bottom boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on bottom boundary.") + call stdout%print(" Defaulting to fix_chivar_p_bottom.") + call stdout%print(" ") + fix_chivar_p_bottom(i) = .true. + endif end do If (no_slip_top .and. strict_L_Conservation) Then @@ -265,7 +382,7 @@ Subroutine Generate_Boundary_Mask() If (.not. full_restart) Then - If (fix_tvar_top) Then + If (fix_tvar_top .or. couple_tvar_top) Then if (trim(T_top_file) .eq. '__nothing__') then bc_val= T_Top*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, uind) @@ -274,7 +391,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_tvar_bottom) Then + If (fix_tvar_bottom .or. couple_tvar_bottom) Then if (trim(T_bottom_file) .eq. '__nothing__') then bc_val= T_bottom*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, lind) @@ -283,7 +400,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_dtdr_top) Then + If (fix_dtdr_top .or. couple_dtdr_top) Then if (trim(dTdr_top_file) .eq. '__nothing__') then bc_val= dtdr_top*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, uind) @@ -292,7 +409,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_dtdr_bottom) Then + If (fix_dtdr_bottom .or. couple_dtdr_bottom) Then if (trim(dTdr_bottom_file) .eq. '__nothing__') then bc_val= dtdr_bottom*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, lind) @@ -322,7 +439,7 @@ Subroutine Generate_Boundary_Mask() Endif do i = 1, n_active_scalars - If (fix_chivar_a_top(i)) Then + If (fix_chivar_a_top(i) .or. couple_chivar_a_top(i)) Then if (trim(chi_a_top_file(i)) .eq. '__nothing__') then bc_val= chi_a_Top(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, uind) @@ -344,7 +461,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (fix_chivar_a_bottom(i)) Then + If (fix_chivar_a_bottom(i) .or. couple_chivar_a_bottom(i)) Then if (trim(chi_a_bottom_file(i)) .eq. '__nothing__') then bc_val= chi_a_bottom(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, lind) @@ -366,7 +483,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (Fix_dchidr_a_Top(i)) Then + If (Fix_dchidr_a_Top(i) .or. couple_dchidr_a_top(i)) Then if (trim(dchidr_a_top_file(i)) .eq. '__nothing__') then bc_val= dchidr_a_top(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, uind) @@ -388,7 +505,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (Fix_dchidr_a_Bottom(i)) Then + If (Fix_dchidr_a_Bottom(i) .or. couple_dchidr_a_bottom(i)) Then if (trim(dchidr_a_bottom_file(i)) .eq. '__nothing__') then bc_val= dchidr_a_bottom(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, lind) @@ -496,8 +613,8 @@ End Subroutine Transport_Dependencies Subroutine Restore_BoundaryCondition_Defaults() Implicit None - Fix_Tvar_Top = .True. - Fix_Tvar_Bottom = .True. + Fix_Tvar_Top = .False. + Fix_Tvar_Bottom = .False. Fix_dTdr_Top = .False. Fix_dTdr_Bottom = .False. Fix_divrt_top = .False. @@ -508,7 +625,7 @@ Subroutine Restore_BoundaryCondition_Defaults() Fix_poloidalfield_bottom = .False. Impose_Dipole_Field = .False. - T_Bottom = 1.0d0 + T_Bottom = 0.0d0 T_Top = 0.0d0 dTdr_Top = 0.0d0 dTdr_Bottom = 0.0d0 diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index f8b37e7c..fb35ad02 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -521,7 +521,7 @@ Subroutine Set_Boundary_Conditions(mode_ind) Implicit None Real*8 :: samp,one Integer, Intent(In) :: mode_ind - Integer :: l, r,lp, i + Integer :: l, r,lp, i, j one = 1.0d0 lp = mode_ind @@ -558,6 +558,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_top) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_top) Then + do i = 1, n_active_scalars + samp = -T_chi_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_top) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif r = N_R If (fix_tvar_bottom) Then @@ -566,8 +590,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_bottom) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_bottom) Then + do i = 1, n_active_scalars + samp = -T_chi_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_bottom) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif - ! PASSIVE + ! active scalars do i = 1, n_active_scalars r = 1 If (fix_chivar_a_top(i)) Then @@ -576,6 +624,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_top(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_top(i)) Then + do j = 1, n_active_scalars + samp = -chi_a_chi_a_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_a_dchidr_a_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_a_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_a_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_top(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_a_chi_a_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_a_dchidr_a_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_a_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_a_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif r = N_R If (fix_chivar_a_bottom(i)) Then @@ -584,8 +658,35 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_bottom(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -chi_a_chi_a_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_a_dchidr_a_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_a_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_a_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_a_chi_a_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_a_dchidr_a_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_a_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_a_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif end do + ! passive scalars do i = 1, n_passive_scalars r = 1 If (fix_chivar_p_top(i)) Then @@ -663,6 +764,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_top) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_top) Then + do i = 1, n_active_scalars + samp = -T_chi_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_top) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_a_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif r = N_R If (fix_tvar_bottom) Then @@ -671,6 +796,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_bottom) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_bottom) Then + do i = 1, n_active_scalars + samp = -T_chi_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_bottom) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_a_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif ! chivar BC do i = 1, n_active_scalars @@ -681,6 +830,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_top(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_top(i)) Then + do j = 1, n_active_scalars + samp = -chi_a_chi_a_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_a_dchidr_a_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_a_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_a_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_top(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_a_chi_a_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_a_dchidr_a_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_a_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_a_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif r = N_R If (fix_chivar_a_bottom(i)) Then @@ -689,6 +864,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_bottom(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -chi_a_chi_a_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_a_dchidr_a_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_a_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_a_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_a_chi_a_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_a_dchidr_a_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_a_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_a_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif end do do i = 1, n_passive_scalars diff --git a/tests/coupled_bcs/const/generate_input.py b/tests/coupled_bcs/const/generate_input.py new file mode 100644 index 00000000..c85ad352 --- /dev/null +++ b/tests/coupled_bcs/const/generate_input.py @@ -0,0 +1,14 @@ +### Use rayleigh_spectral_input.py to generate generic input. + +from rayleigh_spectral_input import * + +rmin = 0.5; rmax = 1.0 + +si = SpectralInput(n_theta=1,n_r=48) +si.transform_from_rtp_function(lambda radius: 2.0 - (rmax/radius)*(radius - rmin)/(rmax - rmin), rmin=rmin, rmax=rmax) +si.write('radial_prof_init') + +si = SpectralInput() +si.add_mode(0.0+0.j, 0, 0, 0) +si.write('constant_init') + diff --git a/tests/coupled_bcs/const/main_input b/tests/coupled_bcs/const/main_input new file mode 100644 index 00000000..988215e0 --- /dev/null +++ b/tests/coupled_bcs/const/main_input @@ -0,0 +1,88 @@ +! This is a smaller version of the input_examples/benchmark_diagnostics_input +! model that only runs for two time-steps so that we can check the boundary conditions. + +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + rmin = 0.5 + rmax = 1.0 +/ +&numerical_controls_namelist + chebyshev = .true. +/ +&physical_controls_namelist + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 2 +/ +&temporal_controls_namelist + max_time_step = 1.0d-4 + max_iterations = 2 + alpha_implicit = 0.50001d0 + checkpoint_interval = 500000 + quicksave_interval = 1000000 + num_quicksaves = -1 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +full3d_values = 501 507 10001 10004 10201 10204 +full3d_frequency = 2 + +shellavg_values = 501 507 10001 10004 10201 10204 +shellavg_frequency = 2 +shellavg_nrec = 1 + +shellslice_levels = 1,48 +shellslice_values = 501 507 10001 10004 10201 10204 ! T dTdr chi1 dchi1dr chi2 dchi2dr +shellslice_frequency = 2 +shellslice_nrec = 1 +/ + +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +dTdr_Top = 32.31 +T_Bottom = 17.865 +couple_dtdr_top = .true. +dTdr_chi_a_coeff_top(1) = 70.2 +couple_tvar_bottom = .true. +T_chi_a_coeff_bottom(1) = 34.1 +T_dchidr_a_coeff_bottom(1) = -51.87 +chi_a_Top(1) = 1.0d0 +chi_a_Bottom(1) = 2.0d0 +fix_chivar_a_top(1) = .true. +fix_chivar_a_bottom(1) = .true. +chi_a_Top(2) = 1.0342d0 +dchidr_a_Bottom(2) = -27.6d0 +couple_chivar_a_top(2) = .true. +chi_a_chi_a_coeff_top(2,1) = 107.01 +chi_a_dchidr_a_coeff_top(2,1) = 7.064 +couple_dchidr_a_bottom(2) = .true. +dchidr_a_chi_a_coeff_bottom(2,1) = 9.872 +/ +&Initial_Conditions_Namelist +init_type = 8 ! File init +t_init_file = 'constant_init' +chi_a_init_file(1) = 'radial_prof_init' +chi_a_init_file(2) = 'constant_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 0.0 +Prandtl_Number = 1.0d0 +Magnetic_Prandtl_Number = 5.0d0 +reference_type = 1 +heating_type = 0 ! No heating +gravity_power = 1.0d0 ! g ~ radius +/ +&Transport_Namelist +/ diff --git a/tests/coupled_bcs/run_test.sh b/tests/coupled_bcs/run_test.sh new file mode 100755 index 00000000..40e09c34 --- /dev/null +++ b/tests/coupled_bcs/run_test.sh @@ -0,0 +1,13 @@ +#!/usr/bin/env bash + +cd tests/coupled_bcs + +# we run a const case, where the bcs on the field that everything is coupled to are just constant +cd const +PYTHONPATH=../../../pre_processing:$PYTHONPATH python3 generate_input.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + +# after running, we test the output for errors +PYTHONPATH=../../post_processing:../../pre_processing:$PYTHONPATH python3 test_output.py + diff --git a/tests/coupled_bcs/test_output.py b/tests/coupled_bcs/test_output.py new file mode 100644 index 00000000..ea995193 --- /dev/null +++ b/tests/coupled_bcs/test_output.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 + +from rayleigh_diagnostics import Shell_Slices, main_input +import numpy as np +import sys, os + +error = False + +def test_bcs(case): + """ + For a given case check if the values contained in the shell slices match the expected values based on the coefficients in + main_input + """ + error = False + + mi = main_input(os.path.join(case, 'main_input')) + bcs = mi.vals['boundary_conditions'] + slices = Shell_Slices('00000002', path=os.path.join(case, 'Shell_Slices/')) + + # values produced by the code + chi1_rmin = slices.vals[:,:,1,slices.lut[10001],0] + chi1_rmax = slices.vals[:,:,0,slices.lut[10001],0] + T_rmin = slices.vals[:,:,1,slices.lut[501],0] + T_rmax = slices.vals[:,:,0,slices.lut[501],0] + chi2_rmin = slices.vals[:,:,1,slices.lut[10201],0] + chi2_rmax = slices.vals[:,:,0,slices.lut[10201],0] + + # derivatives produced by the code + dchi1dr_rmin = slices.vals[:,:,1,slices.lut[10004],0] + dchi1dr_rmax = slices.vals[:,:,0,slices.lut[10004],0] + dTdr_rmin = slices.vals[:,:,1,slices.lut[507],0] + dTdr_rmax = slices.vals[:,:,0,slices.lut[507],0] + dchi2dr_rmin = slices.vals[:,:,1,slices.lut[10204],0] + dchi2dr_rmax = slices.vals[:,:,0,slices.lut[10204],0] + + # coefficients input + chi1bc_rmax = float(bcs['chi_a_top(1)'.lower()].replace('d', 'e')) + chi1bc_rmin = float(bcs['chi_a_bottom(1)'.lower()].replace('d', 'e')) + + dTdrbc_rmax = float(bcs['dTdr_top'.lower()].replace('d', 'e')) + dtdr_chi1_coeff_rmax = float(bcs['dTdr_chi_a_coeff_top(1)'.lower()].replace('d', 'e')) + + Tbc_rmin = float(bcs['T_bottom'.lower()].replace('d', 'e')) + t_chi1_coeff_rmin = float(bcs['T_chi_a_coeff_bottom(1)'.lower()].replace('d', 'e')) + t_dchi1dr_coeff_rmin = float(bcs['T_dchidr_a_coeff_bottom(1)'.lower()].replace('d', 'e')) + + chi2bc_rmax = float(bcs['chi_a_top(2)'.lower()].replace('d', 'e')) + chi2_chi1_coeff_rmax = float(bcs['chi_a_chi_a_coeff_top(2,1)'.lower()].replace('d', 'e')) + chi2_dchi1dr_coeff_rmax = float(bcs['chi_a_dchidr_a_coeff_top(2,1)'.lower()].replace('d', 'e')) + + dchi2drbc_rmin = float(bcs['dchidr_a_bottom(2)'.lower()].replace('d', 'e')) + dchi2dr_chi1_coeff_rmin = float(bcs['dchidr_a_chi_a_coeff_bottom(2,1)'.lower()].replace('d', 'e')) + + # check if values produced match expected values based on input coefficients + def failed(values, expected, tol=1.e-10): + return np.any(np.abs(values - expected) > tol) + + error = error or failed(chi1_rmax, chi1bc_rmax) + error = error or failed(chi1_rmin, chi1bc_rmin) + + error = error or failed(dTdr_rmax, dtdr_chi1_coeff_rmax*chi1_rmax + dTdrbc_rmax) + error = error or failed(T_rmin, t_chi1_coeff_rmin*chi1_rmin + t_dchi1dr_coeff_rmin*dchi1dr_rmin + Tbc_rmin) + + error = error or failed(chi2_rmax, chi2_chi1_coeff_rmax*chi1_rmax + chi2_dchi1dr_coeff_rmax*dchi1dr_rmax + chi2bc_rmax) + error = error or failed(dchi2dr_rmin, dchi2dr_chi1_coeff_rmin*chi1_rmin + dchi2drbc_rmin) + + return error + +lerror = test_bcs('const') +if lerror: + print("ERROR: const case produced unexpected output relative to input coefficients!") +error = error or lerror + +if error: + sys.exit(1) +sys.exit(0) +