2020import matplotlib .pyplot as plt
2121import numpy as np
2222import gotranx
23- import adios4dolfinx
23+ import io4dolfinx
2424import cardiac_geometries
2525import cardiac_geometries .geometry
2626import pulse
3030def mmHg_to_kPa (x ):
3131 return x * 0.133322
3232
33+
3334def custom_json (obj ):
3435 if isinstance (obj , np .float64 ):
3536 return float (obj )
@@ -141,8 +142,6 @@ def run_zenker(outdir: Path):
141142 history ["before_index" ] = before_index
142143 history ["after_index" ] = after_index
143144
144-
145-
146145 Path (zenker_file ).write_text (json .dumps (history , indent = 4 , default = custom_json ))
147146
148147 return json .loads (zenker_file .read_text ())
@@ -202,7 +201,6 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
202201 monitor = mon (ti , states , p )
203202 Tas [i ] = monitor [Ta_index ]
204203
205-
206204 nbeats = 10 if os .environ .get ("CI" ) else 200
207205 times = np .arange (0 , T , dt_cell )
208206 all_times = np .arange (0 , T * nbeats , dt_cell )
@@ -261,7 +259,6 @@ def run_3D_model(
261259 mesh_unit : str = "m" ,
262260 volume2ml : float = 1000.0 ,
263261):
264-
265262 geometry = pulse .HeartGeometry .from_cardiac_geometries (
266263 geo ,
267264 metadata = {"quadrature_degree" : 6 },
@@ -274,7 +271,8 @@ def run_3D_model(
274271 # We use an active stress approach with 30% transverse active stress (see {py:meth}`pulse.active_stress.transversely_active_stress`)
275272
276273 Ta = pulse .Variable (
277- dolfinx .fem .Constant (geometry .mesh , dolfinx .default_scalar_type (0.0 )), "kPa" ,
274+ dolfinx .fem .Constant (geometry .mesh , dolfinx .default_scalar_type (0.0 )),
275+ "kPa" ,
278276 )
279277 active_model = pulse .ActiveStress (geo .f0 , activation = Ta )
280278
@@ -312,11 +310,15 @@ def run_3D_model(
312310 robin = [robin_epi , robin_epi_perp , robin_base ]
313311
314312 lvv_initial = comm .allreduce (geometry .volume ("LV" ), op = MPI .SUM )
315- lv_volume = dolfinx .fem .Constant (geometry .mesh , dolfinx .default_scalar_type (lvv_initial ))
313+ lv_volume = dolfinx .fem .Constant (
314+ geometry .mesh , dolfinx .default_scalar_type (lvv_initial ),
315+ )
316316 lv_cavity = pulse .problem .Cavity (marker = "LV" , volume = lv_volume )
317317
318318 rvv_initial = comm .allreduce (geometry .volume ("RV" ), op = MPI .SUM )
319- rv_volume = dolfinx .fem .Constant (geometry .mesh , dolfinx .default_scalar_type (rvv_initial ))
319+ rv_volume = dolfinx .fem .Constant (
320+ geometry .mesh , dolfinx .default_scalar_type (rvv_initial ),
321+ )
320322 rv_cavity = pulse .problem .Cavity (marker = "RV" , volume = rv_volume )
321323
322324 print ("Initial volumes" , lvv_initial * volume2ml , rvv_initial * volume2ml )
@@ -325,7 +327,11 @@ def run_3D_model(
325327 parameters = {"base_bc" : pulse .problem .BaseBC .free , "mesh_unit" : mesh_unit }
326328 bcs = pulse .BoundaryConditions (robin = robin )
327329 problem = pulse .problem .StaticProblem (
328- model = model , geometry = geometry , bcs = bcs , cavities = cavities , parameters = parameters ,
330+ model = model ,
331+ geometry = geometry ,
332+ bcs = bcs ,
333+ cavities = cavities ,
334+ parameters = parameters ,
329335 )
330336
331337 problem .solve ()
@@ -339,7 +345,7 @@ def run_3D_model(
339345 vtx .write (0.0 )
340346
341347 filename = outdir / Path (f"function_checkpoint_{ label } .bp" )
342- adios4dolfinx .write_mesh (filename , geometry .mesh )
348+ io4dolfinx .write_mesh (filename , geometry .mesh )
343349 output_file = outdir / f"output_{ label } .json"
344350
345351 Ta_history = []
@@ -348,9 +354,9 @@ def callback(model, i: int, t: float, save=True):
348354 Ta_history .append (get_activation (t ))
349355
350356 if save :
351- adios4dolfinx .write_function (filename , problem .u , time = t , name = "displacement" )
357+ io4dolfinx .write_function (filename , problem .u , time = t , name = "displacement" )
352358 vtx .write (t )
353- out = {k : v [:i + 1 ] for k , v in model .history .items ()}
359+ out = {k : v [: i + 1 ] for k , v in model .history .items ()}
354360 out ["Ta" ] = Ta_history
355361 if comm .rank == 0 :
356362 output_file .write_text (json .dumps (out , indent = 4 , default = custom_json ))
@@ -396,7 +402,9 @@ def callback(model, i: int, t: float, save=True):
396402 def p_BiV_func (V_LV , V_RV , t ):
397403 print ("Calculating pressure at time" , t )
398404 value = get_activation (t )
399- print (f"Time{ t } with activation: { value } and volumes: { V_LV } mL (LV) { V_RV } mL (RV)" )
405+ print (
406+ f"Time{ t } with activation: { value } and volumes: { V_LV } mL (LV) { V_RV } mL (RV)" ,
407+ )
400408 old_Ta = Ta .value .value
401409 dTa = value - old_Ta
402410
@@ -422,8 +430,12 @@ def p_BiV_func(V_LV, V_RV, t):
422430 Ta .assign (value )
423431 problem .solve ()
424432
425- lv_pendo_mmHg = circulation .units .kPa_to_mmHg (problem .cavity_pressures [0 ].x .array [0 ] * 1e-3 )
426- rv_pendo_mmHg = circulation .units .kPa_to_mmHg (problem .cavity_pressures [1 ].x .array [0 ] * 1e-3 )
433+ lv_pendo_mmHg = circulation .units .kPa_to_mmHg (
434+ problem .cavity_pressures [0 ].x .array [0 ] * 1e-3 ,
435+ )
436+ rv_pendo_mmHg = circulation .units .kPa_to_mmHg (
437+ problem .cavity_pressures [1 ].x .array [0 ] * 1e-3 ,
438+ )
427439 print (f"Compute pressures: { lv_pendo_mmHg } mmHg (LV) { rv_pendo_mmHg } mmHg (RV)" )
428440 return lv_pendo_mmHg , rv_pendo_mmHg
429441
@@ -432,19 +444,23 @@ def p_BiV_func(V_LV, V_RV, t):
432444 s = circulation .units .ureg ("s" )
433445 add_units = False
434446 lvv_init = (
435- geo .mesh .comm .allreduce (geometry .volume ("LV" , u = problem .u ), op = MPI .SUM ) * 1e6 * 1.0
447+ geo .mesh .comm .allreduce (geometry .volume ("LV" , u = problem .u ), op = MPI .SUM )
448+ * 1e6
449+ * 1.0
436450 )
437451 rvv_init = (
438- geo .mesh .comm .allreduce (geometry .volume ("RV" , u = problem .u ), op = MPI .SUM ) * 1e6 * 1.0
452+ geo .mesh .comm .allreduce (geometry .volume ("RV" , u = problem .u ), op = MPI .SUM )
453+ * 1e6
454+ * 1.0
439455 )
440456 print (f"Initial volume (LV): { lvv_init } mL and (RV): { rvv_init } mL" )
441457 init_state = {
442458 "V_LV" : lvv_initial * 1e6 * mL ,
443459 "V_RV" : rvv_initial * 1e6 * mL ,
444- ' p_AR_PUL' : p_AR_PUL * mmHg ,
445- ' p_AR_SYS' : p_AR_SYS * mmHg ,
446- ' p_VEN_PUL' : p_VEN_PUL * mmHg ,
447- ' p_VEN_SYS' : p_VEN_SYS * mmHg ,
460+ " p_AR_PUL" : p_AR_PUL * mmHg ,
461+ " p_AR_SYS" : p_AR_SYS * mmHg ,
462+ " p_VEN_PUL" : p_VEN_PUL * mmHg ,
463+ " p_VEN_SYS" : p_VEN_SYS * mmHg ,
448464 }
449465
450466 regazzoni_parmeters = circulation .regazzoni2020 .Regazzoni2020 .default_parameters ()
@@ -475,7 +491,9 @@ def p_BiV_func(V_LV, V_RV, t):
475491 )
476492 # Set end time for early stopping if running in CI
477493 end_time = 2 * dt if os .getenv ("CI" ) else None
478- regazzoni .solve (num_beats = num_beats , initial_state = init_state , dt = dt , T = end_time ) #, checkpoint=RR)
494+ regazzoni .solve (
495+ num_beats = num_beats , initial_state = init_state , dt = dt , T = end_time ,
496+ ) # , checkpoint=RR)
479497 regazzoni .print_info ()
480498
481499
@@ -534,10 +552,20 @@ def p_BiV_func(V_LV, V_RV, t):
534552print (f"HR before: { HR_before } , HR after: { HR_after } " )
535553
536554get_activation_before = run_TorOrdLand (
537- comm , outdir , HR_before , Ta_factor = 1 , label = "before" , dt = dt ,
555+ comm ,
556+ outdir ,
557+ HR_before ,
558+ Ta_factor = 1 ,
559+ label = "before" ,
560+ dt = dt ,
538561)
539562get_activation_after = run_TorOrdLand (
540- comm , outdir , HR_after , Ta_factor = C_PRSW_factor , label = "after" , dt = dt ,
563+ comm ,
564+ outdir ,
565+ HR_after ,
566+ Ta_factor = C_PRSW_factor ,
567+ label = "after" ,
568+ dt = dt ,
541569)
542570
543571run_3D_model (
0 commit comments