Coverage for src / gwtransport / fronttracking / handlers.py: 84%

1""" 

2Event Handlers for Front Tracking. 

3 

4==================================== 

5 

6This module provides handlers for all wave interaction events in the front 

7tracking algorithm. Each handler receives waves involved in an event and 

8returns new waves created by the interaction. 

9 

10All handlers enforce physical correctness: 

11- Mass conservation (Rankine-Hugoniot condition) 

12- Entropy conditions (Lax condition for shocks) 

13- Causality (no backward-traveling information) 

14 

15Handlers modify wave states in-place by deactivating parent waves and 

16creating new child waves. 

17""" 

18 

19from gwtransport.fronttracking.math import FreundlichSorption, SorptionModel, characteristic_velocity 

20from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave 

21 

22# Numerical tolerance constants 

23EPSILON_CONCENTRATION = 1e-15 # Tolerance for checking if concentration change is negligible 

24 

25 

26def handle_characteristic_collision( 

27 char1: CharacteristicWave, 

28 char2: CharacteristicWave, 

29 t_event: float, 

30 v_event: float, 

31) -> list[ShockWave | RarefactionWave]: 

32 """ 

33 Handle collision of two characteristics → create shock or rarefaction. 

34 

35 When two characteristics with different concentrations intersect, they 

36 form a shock discontinuity. The faster characteristic (lower concentration 

37 for n>1) catches the slower one from behind. 

38 

39 Parameters 

40 ---------- 

41 char1 : CharacteristicWave 

42 First characteristic 

43 char2 : CharacteristicWave 

44 Second characteristic 

45 t_event : float 

46 Time of collision [days] 

47 v_event : float 

48 Position of collision [m³] 

49 

50 Returns 

51 ------- 

52 list of ShockWave or RarefactionWave 

53 Single shock or rarefaction wave created at collision point 

54 

55 Raises 

56 ------ 

57 RuntimeError 

58 If the characteristic collision creates a non-entropic shock, 

59 indicating a bug in the collision detection logic. 

60 

61 Notes 

62 ----- 

63 The shock has: 

64 - c_left: concentration from faster (upstream) characteristic 

65 - c_right: concentration from slower (downstream) characteristic 

66 - velocity: computed from Rankine-Hugoniot condition 

67 

68 The parent characteristics are deactivated. 

69 

70 Examples 

71 -------- 

72 :: 

73 

74 shock = handle_characteristic_collision(char1, char2, t=15.0, v=100.0) 

75 assert shock.satisfies_entropy() 

76 assert not char1.is_active # Parent deactivated 

77 """ 

78 # Get c_min from sorption to determine concentration threshold 

79 c_min = getattr(char1.sorption, "c_min", 0.0) 

80 is_n_lt_1 = isinstance(char1.sorption, FreundlichSorption) and char1.sorption.n < 1.0 

81 

82 # Special case: if one characteristic has C near c_min 

83 # Need to determine if this is: 

84 # 1. C≈c_min from initial condition being overtaken by C>0 → keep C>0 

85 # 2. C≈c_min from inlet (clean water) catching C>0 → analyze velocities 

86 # Only use special handling for n<1 with c_min=0 where R(0)=1 

87 if char1.concentration <= c_min and char2.concentration > c_min and is_n_lt_1 and c_min == 0: 

88 # char1 is C≈0, char2 is C>0 

89 # Check velocities to determine who is catching whom 

90 vel1 = characteristic_velocity(char1.concentration, char1.flow, char1.sorption) 

91 vel2 = characteristic_velocity(char2.concentration, char2.flow, char2.sorption) 

92 

93 if vel1 > vel2: 

94 # C=0 is faster → catching C>0 from behind 

95 # For Freundlich n>1: concentration decrease (C>0 → C=0) forms rarefaction 

96 # Physical: clean water (fast) catching contaminated water (slow) 

97 # The rarefaction needs head_velocity > tail_velocity. We just verified 

98 # vel1 > vel2 and vel1 is the head velocity, so the rarefaction is valid. 

99 raref = RarefactionWave( 

100 t_start=t_event, 

101 v_start=v_event, 

102 flow=char1.flow, 

103 c_head=char1.concentration, # C=0 is head (faster) 

104 c_tail=char2.concentration, # C>0 is tail (slower) 

105 sorption=char1.sorption, 

106 ) 

107 char1.is_active = False 

108 char2.is_active = False 

109 return [raref] 

110 # C>0 is faster → C>0 catching C=0 → C=0 is from initial condition 

111 # Just deactivate the C=0 and keep C>0 active 

112 char1.is_active = False 

113 return [] 

114 

115 if char2.concentration <= c_min and char1.concentration > c_min and is_n_lt_1 and c_min == 0: 

116 # char2 is C≈0, char1 is C>0 

117 vel1 = characteristic_velocity(char1.concentration, char1.flow, char1.sorption) 

118 vel2 = characteristic_velocity(char2.concentration, char2.flow, char2.sorption) 

119 

120 if vel2 > vel1: 

121 # C=0 is faster → catching C>0 from behind 

122 # For Freundlich n>1: concentration decrease forms rarefaction 

123 # The rarefaction needs head_velocity > tail_velocity. We just verified 

124 # vel2 > vel1 and vel2 is the head velocity, so the rarefaction is valid. 

125 raref = RarefactionWave( 

126 t_start=t_event, 

127 v_start=v_event, 

128 flow=char1.flow, 

129 c_head=char2.concentration, # C=0 is head (faster) 

130 c_tail=char1.concentration, # C>0 is tail (slower) 

131 sorption=char1.sorption, 

132 ) 

133 char1.is_active = False 

134 char2.is_active = False 

135 return [raref] 

136 # C>0 is faster → C=0 is from initial condition 

137 char2.is_active = False 

138 return [] 

139 

140 # Normal case: analyze velocities to determine wave type 

141 # This now handles all cases for n>1 (higher C travels faster) and 

142 # concentrations above c_min for n<1 (lower C travels faster) 

143 vel1 = characteristic_velocity(char1.concentration, char1.flow, char1.sorption) 

144 vel2 = characteristic_velocity(char2.concentration, char2.flow, char2.sorption) 

145 

146 if vel1 > vel2: 

147 c_left = char1.concentration 

148 c_right = char2.concentration 

149 else: 

150 c_left = char2.concentration 

151 c_right = char1.concentration 

152 

153 # Create shock at collision point 

154 shock = ShockWave( 

155 t_start=t_event, 

156 v_start=v_event, 

157 flow=char1.flow, # Assume same flow (piecewise constant) 

158 c_left=c_left, 

159 c_right=c_right, 

160 sorption=char1.sorption, 

161 ) 

162 

163 # Verify entropy condition 

164 if not shock.satisfies_entropy(): 

165 # This shouldn't happen if characteristics collided correctly 

166 msg = ( 

167 f"Characteristic collision created non-entropic shock at t={t_event:.3f}, V={v_event:.3f}. " 

168 f"c_left={c_left:.3f}, c_right={c_right:.3f}, shock_vel={shock.velocity:.3f}" 

169 ) 

170 raise RuntimeError(msg) 

171 

172 # Deactivate parent characteristics 

173 char1.is_active = False 

174 char2.is_active = False 

175 

176 return [shock] 

177 

178 

179def handle_shock_collision( 

180 shock1: ShockWave, 

181 shock2: ShockWave, 

182 t_event: float, 

183 v_event: float, 

184) -> list[ShockWave]: 

185 """ 

186 Handle collision of two shocks → merge into single shock. 

187 

188 When two shocks collide, they merge into a single shock that connects 

189 the left state of the upstream shock to the right state of the downstream 

190 shock. 

191 

192 Parameters 

193 ---------- 

194 shock1 : ShockWave 

195 First shock 

196 shock2 : ShockWave 

197 Second shock 

198 t_event : float 

199 Time of collision [days] 

200 v_event : float 

201 Position of collision [m³] 

202 

203 Returns 

204 ------- 

205 list of ShockWave 

206 Single merged shock wave 

207 

208 Raises 

209 ------ 

210 RuntimeError 

211 If shock velocities are not set, or if the merged shock violates 

212 the entropy condition. 

213 

214 Notes 

215 ----- 

216 The merged shock has: 

217 - c_left: from the faster (upstream) shock 

218 - c_right: from the slower (downstream) shock 

219 - velocity: recomputed from Rankine-Hugoniot 

220 

221 The parent shocks are deactivated. 

222 

223 Examples 

224 -------- 

225 :: 

226 

227 merged = handle_shock_collision(shock1, shock2, t=20.0, v=200.0) 

228 assert merged.satisfies_entropy() 

229 assert not shock1.is_active # Parents deactivated 

230 """ 

231 # Determine which shock is upstream (faster) 

232 # The shock catching up from behind is upstream 

233 if shock1.velocity is None or shock2.velocity is None: 

234 msg = "Shock velocities should be set in __post_init__" 

235 raise RuntimeError(msg) 

236 if shock1.velocity > shock2.velocity: 

237 c_left = shock1.c_left 

238 c_right = shock2.c_right 

239 else: 

240 c_left = shock2.c_left 

241 c_right = shock1.c_right 

242 

243 # Create merged shock 

244 merged = ShockWave( 

245 t_start=t_event, 

246 v_start=v_event, 

247 flow=shock1.flow, 

248 c_left=c_left, 

249 c_right=c_right, 

250 sorption=shock1.sorption, 

251 ) 

252 

253 # Entropy should be satisfied (both parents were entropic) 

254 if not merged.satisfies_entropy(): 

255 # This can happen if the intermediate state causes issues 

256 # In some cases, the shocks might pass through each other instead 

257 msg = ( 

258 f"Shock merger created non-entropic shock at t={t_event:.3f}. " 

259 f"This may indicate complex wave interaction requiring special handling." 

260 ) 

261 raise RuntimeError(msg) 

262 

263 # Deactivate parent shocks 

264 shock1.is_active = False 

265 shock2.is_active = False 

266 

267 return [merged] 

268 

269 

270def handle_shock_characteristic_collision( 

271 shock: ShockWave, 

272 char: CharacteristicWave, 

273 t_event: float, 

274 v_event: float, 

275) -> list: 

276 """ 

277 Handle shock catching or being caught by characteristic. 

278 

279 When the attempted shock would violate entropy (indicating expansion rather 

280 than compression), a rarefaction wave is created instead to preserve mass 

281 balance. This addresses High Priority #1 from FRONT_TRACKING_REBUILD_PLAN.md. 

282 

283 The outcome depends on which wave is faster: 

284 - If shock is faster: shock catches characteristic, absorbs it 

285 - If characteristic is faster: characteristic catches shock, modifies it 

286 

287 Parameters 

288 ---------- 

289 shock : ShockWave 

290 Shock wave 

291 char : CharacteristicWave 

292 Characteristic wave 

293 t_event : float 

294 Time of collision [days] 

295 v_event : float 

296 Position of collision [m³] 

297 

298 Returns 

299 ------- 

300 list 

301 List containing new wave(s): ShockWave if compression, RarefactionWave 

302 if expansion, or empty list in edge cases 

303 

304 Raises 

305 ------ 

306 RuntimeError 

307 If the shock velocity is not set (i.e., ``shock.velocity`` is None). 

308 

309 Notes 

310 ----- 

311 The characteristic concentration modifies one side of the shock: 

312 - If shock catches char: modifies c_right 

313 - If char catches shock: modifies c_left 

314 

315 If the new shock satisfies entropy → return shock (compression) 

316 If entropy violated → create rarefaction instead (expansion) 

317 

318 Examples 

319 -------- 

320 :: 

321 

322 new_shock = handle_shock_characteristic_collision(shock, char, t=25.0, v=300.0) 

323 if new_shock: 

324 assert new_shock[0].satisfies_entropy() 

325 """ 

326 if shock.velocity is None: 

327 msg = "Shock velocity should be set in __post_init__" 

328 raise RuntimeError(msg) 

329 shock_vel = shock.velocity 

330 char_vel = characteristic_velocity(char.concentration, char.flow, char.sorption) 

331 

332 if shock_vel > char_vel: 

333 # Shock catching characteristic from behind 

334 # Characteristic is on right side of shock 

335 # New shock: c_left unchanged, c_right = char.concentration 

336 new_shock = ShockWave( 

337 t_start=t_event, 

338 v_start=v_event, 

339 flow=shock.flow, 

340 c_left=shock.c_left, 

341 c_right=char.concentration, 

342 sorption=shock.sorption, 

343 ) 

344 else: 

345 # Characteristic catching shock from behind 

346 # Characteristic is on left side of shock 

347 # New shock: c_left = char.concentration, c_right unchanged 

348 new_shock = ShockWave( 

349 t_start=t_event, 

350 v_start=v_event, 

351 flow=shock.flow, 

352 c_left=char.concentration, 

353 c_right=shock.c_right, 

354 sorption=shock.sorption, 

355 ) 

356 

357 # Check entropy condition 

358 if not new_shock.satisfies_entropy(): 

359 # Entropy violated → this is an expansion, not compression 

360 # Create rarefaction wave instead of shock to preserve mass balance 

361 

362 # Determine head and tail concentrations based on velocity ordering 

363 # For a rarefaction: head (faster) follows tail (slower) 

364 if shock_vel > char_vel: 

365 # Shock was catching characteristic 

366 # Expansion between shock.c_left (faster) and char.concentration (slower) 

367 c_head = shock.c_left 

368 c_tail = char.concentration 

369 else: 

370 # Characteristic was catching shock 

371 # Expansion between char.concentration (faster) and shock.c_right (slower) 

372 c_head = char.concentration 

373 c_tail = shock.c_right 

374 

375 # Verify this creates a valid rarefaction (head faster than tail) 

376 head_vel = characteristic_velocity(c_head, shock.flow, shock.sorption) 

377 tail_vel = characteristic_velocity(c_tail, shock.flow, shock.sorption) 

378 

379 if head_vel > tail_vel: 

380 # Valid rarefaction - head_vel > tail_vel guarantees the constructor succeeds. 

381 raref = RarefactionWave( 

382 t_start=t_event, 

383 v_start=v_event, 

384 flow=shock.flow, 

385 c_head=c_head, 

386 c_tail=c_tail, 

387 sorption=shock.sorption, 

388 ) 

389 # Deactivate parent waves 

390 shock.is_active = False 

391 char.is_active = False 

392 return [raref] 

393 # Not a valid rarefaction - waves may pass through each other. 

394 # Edge case (head_vel == tail_vel within machine precision): deactivate and return empty. 

395 shock.is_active = False 

396 char.is_active = False 

397 return [] 

398 

399 # Shock satisfies entropy - return it 

400 # Deactivate parent waves 

401 shock.is_active = False 

402 char.is_active = False 

403 

404 return [new_shock] 

405 

406 

407def handle_shock_rarefaction_collision( 

408 shock: ShockWave, 

409 raref: RarefactionWave, 

410 t_event: float, 

411 v_event: float, 

412 boundary_type: str | None, 

413) -> list: 

414 """ 

415 Handle shock interacting with rarefaction fan with wave splitting. 

416 

417 Implements proper wave splitting for shock-rarefaction interactions, 

418 addressing High Priority #2 from FRONT_TRACKING_REBUILD_PLAN.md. 

419 

420 This is the most complex interaction. A shock can: 

421 

422 - Catch the rarefaction tail: shock penetrates into rarefaction fan, 

423 creating both a modified rarefaction and a continuing shock 

424 - Be caught by rarefaction head: creates compression wave 

425 

426 Parameters 

427 ---------- 

428 shock : ShockWave 

429 Shock wave 

430 raref : RarefactionWave 

431 Rarefaction wave 

432 t_event : float 

433 Time of collision [days] 

434 v_event : float 

435 Position of collision [m³] 

436 boundary_type : str 

437 Which boundary collided: 'head' or 'tail' 

438 

439 Returns 

440 ------- 

441 list 

442 List of new waves created: may include shock and modified rarefaction 

443 for tail collision, or compression shock for head collision 

444 

445 Raises 

446 ------ 

447 RuntimeError 

448 If the shock velocity is not set (i.e., ``shock.velocity`` is None) 

449 when processing a head collision. 

450 

451 Notes 

452 ----- 

453 **Tail collision**: Shock penetrates rarefaction, creating: 

454 - New shock continuing through rarefaction 

455 - Modified rarefaction with compressed tail (if rarefaction not fully overtaken) 

456 

457 **Head collision**: Rarefaction head catches shock, may create compression shock 

458 

459 Examples 

460 -------- 

461 :: 

462 

463 waves = handle_shock_rarefaction_collision( 

464 shock, raref, t=30.0, v=400.0, boundary_type="tail" 

465 ) 

466 """ 

467 if boundary_type == "tail": 

468 # Shock catching rarefaction tail - FULL WAVE SPLITTING 

469 # Shock penetrates into rarefaction fan, need to split waves 

470 

471 # Query rarefaction concentration at collision point 

472 # This tells us where in the rarefaction fan the shock is 

473 raref_c_at_collision = raref.concentration_at_point(v_event, t_event) 

474 

475 if raref_c_at_collision is None: 

476 # Shock not actually inside rarefaction - edge case 

477 # Fall back to simple approach 

478 new_shock = ShockWave( 

479 t_start=t_event, 

480 v_start=v_event, 

481 flow=shock.flow, 

482 c_left=shock.c_left, 

483 c_right=raref.c_tail, 

484 sorption=shock.sorption, 

485 ) 

486 if new_shock.satisfies_entropy(): 

487 raref.is_active = False 

488 shock.is_active = False 

489 return [new_shock] 

490 return [] 

491 

492 # Create shock that continues through rarefaction 

493 # Right state is the rarefaction concentration at collision 

494 new_shock = ShockWave( 

495 t_start=t_event, 

496 v_start=v_event, 

497 flow=shock.flow, 

498 c_left=shock.c_left, 

499 c_right=raref_c_at_collision, 

500 sorption=shock.sorption, 

501 ) 

502 

503 if not new_shock.satisfies_entropy(): 

504 # Complex case - shock doesn't continue 

505 raref.is_active = False 

506 shock.is_active = False 

507 return [] 

508 

509 # Create modified rarefaction with compressed tail 

510 # The portion of rarefaction ahead of shock remains 

511 # New tail is at the collision concentration 

512 c_new_tail = raref_c_at_collision 

513 

514 # Verify head is still faster than new tail 

515 head_vel = characteristic_velocity(raref.c_head, raref.flow, raref.sorption) 

516 tail_vel = characteristic_velocity(c_new_tail, raref.flow, raref.sorption) 

517 

518 if head_vel > tail_vel: 

519 # Create modified rarefaction starting from collision point. 

520 # head_vel > tail_vel guarantees the constructor succeeds. 

521 modified_raref = RarefactionWave( 

522 t_start=t_event, 

523 v_start=v_event, 

524 flow=raref.flow, 

525 c_head=raref.c_head, 

526 c_tail=c_new_tail, 

527 sorption=raref.sorption, 

528 ) 

529 # Deactivate original waves 

530 shock.is_active = False 

531 raref.is_active = False 

532 return [new_shock, modified_raref] 

533 # Rarefaction completely overtaken (head_vel <= tail_vel) - only shock continues 

534 shock.is_active = False 

535 raref.is_active = False 

536 return [new_shock] 

537 

538 # boundary_type == 'head' 

539 # Rarefaction head catching shock 

540 # This creates compression between rarefaction head and shock 

541 # May form new compression shock 

542 

543 # Check if compression forms between rarefaction head and shock 

544 raref_head_vel = characteristic_velocity(raref.c_head, raref.flow, raref.sorption) 

545 if shock.velocity is None: 

546 msg = "Shock velocity should be set in __post_init__" 

547 raise RuntimeError(msg) 

548 shock_vel = shock.velocity 

549 

550 if raref_head_vel > shock_vel: 

551 # Rarefaction head is faster - creates compression 

552 # Try to form shock between rarefaction head and state downstream of original shock 

553 new_shock = ShockWave( 

554 t_start=t_event, 

555 v_start=v_event, 

556 flow=raref.flow, 

557 c_left=raref.c_head, 

558 c_right=shock.c_right, 

559 sorption=raref.sorption, 

560 ) 

561 

562 if new_shock.satisfies_entropy(): 

563 # Compression shock forms 

564 # Deactivate original shock (rarefaction continues) 

565 shock.is_active = False 

566 return [new_shock] 

567 

568 # No compression shock forms - deactivate both for safety 

569 shock.is_active = False 

570 raref.is_active = False 

571 return [] 

572 

573 

574def handle_rarefaction_characteristic_collision( 

575 raref: RarefactionWave, 

576 char: CharacteristicWave, 

577 t_event: float, 

578 v_event: float, 

579 boundary_type: str | None, 

580) -> list: 

581 """ 

582 Handle rarefaction boundary intersecting with characteristic. 

583 

584 Implements the safe option (b) of the front tracking rebuild plan: when a 

585 characteristic intersects either boundary of a rarefaction, the 

586 characteristic is absorbed into the rarefaction provided that the 

587 characteristic's concentration matches the boundary concentration to 

588 within a tight tolerance. If the concentrations differ by more than the 

589 tolerance, deactivating the characteristic would silently destroy mass, 

590 so an informative ``RuntimeError`` is raised instead. 

591 

592 The matching tolerance is based on the rarefaction's own concentration 

593 range; characteristics that are tangent to (and therefore physically 

594 indistinguishable from) the rarefaction boundary pass the check, while 

595 truly distinct characteristics are flagged. 

596 

597 Parameters 

598 ---------- 

599 raref : RarefactionWave 

600 Rarefaction wave whose boundary the characteristic crosses. 

601 char : CharacteristicWave 

602 Characteristic wave that intersects the rarefaction boundary. 

603 t_event : float 

604 Time of collision [days]. 

605 v_event : float 

606 Position of collision [m^3]. 

607 boundary_type : str or None 

608 Which boundary collided: ``'head'`` or ``'tail'``. 

609 

610 Returns 

611 ------- 

612 list 

613 Empty list -- the characteristic is absorbed; no new waves are 

614 created. 

615 

616 Raises 

617 ------ 

618 RuntimeError 

619 If the characteristic's concentration does not match the colliding 

620 rarefaction boundary concentration within tolerance, indicating that 

621 absorption would silently violate mass balance. 

622 

623 Notes 

624 ----- 

625 Future enhancement: properly split the rarefaction at the collision 

626 point and create new wave(s) representing the post-interaction state. 

627 """ 

628 # Tolerance for treating the characteristic as tangent to the rarefaction 

629 # boundary. We use a fraction of the rarefaction's concentration jump, 

630 # falling back to an absolute tolerance when the rarefaction is very 

631 # narrow. 

632 rel_tol = 1e-9 

633 abs_tol = 1e-12 

634 raref_range = abs(raref.c_head - raref.c_tail) 

635 tol = max(rel_tol * raref_range, abs_tol) 

636 

637 if boundary_type == "head": 

638 boundary_c = raref.c_head 

639 elif boundary_type == "tail": 

640 boundary_c = raref.c_tail 

641 else: 

642 msg = f"handle_rarefaction_characteristic_collision: unknown boundary_type {boundary_type!r}" 

643 raise RuntimeError(msg) 

644 

645 if abs(char.concentration - boundary_c) > tol: 

646 msg = ( 

647 f"Rarefaction-characteristic collision at t={t_event:.6f}, V={v_event:.6f} would silently " 

648 f"destroy mass: characteristic concentration {char.concentration:.6g} differs from " 

649 f"rarefaction {boundary_type} concentration {boundary_c:.6g} by " 

650 f"{abs(char.concentration - boundary_c):.3g} (tolerance {tol:.3g}). " 

651 f"Proper wave splitting at the rarefaction boundary is required for this case." 

652 ) 

653 raise RuntimeError(msg) 

654 

655 # Characteristic is tangent to rarefaction boundary -> safe to absorb. 

656 char.is_active = False 

657 return [] 

658 

659 

660def handle_rarefaction_rarefaction_collision( 

661 raref1: RarefactionWave, 

662 raref2: RarefactionWave, 

663 t_event: float, 

664 v_event: float, 

665 boundary_type: str | None, 

666) -> list: 

667 """Handle collision between two rarefaction boundaries. 

668 

669 This handler is intentionally conservative: it records the fact that two 

670 rarefaction fans have intersected but does not yet modify the wave 

671 topology. Full entropic treatment of rarefaction-rarefaction interactions 

672 (potentially involving wave splitting) is reserved for a dedicated 

673 future enhancement. 

674 

675 Parameters 

676 ---------- 

677 raref1 : RarefactionWave 

678 First rarefaction wave in the collision. 

679 raref2 : RarefactionWave 

680 Second rarefaction wave in the collision. 

681 t_event : float 

682 Time of the boundary intersection [days]. 

683 v_event : float 

684 Position of the intersection [m³]. 

685 boundary_type : str 

686 Boundary of the first rarefaction that intersected: 'head' or 'tail'. 

687 

688 Returns 

689 ------- 

690 list 

691 Empty list; no new waves are created at this stage. 

692 

693 Notes 

694 ----- 

695 - Waves remain active so that concentration queries remain valid. 

696 - The FrontTracker records the event in its diagnostics history. 

697 - This is consistent with the design goal of exact analytical 

698 computation while deferring complex topology changes. 

699 """ 

700 # No topology changes yet; keep both rarefactions active. 

701 _ = (raref1, raref2, t_event, v_event, boundary_type) 

702 return [] 

703 

704 

705def handle_outlet_crossing(wave, t_event: float, v_outlet: float) -> dict: 

706 """ 

707 Handle wave crossing outlet boundary. 

708 

709 The wave exits the domain. It remains in the wave list for querying 

710 concentration at earlier times but is marked for different handling. 

711 

712 Parameters 

713 ---------- 

714 wave : Wave 

715 Any wave type (Characteristic, Shock, or Rarefaction) 

716 t_event : float 

717 Time when wave exits [days] 

718 v_outlet : float 

719 Outlet position [m³] 

720 

721 Returns 

722 ------- 

723 dict 

724 Event record with details about the crossing 

725 

726 Notes 

727 ----- 

728 Waves are NOT deactivated when they cross the outlet. They remain active 

729 for concentration queries at points between their origin and outlet. 

730 

731 The event record includes: 

732 - time: crossing time 

733 - type: 'outlet_crossing' 

734 - wave: reference to wave object 

735 - concentration_left: upstream concentration 

736 - concentration_right: downstream concentration 

737 

738 Examples 

739 -------- 

740 :: 

741 

742 event = handle_outlet_crossing(shock, t=50.0, v_outlet=500.0) 

743 print(f"Wave exited at t={event['time']:.2f}") 

744 """ 

745 return { 

746 "time": t_event, 

747 "type": "outlet_crossing", 

748 "wave": wave, 

749 "location": v_outlet, 

750 "concentration_left": wave.concentration_left(), 

751 "concentration_right": wave.concentration_right(), 

752 } 

753 

754 

755def recreate_characteristic_with_new_flow( 

756 char: CharacteristicWave, 

757 t_change: float, 

758 flow_new: float, 

759) -> CharacteristicWave: 

760 """ 

761 Create new characteristic at current position with new flow. 

762 

763 When flow changes, existing characteristics must be recreated with updated 

764 velocities. The concentration remains constant, but velocity becomes 

765 flow_new / R(concentration). 

766 

767 Parameters 

768 ---------- 

769 char : CharacteristicWave 

770 Existing characteristic to recreate 

771 t_change : float 

772 Time of flow change [days] 

773 flow_new : float 

774 New flow rate [m³/day] 

775 

776 Returns 

777 ------- 

778 CharacteristicWave 

779 New characteristic at current position with new flow 

780 

781 Raises 

782 ------ 

783 ValueError 

784 If the characteristic is not yet active at ``t_change``. 

785 

786 Notes 

787 ----- 

788 The parent characteristic should be deactivated by the caller. 

789 

790 Examples 

791 -------- 

792 :: 

793 

794 char_old = CharacteristicWave( 

795 t_start=0.0, v_start=0.0, flow=100.0, concentration=5.0, sorption=sorption 

796 ) 

797 char_new = recreate_characteristic_with_new_flow( 

798 char_old, t_change=10.0, flow_new=200.0 

799 ) 

800 assert char_new.flow == 200.0 

801 assert char_new.concentration == 5.0 # Concentration unchanged 

802 """ 

803 v_at_change = char.position_at_time(t_change) 

804 

805 if v_at_change is None: 

806 msg = f"Characteristic not yet active at t={t_change}" 

807 raise ValueError(msg) 

808 

809 return CharacteristicWave( 

810 t_start=t_change, 

811 v_start=v_at_change, 

812 flow=flow_new, 

813 concentration=char.concentration, 

814 sorption=char.sorption, 

815 is_active=True, 

816 ) 

817 

818 

819def recreate_shock_with_new_flow( 

820 shock: ShockWave, 

821 t_change: float, 

822 flow_new: float, 

823) -> ShockWave: 

824 """ 

825 Create new shock at current position with new flow. 

826 

827 When flow changes, shock velocity must be recomputed using Rankine-Hugoniot 

828 condition with the new flow: s = flow*(c_R - c_L) / (C_total(c_R) - C_total(c_L)). 

829 

830 Parameters 

831 ---------- 

832 shock : ShockWave 

833 Existing shock to recreate 

834 t_change : float 

835 Time of flow change [days] 

836 flow_new : float 

837 New flow rate [m³/day] 

838 

839 Returns 

840 ------- 

841 ShockWave 

842 New shock at current position with updated velocity 

843 

844 Raises 

845 ------ 

846 ValueError 

847 If the shock is not yet active at ``t_change``. 

848 

849 Notes 

850 ----- 

851 The parent shock should be deactivated by the caller. 

852 Shock velocity is automatically recomputed in ShockWave.__post_init__. 

853 

854 Examples 

855 -------- 

856 :: 

857 

858 shock_old = ShockWave( 

859 t_start=0.0, 

860 v_start=0.0, 

861 flow=100.0, 

862 c_left=10.0, 

863 c_right=2.0, 

864 sorption=sorption, 

865 ) 

866 shock_new = recreate_shock_with_new_flow( 

867 shock_old, t_change=10.0, flow_new=200.0 

868 ) 

869 assert shock_new.flow == 200.0 

870 assert ( 

871 shock_new.velocity == 2 * shock_old.velocity 

872 ) # Velocity scales linearly with flow 

873 """ 

874 v_at_change = shock.position_at_time(t_change) 

875 

876 if v_at_change is None: 

877 msg = f"Shock not yet active at t={t_change}" 

878 raise ValueError(msg) 

879 

880 return ShockWave( 

881 t_start=t_change, 

882 v_start=v_at_change, 

883 flow=flow_new, 

884 c_left=shock.c_left, 

885 c_right=shock.c_right, 

886 sorption=shock.sorption, 

887 is_active=True, 

888 ) 

889 

890 

891def recreate_rarefaction_with_new_flow( 

892 raref: RarefactionWave, 

893 t_change: float, 

894 flow_new: float, 

895) -> RarefactionWave: 

896 """ 

897 Create new rarefaction at current position with new flow. 

898 

899 When flow changes, rarefaction head and tail velocities are updated. 

900 The fan structure (c_head, c_tail) is preserved, but the self-similar 

901 solution pivots at the flow change point. 

902 

903 Parameters 

904 ---------- 

905 raref : RarefactionWave 

906 Existing rarefaction to recreate 

907 t_change : float 

908 Time of flow change [days] 

909 flow_new : float 

910 New flow rate [m³/day] 

911 

912 Returns 

913 ------- 

914 RarefactionWave 

915 New rarefaction at current position with updated velocities 

916 

917 Raises 

918 ------ 

919 ValueError 

920 If the rarefaction is not yet active at ``t_change``. 

921 

922 Notes 

923 ----- 

924 The parent rarefaction should be deactivated by the caller. 

925 The rarefaction fan "pivots" at (v_at_change, t_change). 

926 

927 Before: R(C) = flow_old * (t - t_start_old) / (v - v_start_old) 

928 After: R(C) = flow_new * (t - t_change) / (v - v_at_change) 

929 

930 Examples 

931 -------- 

932 :: 

933 

934 raref_old = RarefactionWave( 

935 t_start=0.0, 

936 v_start=0.0, 

937 flow=100.0, 

938 c_head=10.0, 

939 c_tail=2.0, 

940 sorption=sorption, 

941 ) 

942 raref_new = recreate_rarefaction_with_new_flow( 

943 raref_old, t_change=10.0, flow_new=200.0 

944 ) 

945 assert raref_new.flow == 200.0 

946 assert raref_new.c_head == 10.0 # Concentrations unchanged 

947 """ 

948 v_at_change = raref.position_at_time(t_change) 

949 

950 if v_at_change is None: 

951 msg = f"Rarefaction not yet active at t={t_change}" 

952 raise ValueError(msg) 

953 

954 return RarefactionWave( 

955 t_start=t_change, 

956 v_start=v_at_change, 

957 flow=flow_new, 

958 c_head=raref.c_head, 

959 c_tail=raref.c_tail, 

960 sorption=raref.sorption, 

961 is_active=True, 

962 ) 

963 

964 

965def handle_flow_change( 

966 t_change: float, 

967 flow_new: float, 

968 active_waves: list, 

969) -> list: 

970 """ 

971 Handle flow change event by recreating all active waves with new flow. 

972 

973 When flow changes, all existing waves must be recreated at their current 

974 positions with updated velocities. This maintains exact analytical computation 

975 while correctly handling time-varying flow. 

976 

977 Parameters 

978 ---------- 

979 t_change : float 

980 Time of flow change [days] 

981 flow_new : float 

982 New flow rate [m³/day] 

983 active_waves : list 

984 All currently active waves 

985 

986 Returns 

987 ------- 

988 list 

989 New waves created at current positions with new flow 

990 

991 Raises 

992 ------ 

993 TypeError 

994 If an unrecognized wave type is encountered in ``active_waves``. 

995 

996 Notes 

997 ----- 

998 Parent waves are deactivated by this handler. 

999 

1000 Physical interpretation: 

1001 - Characteristics: velocity changes from flow_old/R(c) to flow_new/R(c) 

1002 - Shocks: Rankine-Hugoniot velocity recomputed with new flow 

1003 - Rarefactions: fan pivots at (v_change, t_change) 

1004 

1005 Examples 

1006 -------- 

1007 :: 

1008 

1009 new_waves = handle_flow_change( 

1010 t_change=10.0, flow_new=200.0, active_waves=[char1, shock1, raref1] 

1011 ) 

1012 assert len(new_waves) == 3 

1013 assert all(w.flow == 200.0 for w in new_waves) 

1014 """ 

1015 new_waves = [] 

1016 

1017 for wave in active_waves: 

1018 if not wave.is_active: 

1019 continue 

1020 

1021 # Create replacement wave with new flow BEFORE deactivating parent 

1022 # (position_at_time requires wave to be active) 

1023 if isinstance(wave, CharacteristicWave): 

1024 new_wave = recreate_characteristic_with_new_flow(wave, t_change, flow_new) 

1025 elif isinstance(wave, ShockWave): 

1026 new_wave = recreate_shock_with_new_flow(wave, t_change, flow_new) 

1027 elif isinstance(wave, RarefactionWave): 

1028 new_wave = recreate_rarefaction_with_new_flow(wave, t_change, flow_new) 

1029 else: 

1030 msg = f"Unknown wave type: {type(wave)}" 

1031 raise TypeError(msg) 

1032 

1033 new_waves.append(new_wave) 

1034 

1035 # Deactivate parent wave AFTER recreation 

1036 wave.is_active = False 

1037 

1038 return new_waves 

1039 

1040 

1041def create_inlet_waves_at_time( 

1042 c_prev: float, 

1043 c_new: float, 

1044 t: float, 

1045 flow: float, 

1046 sorption: SorptionModel, 

1047 v_inlet: float = 0.0, 

1048) -> list: 

1049 """ 

1050 Create appropriate waves when inlet concentration changes. 

1051 

1052 Analyzes the concentration change and creates the physically correct 

1053 wave type based on characteristic velocities. 

1054 

1055 Parameters 

1056 ---------- 

1057 c_prev : float 

1058 Previous concentration [mass/volume] 

1059 c_new : float 

1060 New concentration [mass/volume] 

1061 t : float 

1062 Time of concentration change [days] 

1063 flow : float 

1064 Flow rate [m³/day] 

1065 sorption : FreundlichSorption or ConstantRetardation 

1066 Sorption parameters 

1067 v_inlet : float, optional 

1068 Inlet position [m³], default 0.0 

1069 

1070 Returns 

1071 ------- 

1072 list 

1073 List of newly created waves (typically 1 wave per concentration change) 

1074 

1075 Notes 

1076 ----- 

1077 Wave type logic: 

1078 - vel_new > vel_prev: Compression → create ShockWave 

1079 - vel_new < vel_prev: Expansion → create RarefactionWave 

1080 - vel_new == vel_prev: Contact discontinuity → create CharacteristicWave 

1081 

1082 For shocks, verifies entropy condition before creation. 

1083 

1084 Examples 

1085 -------- 

1086 :: 

1087 

1088 # Step increase from zero creates characteristic 

1089 waves = create_inlet_waves_at_time( 

1090 c_prev=0.0, c_new=10.0, t=10.0, flow=100.0, sorption=sorption 

1091 ) 

1092 assert isinstance(waves[0], CharacteristicWave) 

1093 # Step between nonzero values creates shock for n>1 (compression)