test_stardis_smeteo_ground_temperature.sh (6341B)
1 #!/bin/sh 2 3 # Copyright (C) 2025 |Méso|Star> (contact@meso-star.com) 4 # 5 # This program is free software: you can redismeteobute it and/or modify 6 # it under the terms of the GNU General Public License as published by 7 # the Free Software Foundation, either version 3 of the License, or 8 # (at your option) any later version. 9 # 10 # This program is dismeteobuted in the hope that it will be useful, 11 # but WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 # GNU General Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License 16 # along with this program. If not, see <http://www.gnu.org/licenses/>. 17 18 set -e 19 20 if [ $# -lt 1 ]; then 21 >&2 printf 'usage: %s smeteo-file\n' "${0##*/}" 22 exit 1 23 fi 24 25 smeteo="$1" 26 27 ######################################################################## 28 # Ground geometry 29 ######################################################################## 30 bound=50000 # [m] 31 depth=8 # [m] 32 33 x="-${bound}"; y="-${bound}"; z="-${depth}" 34 X="+${bound}"; Y="+${bound}"; Z="0" 35 36 vertices="\ 37 ${x} ${y} ${z} 38 ${X} ${y} ${z} 39 ${X} ${Y} ${z} 40 ${x} ${Y} ${z} 41 ${x} ${y} ${Z} 42 ${X} ${y} ${Z} 43 ${X} ${Y} ${Z} 44 ${x} ${Y} ${Z}" 45 46 indices="\ 47 1 4 2 48 2 4 3 49 8 5 6 50 8 6 7 51 1 2 5 52 5 2 6 53 6 2 3 54 6 3 7 55 7 3 8 56 8 3 4 57 8 4 5 58 5 4 1" 59 60 ######################################################################## 61 # Helper functions 62 ######################################################################## 63 facet() # "index0 index1 index2" 64 { 65 printf '\tfacet normal 0 0 1\n' 66 printf '\t\touter loop\n' 67 68 printf '%s\n' "$1" \ 69 | sed -e 's/^[[:space:]]\{1,\}//g' \ 70 -e 's/[[:space:]]\{1,\}$//g' \ 71 -e 's/[[:space:]]\{1,\}/\n/g' \ 72 | while read -r i; do 73 printf '\t\t\tvertex ' 74 printf '%s\n' "${vertices}" | sed -n "${i}p" 75 done 76 77 printf '\t\tendloop\n' 78 printf '\tendfacet\n' 79 } 80 81 ground() # list of facets id 82 { 83 84 if [ $# -eq 0 ]; then 85 facets="1,\$p" 86 else 87 facets="$(printf '%s ' "$@" | sed \ 88 -e 's/^[[:space:]]\{1,\}//g' \ 89 -e 's/[[:space:]]\{1,\}$/p/g' \ 90 -e 's/[[:space:]]\{1,\}/p;/g')" 91 fi 92 93 printf 'solid ground\n' 94 95 printf '%s\n' "${indices}" \ 96 | sed -n "${facets}" \ 97 | while read -r facet; do 98 facet "${facet}" 99 done 100 101 printf 'endsolid\n' 102 } 103 104 stardis_input_imposed_flux() 105 { 106 # The plugin 107 printf 'PROGRAM Meteo libstardis_smeteo.so %s\n' "${smeteo}" 108 echo '' 109 110 # Media 111 echo 'SOLID ground 1 1500 1500 0.02 281.85 UNKNOWN 0 FRONT ground.stl' 112 echo '' 113 114 # Limit condition 115 echo 'H_BOUNDARY_FOR_SOLID adiabatic 0 0 0 0 0 ground_xXyY.stl' 116 echo 'T_BOUNDARY_FOR_SOLID underground 284.0 ground_z.stl' 117 echo 'HF_BOUNDARY_FOR_SOLID_PROG atm Meteo ground_Z.stl PROG_PARAMS -hls' 118 echo 'TRAD_PROG Meteo' 119 } 120 121 stardis_input_external_source() 122 { 123 # The plugin 124 printf 'PROGRAM Meteo libstardis_smeteo.so %s PROG_PARAMS -a meeus\n' "${smeteo}" 125 echo '' 126 127 # Media 128 echo 'SOLID ground 1 1500 1500 0.02 281.85 UNKNOWN 0 FRONT ground.stl' 129 echo 'FLUID_PROG sky Meteo BACK ground.stl' 130 echo '' 131 132 # Connection 133 echo 'F_SOLID_FLUID_CONNECTION_PROG ground_sky_connect Meteo ground_Z.stl PROG_PARAMS -hl' 134 echo 'SOLID_FLUID_CONNECTION adiabatic 300 0 0 0 ground_xXyY.stl' 135 echo '' 136 137 # Limit condition 138 echo 'SOLID_FLUID_CONNECTION_PROG underground Meteo ground_z.stl PROG_PARAMS -t 284.0' 139 echo '' 140 echo 'SPHERICAL_SOURCE_PROG 6.987e8 Meteo' 141 echo '' 142 echo 'TRAD_PROG Meteo' 143 } 144 145 run_stardis() # model, time 146 { 147 stardis -a wos -M "$1" -s ground_Z.stl,"$2" 148 } 149 150 ######################################################################## 151 # The test 152 ######################################################################## 153 # Generate geometry 154 ground > ground.stl 155 ground 1 2 > ground_z.stl 156 ground 3 4 > ground_Z.stl 157 ground 5 6 7 8 9 10 11 12 > ground_xXyY.stl 158 159 # Generate the Stardis input file 160 stardis_input_imposed_flux > stardis_model_imposed_flux.txt 161 stardis_input_external_source > stardis_model_external_source.txt 162 163 # Use Stardis to calculate the ground temperature 164 date="30-AUG-1850 01:30:00" # 1st date on which temperature is calculated 165 ndates="32" # Overall number of dates 166 167 # Check that at least the 1st date exists in the meteorlogical file 168 if ! grep -qe "^${date}" "${smeteo}"; then 169 >&2 printf '%s: unable to find the date "%s" in "%s" file\n' \ 170 "${0##*/}" "${date}" "${smeteo}" 171 exit 1 172 fi 173 174 # Define the sed commands to extract the ${ndates} dates from the 175 # meteorological file 176 i=1; # The corresponding date... 177 n=""; 178 while [ "${i}" -lt "${ndates}" ]; do # ... followed by ${ndates}-1 dates 179 n="${n}N;" 180 i="$((i+1))" 181 done 182 183 isimu=0; 184 185 # Run the simulations 186 >&2 printf 'Running 0 %%\n' 187 sed -n "/^${date}/{${n}p;}" "${smeteo}" \ 188 | while read -r i; do 189 190 # Remove duplicate spaces in the line 191 entry="$(echo "${i}" |sed 's/[[:space:]]\{1,\}/ /g')" 192 193 # Obtain the surface temperature [K]. This will serve as a reference 194 # against which the Monte Carlo estimate of the ground temperature 195 # will be compared. 196 Tsrf="$(echo "${entry}" | cut -d' ' -f3)" 197 198 # Get the date 199 date="$(echo "${entry}" | cut -d' ' -f1-2)" 200 201 # Get the time in seconds since January 1, 1850 202 time="$(echo "${entry}" | cut -d' ' -f14)" # [day/1850] 203 time="$(echo "${time} * 24*3600" | bc -l | cut -d'.' -f1)" # [s] 204 205 # Run stardis to retrieve the expected value and standard error of the 206 # ground temperature estimate. 207 mc0="$(run_stardis stardis_model_imposed_flux.txt "${time}" \ 208 | cut -d' ' -f1-2)" 209 mc1="$(run_stardis stardis_model_external_source.txt "${time}" \ 210 | cut -d' ' -f1-2)" 211 212 # Print the calculation result 213 printf '%s %s %s %s\n' "${date}" "${Tsrf}" "${mc0}" "${mc1}" 214 215 # Print progression bar 216 isimu="$((isimu+1))" 217 pcent="$(echo "${isimu} / ${ndates} * 100" | bc -l | cut -d'.' -f1)" 218 >&2 printf 'Running %3d %%\n' "${pcent}" 219 done 220 221 # TODO: check the Stardis result against a reference, most likely the 222 # surface temperature provided by the Meteorological file at the time of 223 # observation. The test has not yet been performed because the system 224 # parameters may not correspond to the system used to generate the input 225 # meteorlogical data. For example the thermophysical properties of the 226 # ground still need to be verified to ensure that they match those used 227 # to calculate the surface temperature in the smeteo file.