hedge.timestep.multirate

Source Code for Module hedge.timestep.multirate_ab.methods

1 """Multirate-AB ODE solver.""" 2 3 from __future__ import division 4 5 __copyright__ = "Copyright (C) 2009 Andreas Stock, Andreas Kloeckner" 6 7 __license__ = """ 8 This program is free software: you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation, either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see U{http://www.gnu.org/licenses/}. 20 """ 21 22 23 24 25 from pytools import Record 26 from pymbolic import var 27 28 29 30 31 # symbols --------------------------------------------------------------------- 32 # components

33 -class CO_FAST: pass

34 -class CO_SLOW: pass

35 36 # histories:

37 -class HIST_F2F: pass

38 -class HIST_F2S: pass

39 -class HIST_S2F: pass

40 -class HIST_S2S: pass

41 42 HIST_NAMES = [HIST_F2F, HIST_S2F, HIST_F2S, HIST_S2S] 43 44 45 46 47 # actual method descriptions -------------------------------------------------- 48 n = substep_count = var("substep_count") 49 i = substep_index = var("substep_index") 50 51 52 53 54 # method building blocks ------------------------------------------------------

55 -class IntegrateInTime(Record):

56 __slots__ = ["start", "end", "component", "result_name"] 57

58 - def visit(self, processor):

59 processor.integrate_in_time(self)

60

61 -class HistoryUpdate(Record):

62 __slots__ = ["which", "slow_arg", "fast_arg"] 63

64 - def visit(self, processor):

65 processor.history_update(self)

66

67 -class StartSubstepLoop(Record):

68 # everything from this point gets executed substep_count times

69 - def visit(self, processor):

70 processor.start_substep_loop(self)

71

72 -class EndSubstepLoop(Record):

73 __slots__ = ["loop_end"] 74 # everything up to this point gets executed substep_count times 75

76 - def __init__(self, loop_end=n):

77 Record.__init__(self, loop_end=loop_end)

78

79 - def visit(self, processor):

80 processor.end_substep_loop(self)

81 82 83 84 85 # actual method descriptions --------------------------------------------------

86 -class MRABMethod(Record):

87 __slots__ = ["steps", "s2f_hist_is_fast", "result_slow", "result_fast"]

88 89 methods = { 90 "fastest_first_1a": MRABMethod(s2f_hist_is_fast=False, 91 steps=[ 92 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 93 result_name="y_s"), 94 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 95 result_name="y_f"), 96 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 97 which=HIST_F2F), 98 EndSubstepLoop(), 99 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 100 which=HIST_S2F), 101 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 102 which=HIST_F2S), 103 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 104 which=HIST_S2S), 105 ], 106 result_slow="y_s", 107 result_fast="y_f"), 108 "fastest_first_1b": MRABMethod(s2f_hist_is_fast=True, 109 steps=[ 110 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 111 result_name="y_s"), 112 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 113 result_name="y_f"), 114 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 115 which=HIST_F2F), 116 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 117 which=HIST_S2F), 118 EndSubstepLoop(), 119 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 120 which=HIST_F2S), 121 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 122 which=HIST_S2S), 123 ], 124 result_slow="y_s", 125 result_fast="y_f"), 126 "slowest_first_1": MRABMethod(s2f_hist_is_fast=False, 127 steps=[ 128 IntegrateInTime(start=0, end=n, component=CO_SLOW, 129 result_name="\\tilde y_s"), 130 IntegrateInTime(start=0, end=n, component=CO_FAST, 131 result_name="\\tilde y_f"), 132 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 133 which=HIST_S2S), 134 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 135 which=HIST_S2F), 136 StartSubstepLoop(), 137 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 138 result_name="y_s"), 139 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 140 result_name="y_f"), 141 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 142 which=HIST_F2F), 143 EndSubstepLoop(), 144 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 145 which=HIST_F2S), 146 ], 147 result_slow="y_s", 148 result_fast="y_f"), 149 "slowest_first_2a": MRABMethod(s2f_hist_is_fast=False, 150 steps=[ 151 IntegrateInTime(start=0, end=n, component=CO_SLOW, 152 result_name="\\tilde y_s"), 153 IntegrateInTime(start=0, end=n, component=CO_FAST, 154 result_name="\\tilde y_f"), 155 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 156 which=HIST_S2S), 157 StartSubstepLoop(), 158 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 159 result_name="y_s"), 160 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 161 result_name="y_f"), 162 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 163 which=HIST_F2F), 164 EndSubstepLoop(), 165 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 166 which=HIST_F2S), 167 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 168 which=HIST_S2F), 169 ], 170 result_slow="y_s", 171 result_fast="y_f"), 172 "slowest_first_2b": MRABMethod(s2f_hist_is_fast=True, 173 steps=[ 174 IntegrateInTime(start=0, end=n, component=CO_SLOW, 175 result_name="\\tilde y_s"), 176 IntegrateInTime(start=0, end=n, component=CO_FAST, 177 result_name="\\tilde y_f"), 178 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 179 which=HIST_S2S), 180 StartSubstepLoop(), 181 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 182 result_name="y_s"), 183 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 184 result_name="y_f"), 185 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 186 which=HIST_F2F), 187 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 188 which=HIST_S2F), 189 EndSubstepLoop(), 190 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 191 which=HIST_F2S), 192 ], 193 result_slow="y_s", 194 result_fast="y_f"), 195 "slowest_first_3a": MRABMethod(s2f_hist_is_fast=False, 196 steps=[ 197 IntegrateInTime(start=0, end=n, component=CO_SLOW, 198 result_name="\\tilde y_s"), 199 IntegrateInTime(start=0, end=n, component=CO_FAST, 200 result_name="\\tilde y_f"), 201 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 202 which=HIST_S2S), 203 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 204 which=HIST_F2S), 205 StartSubstepLoop(), 206 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 207 result_name="y_s"), 208 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 209 result_name="y_f"), 210 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 211 which=HIST_F2F), 212 EndSubstepLoop(), 213 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 214 which=HIST_S2F), 215 ], 216 result_slow="y_s", 217 result_fast="y_f"), 218 "slowest_first_3b": MRABMethod(s2f_hist_is_fast=True, 219 steps=[ 220 IntegrateInTime(start=0, end=n, component=CO_SLOW, 221 result_name="\\tilde y_s"), 222 IntegrateInTime(start=0, end=n, component=CO_FAST, 223 result_name="\\tilde y_f"), 224 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 225 which=HIST_S2S), 226 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 227 which=HIST_F2S), 228 StartSubstepLoop(), 229 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 230 result_name="y_s"), 231 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 232 result_name="y_f"), 233 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 234 which=HIST_F2F), 235 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 236 which=HIST_S2F), 237 EndSubstepLoop(), 238 ], 239 result_slow="y_s", 240 result_fast="y_f"), 241 "slowest_first_4": MRABMethod(s2f_hist_is_fast=False, 242 steps=[ 243 IntegrateInTime(start=0, end=n, component=CO_SLOW, 244 result_name="\\tilde y_s"), 245 IntegrateInTime(start=0, end=n, component=CO_FAST, 246 result_name="\\tilde y_f"), 247 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 248 which=HIST_S2S), 249 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 250 which=HIST_F2S), 251 HistoryUpdate(slow_arg="\\tilde y_s", fast_arg="\\tilde y_f", 252 which=HIST_S2F), 253 StartSubstepLoop(), 254 IntegrateInTime(start=0, end=i+1, component=CO_SLOW, 255 result_name="y_s"), 256 IntegrateInTime(start=i, end=i+1, component=CO_FAST, 257 result_name="y_f"), 258 HistoryUpdate(slow_arg="y_s", fast_arg="y_f", 259 which=HIST_F2F), 260 EndSubstepLoop(), 261 ], 262 result_slow="y_s", 263 result_fast="y_f") 264 } 265 266 267 268

269 -def _remove_last_yslow_evaluation_from_slowest_first(method):

270 m_steps = method.steps 271 try: 272 loop_start_index = m_steps.index(StartSubstepLoop()) 273 except ValueError: 274 loop_start_index = 0 275 276 loop_end_index = m_steps.index(EndSubstepLoop()) 277 278 before_loop = m_steps[:loop_start_index] 279 loop_body = m_steps[loop_start_index+1:loop_end_index] 280 after_loop = m_steps[loop_end_index+1:] 281 282 new_steps = (before_loop 283 + [StartSubstepLoop()] 284 + loop_body 285 + [EndSubstepLoop(m_steps[loop_end_index].loop_end-1)] 286 ) 287 288 for lb_entry in loop_body+after_loop: 289 if isinstance(lb_entry, IntegrateInTime): 290 if lb_entry.result_name != "y_s": 291 new_steps.append(lb_entry) 292 elif isinstance(lb_entry, HistoryUpdate): 293 assert lb_entry.slow_arg == "y_s" 294 new_steps.append(lb_entry.copy( 295 slow_arg="\\tilde y_s")) 296 else: 297 raise NotImplementedError 298 299 return method.copy(steps=new_steps, result_slow="\\tilde y_s")

300 301 302

303 -def _add_slowest_first_variants(methods):

304 result = {} 305 for name, method in methods.iteritems(): 306 result[name] = method 307 if name.startswith("slowest_first"): 308 result[name+"_no_yslow_reeval"] = \ 309 _remove_last_yslow_evaluation_from_slowest_first(method) 310 311 return result

312 313 314 315 316 methods = _add_slowest_first_variants(methods) 317