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

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 ConstantRetardation, FreundlichSorption, 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 Notes 

56 ----- 

57 The shock has: 

58 - c_left: concentration from faster (upstream) characteristic 

59 - c_right: concentration from slower (downstream) characteristic 

60 - velocity: computed from Rankine-Hugoniot condition 

61 

62 The parent characteristics are deactivated. 

63 

64 Examples 

65 -------- 

66 :: 

67 

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

69 assert shock.satisfies_entropy() 

70 assert not char1.is_active # Parent deactivated 

71 """ 

72 # Get c_min from sorption to determine concentration threshold 

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

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

75 

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

77 # Need to determine if this is: 

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

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

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

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

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

83 # Check velocities to determine who is catching whom 

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

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

86 

87 if vel1 > vel2: 

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

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

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

91 try: 

92 raref = RarefactionWave( 

93 t_start=t_event, 

94 v_start=v_event, 

95 flow=char1.flow, 

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

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

98 sorption=char1.sorption, 

99 ) 

100 except ValueError: 

101 # Rarefaction creation failed - just keep C>0, deactivate C=0 

102 char1.is_active = False 

103 return [] 

104 else: 

105 char1.is_active = False 

106 char2.is_active = False 

107 return [raref] 

108 else: 

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

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

111 char1.is_active = False 

112 return [] 

113 

114 elif char2.concentration <= c_min and char1.concentration > c_min and is_n_lt_1 and c_min == 0: 

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

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

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

118 

119 if vel2 > vel1: 

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

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

122 try: 

123 raref = RarefactionWave( 

124 t_start=t_event, 

125 v_start=v_event, 

126 flow=char1.flow, 

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

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

129 sorption=char1.sorption, 

130 ) 

131 except ValueError: 

132 # Rarefaction creation failed 

133 char2.is_active = False 

134 return [] 

135 else: 

136 char1.is_active = False 

137 char2.is_active = False 

138 return [raref] 

139 else: 

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

141 char2.is_active = False 

142 return [] 

143 

144 # Normal case: analyze velocities to determine wave type 

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

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

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

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

149 

150 if vel1 > vel2: 

151 c_left = char1.concentration 

152 c_right = char2.concentration 

153 else: 

154 c_left = char2.concentration 

155 c_right = char1.concentration 

156 

157 # Create shock at collision point 

158 shock = ShockWave( 

159 t_start=t_event, 

160 v_start=v_event, 

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

162 c_left=c_left, 

163 c_right=c_right, 

164 sorption=char1.sorption, 

165 ) 

166 

167 # Verify entropy condition 

168 if not shock.satisfies_entropy(): 

169 # This shouldn't happen if characteristics collided correctly 

170 msg = ( 

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

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

173 ) 

174 raise RuntimeError(msg) 

175 

176 # Deactivate parent characteristics 

177 char1.is_active = False 

178 char2.is_active = False 

179 

180 return [shock] 

181 

182 

183def handle_shock_collision( 

184 shock1: ShockWave, 

185 shock2: ShockWave, 

186 t_event: float, 

187 v_event: float, 

188) -> list[ShockWave]: 

189 """ 

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

191 

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

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

194 shock. 

195 

196 Parameters 

197 ---------- 

198 shock1 : ShockWave 

199 First shock 

200 shock2 : ShockWave 

201 Second shock 

202 t_event : float 

203 Time of collision [days] 

204 v_event : float 

205 Position of collision [m³] 

206 

207 Returns 

208 ------- 

209 list of ShockWave 

210 Single merged shock wave 

211 

212 Notes 

213 ----- 

214 The merged shock has: 

215 - c_left: from the faster (upstream) shock 

216 - c_right: from the slower (downstream) shock 

217 - velocity: recomputed from Rankine-Hugoniot 

218 

219 The parent shocks are deactivated. 

220 

221 Examples 

222 -------- 

223 :: 

224 

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

226 assert merged.satisfies_entropy() 

227 assert not shock1.is_active # Parents deactivated 

228 """ 

229 # Determine which shock is upstream (faster) 

230 # The shock catching up from behind is upstream 

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

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

233 raise RuntimeError(msg) 

234 if shock1.velocity > shock2.velocity: 

235 c_left = shock1.c_left 

236 c_right = shock2.c_right 

237 else: 

238 c_left = shock2.c_left 

239 c_right = shock1.c_right 

240 

241 # Create merged shock 

242 merged = ShockWave( 

243 t_start=t_event, 

244 v_start=v_event, 

245 flow=shock1.flow, 

246 c_left=c_left, 

247 c_right=c_right, 

248 sorption=shock1.sorption, 

249 ) 

250 

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

252 if not merged.satisfies_entropy(): 

253 # This can happen if the intermediate state causes issues 

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

255 msg = ( 

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

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

258 ) 

259 raise RuntimeError(msg) 

260 

261 # Deactivate parent shocks 

262 shock1.is_active = False 

263 shock2.is_active = False 

264 

265 return [merged] 

266 

267 

268def handle_shock_characteristic_collision( 

269 shock: ShockWave, 

270 char: CharacteristicWave, 

271 t_event: float, 

272 v_event: float, 

273) -> list: 

274 """ 

275 Handle shock catching or being caught by characteristic. 

276 

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

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

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

280 

281 The outcome depends on which wave is faster: 

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

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

284 

285 Parameters 

286 ---------- 

287 shock : ShockWave 

288 Shock wave 

289 char : CharacteristicWave 

290 Characteristic wave 

291 t_event : float 

292 Time of collision [days] 

293 v_event : float 

294 Position of collision [m³] 

295 

296 Returns 

297 ------- 

298 list 

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

300 if expansion, or empty list in edge cases 

301 

302 Notes 

303 ----- 

304 The characteristic concentration modifies one side of the shock: 

305 - If shock catches char: modifies c_right 

306 - If char catches shock: modifies c_left 

307 

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

309 If entropy violated → create rarefaction instead (expansion) 

310 

311 Examples 

312 -------- 

313 :: 

314 

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

316 if new_shock: 

317 assert new_shock[0].satisfies_entropy() 

318 """ 

319 if shock.velocity is None: 

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

321 raise RuntimeError(msg) 

322 shock_vel = shock.velocity 

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

324 

325 if shock_vel > char_vel: 

326 # Shock catching characteristic from behind 

327 # Characteristic is on right side of shock 

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

329 new_shock = ShockWave( 

330 t_start=t_event, 

331 v_start=v_event, 

332 flow=shock.flow, 

333 c_left=shock.c_left, 

334 c_right=char.concentration, 

335 sorption=shock.sorption, 

336 ) 

337 else: 

338 # Characteristic catching shock from behind 

339 # Characteristic is on left side of shock 

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

341 new_shock = ShockWave( 

342 t_start=t_event, 

343 v_start=v_event, 

344 flow=shock.flow, 

345 c_left=char.concentration, 

346 c_right=shock.c_right, 

347 sorption=shock.sorption, 

348 ) 

349 

350 # Check entropy condition 

351 if not new_shock.satisfies_entropy(): 

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

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

354 

355 # Determine head and tail concentrations based on velocity ordering 

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

357 if shock_vel > char_vel: 

358 # Shock was catching characteristic 

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

360 c_head = shock.c_left 

361 c_tail = char.concentration 

362 else: 

363 # Characteristic was catching shock 

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

365 c_head = char.concentration 

366 c_tail = shock.c_right 

367 

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

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

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

371 

372 if head_vel > tail_vel: 

373 # Valid rarefaction - create it 

374 try: 

375 raref = RarefactionWave( 

376 t_start=t_event, 

377 v_start=v_event, 

378 flow=shock.flow, 

379 c_head=c_head, 

380 c_tail=c_tail, 

381 sorption=shock.sorption, 

382 ) 

383 except ValueError: 

384 # Rarefaction validation failed - edge case 

385 # Deactivate waves and return empty 

386 shock.is_active = False 

387 char.is_active = False 

388 return [] 

389 else: 

390 # Deactivate parent waves 

391 shock.is_active = False 

392 char.is_active = False 

393 return [raref] 

394 else: 

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

396 # This is an edge case - deactivate and return empty 

397 shock.is_active = False 

398 char.is_active = False 

399 return [] 

400 

401 # Shock satisfies entropy - return it 

402 # Deactivate parent waves 

403 shock.is_active = False 

404 char.is_active = False 

405 

406 return [new_shock] 

407 

408 

409def handle_shock_rarefaction_collision( 

410 shock: ShockWave, 

411 raref: RarefactionWave, 

412 t_event: float, 

413 v_event: float, 

414 boundary_type: str | None, 

415) -> list: 

416 """ 

417 Handle shock interacting with rarefaction fan with wave splitting. 

418 

419 Implements proper wave splitting for shock-rarefaction interactions, 

420 addressing High Priority #2 from FRONT_TRACKING_REBUILD_PLAN.md. 

421 

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

423 

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

425 creating both a modified rarefaction and a continuing shock 

426 - Be caught by rarefaction head: creates compression wave 

427 

428 Parameters 

429 ---------- 

430 shock : ShockWave 

431 Shock wave 

432 raref : RarefactionWave 

433 Rarefaction wave 

434 t_event : float 

435 Time of collision [days] 

436 v_event : float 

437 Position of collision [m³] 

438 boundary_type : str 

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

440 

441 Returns 

442 ------- 

443 list 

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

445 for tail collision, or compression shock for head collision 

446 

447 Notes 

448 ----- 

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

450 - New shock continuing through rarefaction 

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

452 

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

454 

455 Examples 

456 -------- 

457 :: 

458 

459 waves = handle_shock_rarefaction_collision( 

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

461 ) 

462 """ 

463 if boundary_type == "tail": 

464 # Shock catching rarefaction tail - FULL WAVE SPLITTING 

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

466 

467 # Query rarefaction concentration at collision point 

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

469 raref_c_at_collision = raref.concentration_at_point(v_event, t_event) 

470 

471 if raref_c_at_collision is None: 

472 # Shock not actually inside rarefaction - edge case 

473 # Fall back to simple approach 

474 new_shock = ShockWave( 

475 t_start=t_event, 

476 v_start=v_event, 

477 flow=shock.flow, 

478 c_left=shock.c_left, 

479 c_right=raref.c_tail, 

480 sorption=shock.sorption, 

481 ) 

482 if new_shock.satisfies_entropy(): 

483 raref.is_active = False 

484 shock.is_active = False 

485 return [new_shock] 

486 return [] 

487 

488 # Create shock that continues through rarefaction 

489 # Right state is the rarefaction concentration at collision 

490 new_shock = ShockWave( 

491 t_start=t_event, 

492 v_start=v_event, 

493 flow=shock.flow, 

494 c_left=shock.c_left, 

495 c_right=raref_c_at_collision, 

496 sorption=shock.sorption, 

497 ) 

498 

499 if not new_shock.satisfies_entropy(): 

500 # Complex case - shock doesn't continue 

501 raref.is_active = False 

502 shock.is_active = False 

503 return [] 

504 

505 # Create modified rarefaction with compressed tail 

506 # The portion of rarefaction ahead of shock remains 

507 # New tail is at the collision concentration 

508 c_new_tail = raref_c_at_collision 

509 

510 # Verify head is still faster than new tail 

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

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

513 

514 if head_vel > tail_vel: 

515 # Create modified rarefaction starting from collision point 

516 try: 

517 modified_raref = RarefactionWave( 

518 t_start=t_event, 

519 v_start=v_event, 

520 flow=raref.flow, 

521 c_head=raref.c_head, 

522 c_tail=c_new_tail, 

523 sorption=raref.sorption, 

524 ) 

525 except ValueError: 

526 # Rarefaction validation failed 

527 shock.is_active = False 

528 raref.is_active = False 

529 return [new_shock] 

530 else: 

531 # Deactivate original waves 

532 shock.is_active = False 

533 raref.is_active = False 

534 

535 return [new_shock, modified_raref] 

536 else: 

537 # Rarefaction completely overtaken - only shock continues 

538 shock.is_active = False 

539 raref.is_active = False 

540 return [new_shock] 

541 

542 # boundary_type == 'head' 

543 # Rarefaction head catching shock 

544 # This creates compression between rarefaction head and shock 

545 # May form new compression shock 

546 

547 # Check if compression forms between rarefaction head and shock 

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

549 if shock.velocity is None: 

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

551 raise RuntimeError(msg) 

552 shock_vel = shock.velocity 

553 

554 if raref_head_vel > shock_vel: 

555 # Rarefaction head is faster - creates compression 

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

557 new_shock = ShockWave( 

558 t_start=t_event, 

559 v_start=v_event, 

560 flow=raref.flow, 

561 c_left=raref.c_head, 

562 c_right=shock.c_right, 

563 sorption=raref.sorption, 

564 ) 

565 

566 if new_shock.satisfies_entropy(): 

567 # Compression shock forms 

568 # Deactivate original shock (rarefaction continues) 

569 shock.is_active = False 

570 return [new_shock] 

571 

572 # No compression shock forms - deactivate both for safety 

573 shock.is_active = False 

574 raref.is_active = False 

575 return [] 

576 

577 

578def handle_rarefaction_characteristic_collision( 

579 raref: RarefactionWave, # noqa: ARG001 

580 char: CharacteristicWave, 

581 t_event: float, # noqa: ARG001 

582 v_event: float, # noqa: ARG001 

583 boundary_type: str | None, # noqa: ARG001 

584) -> list: 

585 """ 

586 Handle rarefaction boundary intersecting with characteristic. 

587 

588 **SIMPLIFIED IMPLEMENTATION**: Just deactivates characteristic. See 

589 FRONT_TRACKING_REBUILD_PLAN.md "Known Issues and Future Improvements" 

590 Medium Priority #5. 

591 

592 Parameters 

593 ---------- 

594 raref : RarefactionWave 

595 Rarefaction wave 

596 char : CharacteristicWave 

597 Characteristic wave 

598 t_event : float 

599 Time of collision [days] 

600 v_event : float 

601 Position of collision [m³] 

602 boundary_type : str 

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

604 

605 Returns 

606 ------- 

607 list 

608 List of new waves created (currently always empty) 

609 

610 Notes 

611 ----- 

612 This is a simplified implementation that deactivates the characteristic 

613 without modifying the rarefaction structure. 

614 

615 Current implementation: deactivates characteristic, leaves rarefaction 

616 unchanged. 

617 

618 **Future enhancement**: Should modify rarefaction head/tail concentration 

619 to properly represent the wave structure instead of just absorbing the 

620 characteristic. 

621 """ 

622 # Simplified: characteristic gets absorbed into rarefaction 

623 # More sophisticated: modify rarefaction boundaries 

624 char.is_active = False 

625 return [] 

626 

627 

628def handle_rarefaction_rarefaction_collision( 

629 raref1: RarefactionWave, 

630 raref2: RarefactionWave, 

631 t_event: float, 

632 v_event: float, 

633 boundary_type: str | None, 

634) -> list: 

635 """Handle collision between two rarefaction boundaries. 

636 

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

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

639 topology. Full entropic treatment of rarefaction-rarefaction interactions 

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

641 future enhancement. 

642 

643 Parameters 

644 ---------- 

645 raref1 : RarefactionWave 

646 First rarefaction wave in the collision. 

647 raref2 : RarefactionWave 

648 Second rarefaction wave in the collision. 

649 t_event : float 

650 Time of the boundary intersection [days]. 

651 v_event : float 

652 Position of the intersection [m³]. 

653 boundary_type : str 

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

655 

656 Returns 

657 ------- 

658 list 

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

660 

661 Notes 

662 ----- 

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

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

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

666 computation while deferring complex topology changes. 

667 """ 

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

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

670 return [] 

671 

672 

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

674 """ 

675 Handle wave crossing outlet boundary. 

676 

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

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

679 

680 Parameters 

681 ---------- 

682 wave : Wave 

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

684 t_event : float 

685 Time when wave exits [days] 

686 v_outlet : float 

687 Outlet position [m³] 

688 

689 Returns 

690 ------- 

691 dict 

692 Event record with details about the crossing 

693 

694 Notes 

695 ----- 

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

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

698 

699 The event record includes: 

700 - time: crossing time 

701 - type: 'outlet_crossing' 

702 - wave: reference to wave object 

703 - concentration_left: upstream concentration 

704 - concentration_right: downstream concentration 

705 

706 Examples 

707 -------- 

708 :: 

709 

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

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

712 """ 

713 return { 

714 "time": t_event, 

715 "type": "outlet_crossing", 

716 "wave": wave, 

717 "location": v_outlet, 

718 "concentration_left": wave.concentration_left(), 

719 "concentration_right": wave.concentration_right(), 

720 } 

721 

722 

723def recreate_characteristic_with_new_flow( 

724 char: CharacteristicWave, 

725 t_change: float, 

726 flow_new: float, 

727) -> CharacteristicWave: 

728 """ 

729 Create new characteristic at current position with new flow. 

730 

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

732 velocities. The concentration remains constant, but velocity becomes 

733 flow_new / R(concentration). 

734 

735 Parameters 

736 ---------- 

737 char : CharacteristicWave 

738 Existing characteristic to recreate 

739 t_change : float 

740 Time of flow change [days] 

741 flow_new : float 

742 New flow rate [m³/day] 

743 

744 Returns 

745 ------- 

746 CharacteristicWave 

747 New characteristic at current position with new flow 

748 

749 Notes 

750 ----- 

751 The parent characteristic should be deactivated by the caller. 

752 

753 Examples 

754 -------- 

755 :: 

756 

757 char_old = CharacteristicWave( 

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

759 ) 

760 char_new = recreate_characteristic_with_new_flow( 

761 char_old, t_change=10.0, flow_new=200.0 

762 ) 

763 assert char_new.flow == 200.0 

764 assert char_new.concentration == 5.0 # Concentration unchanged 

765 """ 

766 v_at_change = char.position_at_time(t_change) 

767 

768 if v_at_change is None: 

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

770 raise ValueError(msg) 

771 

772 return CharacteristicWave( 

773 t_start=t_change, 

774 v_start=v_at_change, 

775 flow=flow_new, 

776 concentration=char.concentration, 

777 sorption=char.sorption, 

778 is_active=True, 

779 ) 

780 

781 

782def recreate_shock_with_new_flow( 

783 shock: ShockWave, 

784 t_change: float, 

785 flow_new: float, 

786) -> ShockWave: 

787 """ 

788 Create new shock at current position with new flow. 

789 

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

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

792 

793 Parameters 

794 ---------- 

795 shock : ShockWave 

796 Existing shock to recreate 

797 t_change : float 

798 Time of flow change [days] 

799 flow_new : float 

800 New flow rate [m³/day] 

801 

802 Returns 

803 ------- 

804 ShockWave 

805 New shock at current position with updated velocity 

806 

807 Notes 

808 ----- 

809 The parent shock should be deactivated by the caller. 

810 Shock velocity is automatically recomputed in ShockWave.__post_init__. 

811 

812 Examples 

813 -------- 

814 :: 

815 

816 shock_old = ShockWave( 

817 t_start=0.0, 

818 v_start=0.0, 

819 flow=100.0, 

820 c_left=10.0, 

821 c_right=2.0, 

822 sorption=sorption, 

823 ) 

824 shock_new = recreate_shock_with_new_flow( 

825 shock_old, t_change=10.0, flow_new=200.0 

826 ) 

827 assert shock_new.flow == 200.0 

828 assert ( 

829 shock_new.velocity == 2 * shock_old.velocity 

830 ) # Velocity scales linearly with flow 

831 """ 

832 v_at_change = shock.position_at_time(t_change) 

833 

834 if v_at_change is None: 

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

836 raise ValueError(msg) 

837 

838 return ShockWave( 

839 t_start=t_change, 

840 v_start=v_at_change, 

841 flow=flow_new, 

842 c_left=shock.c_left, 

843 c_right=shock.c_right, 

844 sorption=shock.sorption, 

845 is_active=True, 

846 ) 

847 

848 

849def recreate_rarefaction_with_new_flow( 

850 raref: RarefactionWave, 

851 t_change: float, 

852 flow_new: float, 

853) -> RarefactionWave: 

854 """ 

855 Create new rarefaction at current position with new flow. 

856 

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

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

859 solution pivots at the flow change point. 

860 

861 Parameters 

862 ---------- 

863 raref : RarefactionWave 

864 Existing rarefaction to recreate 

865 t_change : float 

866 Time of flow change [days] 

867 flow_new : float 

868 New flow rate [m³/day] 

869 

870 Returns 

871 ------- 

872 RarefactionWave 

873 New rarefaction at current position with updated velocities 

874 

875 Notes 

876 ----- 

877 The parent rarefaction should be deactivated by the caller. 

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

879 

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

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

882 

883 Examples 

884 -------- 

885 :: 

886 

887 raref_old = RarefactionWave( 

888 t_start=0.0, 

889 v_start=0.0, 

890 flow=100.0, 

891 c_head=10.0, 

892 c_tail=2.0, 

893 sorption=sorption, 

894 ) 

895 raref_new = recreate_rarefaction_with_new_flow( 

896 raref_old, t_change=10.0, flow_new=200.0 

897 ) 

898 assert raref_new.flow == 200.0 

899 assert raref_new.c_head == 10.0 # Concentrations unchanged 

900 """ 

901 v_at_change = raref.position_at_time(t_change) 

902 

903 if v_at_change is None: 

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

905 raise ValueError(msg) 

906 

907 return RarefactionWave( 

908 t_start=t_change, 

909 v_start=v_at_change, 

910 flow=flow_new, 

911 c_head=raref.c_head, 

912 c_tail=raref.c_tail, 

913 sorption=raref.sorption, 

914 is_active=True, 

915 ) 

916 

917 

918def handle_flow_change( 

919 t_change: float, 

920 flow_new: float, 

921 active_waves: list, 

922) -> list: 

923 """ 

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

925 

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

927 positions with updated velocities. This maintains exact analytical computation 

928 while correctly handling time-varying flow. 

929 

930 Parameters 

931 ---------- 

932 t_change : float 

933 Time of flow change [days] 

934 flow_new : float 

935 New flow rate [m³/day] 

936 active_waves : list 

937 All currently active waves 

938 

939 Returns 

940 ------- 

941 list 

942 New waves created at current positions with new flow 

943 

944 Notes 

945 ----- 

946 Parent waves are deactivated by this handler. 

947 

948 Physical interpretation: 

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

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

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

952 

953 Examples 

954 -------- 

955 :: 

956 

957 new_waves = handle_flow_change( 

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

959 ) 

960 assert len(new_waves) == 3 

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

962 """ 

963 new_waves = [] 

964 

965 for wave in active_waves: 

966 if not wave.is_active: 

967 continue 

968 

969 # Create replacement wave with new flow BEFORE deactivating parent 

970 # (position_at_time requires wave to be active) 

971 if isinstance(wave, CharacteristicWave): 

972 new_wave = recreate_characteristic_with_new_flow(wave, t_change, flow_new) 

973 elif isinstance(wave, ShockWave): 

974 new_wave = recreate_shock_with_new_flow(wave, t_change, flow_new) 

975 elif isinstance(wave, RarefactionWave): 

976 new_wave = recreate_rarefaction_with_new_flow(wave, t_change, flow_new) 

977 else: 

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

979 raise TypeError(msg) 

980 

981 new_waves.append(new_wave) 

982 

983 # Deactivate parent wave AFTER recreation 

984 wave.is_active = False 

985 

986 return new_waves 

987 

988 

989def create_inlet_waves_at_time( 

990 c_prev: float, 

991 c_new: float, 

992 t: float, 

993 flow: float, 

994 sorption: FreundlichSorption | ConstantRetardation, 

995 v_inlet: float = 0.0, 

996) -> list: 

997 """ 

998 Create appropriate waves when inlet concentration changes. 

999 

1000 Analyzes the concentration change and creates the physically correct 

1001 wave type based on characteristic velocities. 

1002 

1003 Parameters 

1004 ---------- 

1005 c_prev : float 

1006 Previous concentration [mass/volume] 

1007 c_new : float 

1008 New concentration [mass/volume] 

1009 t : float 

1010 Time of concentration change [days] 

1011 flow : float 

1012 Flow rate [m³/day] 

1013 sorption : FreundlichSorption or ConstantRetardation 

1014 Sorption parameters 

1015 v_inlet : float, optional 

1016 Inlet position [m³], default 0.0 

1017 

1018 Returns 

1019 ------- 

1020 list 

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

1022 

1023 Notes 

1024 ----- 

1025 Wave type logic: 

1026 - vel_new > vel_prev: Compression → create ShockWave 

1027 - vel_new < vel_prev: Expansion → create RarefactionWave 

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

1029 

1030 For shocks, verifies entropy condition before creation. 

1031 

1032 Examples 

1033 -------- 

1034 :: 

1035 

1036 # Step increase from zero creates characteristic 

1037 waves = create_inlet_waves_at_time( 

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

1039 ) 

1040 assert isinstance(waves[0], CharacteristicWave) 

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

1042 waves = create_inlet_waves_at_time( 

1043 c_prev=2.0, c_new=10.0, t=10.0, flow=100.0, sorption=sorption 

1044 ) 

1045 assert isinstance(waves[0], ShockWave) 

1046 """ 

1047 if abs(c_new - c_prev) < EPSILON_CONCENTRATION: # No change 

1048 return [] 

1049 

1050 # Get c_min from sorption if available (determines when to use special treatment) 

1051 c_min = getattr(sorption, "c_min", 0.0) 

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

1053 

1054 # Special case: c_prev ≈ 0 AND this is n<1 with c_min=0 

1055 # For n<1 (lower C travels faster), R(0)=1 is physically correct 

1056 # The C=0 "water" ahead has a well-defined velocity and represents initial condition 

1057 if c_prev <= c_min and is_n_lt_1 and c_min == 0: 

1058 # Create characteristic wave with new concentration 

1059 # The front propagates at v(c_new), leaving c_new behind and 0 ahead 

1060 char = CharacteristicWave( 

1061 t_start=t, 

1062 v_start=v_inlet, 

1063 flow=flow, 

1064 concentration=c_new, 

1065 sorption=sorption, 

1066 ) 

1067 return [char] 

1068 

1069 # Special case: c_new ≈ 0 AND this is n<1 with c_min=0 

1070 # For n<1 (lower C travels faster), clean water (C=0) has well-defined velocity 

1071 if c_new <= c_min and is_n_lt_1 and c_min == 0: 

1072 # Create characteristic wave with zero concentration 

1073 # This represents clean water entering the domain 

1074 char = CharacteristicWave( 

1075 t_start=t, 

1076 v_start=v_inlet, 

1077 flow=flow, 

1078 concentration=c_new, 

1079 sorption=sorption, 

1080 ) 

1081 return [char] 

1082 

1083 # Normal case: analyze velocities to determine wave type 

1084 # For n>1 (higher C travels faster), even stepping down to c_min creates proper waves 

1085 # The velocity analysis will determine if it's a shock, rarefaction, or characteristic 

1086 vel_prev = characteristic_velocity(c_prev, flow, sorption)