Gaspard Jankowiak 1 month ago
parent
commit
60bcee6213

+ 1
- 1
configs/test_area_constraint.yaml View File

@@ -5,7 +5,7 @@ params:
5 5
     δt: 1e-3 # time step
6 6
 
7 7
     P: 4*6.15e-1 # pressure
8
-    K: 2*5e-2 # membrane elasticity
8
+    K: 2*5e3 # membrane elasticity
9 9
     Ka: 1 # cortex viscosity
10 10
     mu: 10 # area constraint relaxation constant
11 11
     target_area: 1 # cortex target area

+ 2
- 1
configs_metrics/no_nucleus/varia_mu.yaml View File

@@ -64,7 +64,7 @@ flags:
64 64
     weighted_confinement: true
65 65
     write_animation:      true
66 66
     landscape_plot:       false
67
-    plot_drag:            true
67
+    plot_drag:            false
68 68
     circular_wall:        false
69 69
     cortex:               true
70 70
     centrosome:           true
@@ -73,6 +73,7 @@ flags:
73 73
     debug:                false
74 74
     force_cortex_area:    true
75 75
     write_metrics:        true
76
+    explicit_mt_force:    false
76 77
 
77 78
 config:
78 79
     output_prefix: nonuc_variamu

+ 0
- 98
configs_metrics/nucleus/test_variakb_mu20.yaml View File

@@ -1,98 +0,0 @@
1
-params:
2
-    N: 250 # number of points
3
-    M: 15000 # max. number of iterations
4
-
5
-    δt: 2e-4 # time step
6
-
7
-    P: 4*6.15e-1 # pressure
8
-    K: 3e-1 # membrane elasticity
9
-    Ka: 1 # cortex viscosity
10
-    mu: 10 # area constraint relaxation constant
11
-    target_area: 1.8 # cortex target area
12
-
13
-    c: 1e1 # polymerization speed
14
-
15
-    x0_a: 2*0.25 # initial ellipsis width
16
-    x0_b: 2 # initial ellipsis height
17
-    x0_shift: 0.0 # initial vertical shift
18
-
19
-    # Confinement field
20
-
21
-    f_α: 2e1 # sharpness
22
-    f_β: 0.2 # depth
23
-    f_ω0: 8 # pulsation
24
-    f_σ: 1 # direction
25
-    f_nk: 1 # number of Fourier components
26
-        # to approximate a saw-tooth signal
27
-    f_width: 0.4 # mean width
28
-    f_iwidth: 3.0 # inner width, unused
29
-
30
-    drag_gauss_power: 4.0 # drag_gauss_power
31
-    drag_gauss_width: 3.5 # drag_gauss_width
32
-
33
-    mass_gauss_power: 6.0 # mass_gauss_power
34
-    mass_gauss_width: 1.0 # mass_gauss_width
35
-
36
-    polar_shift: 6.0
37
-
38
-    k_MT: 1e-4              #friction of MT
39
-    MT_potential_power: -2
40
-    MT_factor: 1e-2         #strength of MT
41
-    #
42
-    # Nucleus related parameters
43
-    Nnuc: 200   # number of points on the nucleus
44
-    N_P: -1e0 # pressure
45
-    N_mu: 20 # area constraint relaxation
46
-    N_target_area: 0.7
47
-    N_kb: 1e-3 # bending stiffness
48
-    N_ω:  1e-1 # inplane stiffness
49
-    N_W0: 7e-1  # potential offset
50
-    N_kcont: 5e0 # contact force intensity
51
-    N_αcont: 1e1 # contact force sharpness
52
-    N_kc: 1e-3 # centrosome link stiffness
53
-    N_l0c: 5e-1 # centrosome link rest length
54
-    N_r_init: 2.0e-1 # initial nucleus radius
55
-
56
-flags:
57
-    confine:              true
58
-    adjust_drag:          false
59
-    polymerize:           true
60
-    dryrun:               true
61
-    plot:                 true
62
-    pretty:               false
63
-    continuous:           true
64
-    innerloop:            false
65
-    weighted_confinement: true
66
-    write_animation:      true
67
-    landscape_plot:       false
68
-    plot_drag:            true
69
-    circular_wall:        false
70
-    cortex:               true
71
-    centrosome:           true
72
-    nucleus:              true
73
-    follow_cam:           false
74
-    debug:                false
75
-    force_cortex_area:    true
76
-    write_metrics:        true
77
-
78
-config:
79
-    output_prefix: test_ultimaGaspard7bis
80
-
81
-    plot_period: 20
82
-
83
-    metrics:
84
-        start_iteration: 1
85
-        periods: 2
86
-
87
-    load_state:
88
-        do_load: false
89
-        filename: "initial_conditions/torino_width_1_x.csv"
90
-        do_resample: false
91
-        do_recenter: false
92
-
93
-        init_centro: false
94
-        filename_centro: "initial_conditions/torino_width_1_centro.csv"
95
-
96
-        # shift_nucleus: 0.25 # 3/4 of a period = 3/4 * 2π / f.ω0
97
-
98
-   

+ 0
- 100
configs_metrics/nucleus/test_variamu.yaml View File

@@ -1,100 +0,0 @@
1
-params:
2
-    N: 250 # number of points
3
-    M: 80000 # max. number of iterations
4
-
5
-    δt: 2e-4 # time step
6
-
7
-    P: 4*6.15e-1 # pressure
8
-    K: 3e-1 # membrane elasticity
9
-    Ka: 1 # cortex viscosity
10
-    mu: 10 # area constraint relaxation constant
11
-    target_area: 1.8 # cortex target area
12
-
13
-    c: 1e1 # polymerization speed
14
-
15
-    x0_a: 2*0.25 # initial ellipsis width
16
-    x0_b: 2 # initial ellipsis height
17
-    x0_shift: 0.0 # initial vertical shift
18
-
19
-    # Confinement field
20
-
21
-    f_α: 2e1 # sharpness
22
-    f_β: 0.2 # depth
23
-    f_ω0: 8 # pulsation
24
-    f_σ: 1 # direction
25
-    f_nk: 1 # number of Fourier components
26
-        # to approximate a saw-tooth signal
27
-    f_width: 0.4 # mean width
28
-    f_iwidth: 3.0 # inner width, unused
29
-
30
-    drag_gauss_power: 4.0 # drag_gauss_power
31
-    drag_gauss_width: 3.5 # drag_gauss_width
32
-
33
-    mass_gauss_power: 6.0 # mass_gauss_power
34
-    mass_gauss_width: 1.0 # mass_gauss_width
35
-
36
-    polar_shift: 6.0
37
-
38
-    k_MT: 1e-4              #friction of MT
39
-    MT_potential_power: -2
40
-    MT_factor: 1e-2         #strength of MT
41
-    #
42
-    # Nucleus related parameters
43
-    Nnuc: 200   # number of points on the nucleus
44
-    N_P: -1e0 # pressure
45
-    N_mu: 60 # area constraint relaxation
46
-    N_target_area: 0.7
47
-    N_kb: 5e-3 # bending stiffness
48
-    N_ω:  1e-1 # inplane stiffness
49
-    N_W0: 7e-1  # potential offset
50
-    N_kcont: 5e0 # contact force intensity
51
-    N_αcont: 1e1 # contact force sharpness
52
-    N_kc: 1e-3 # centrosome link stiffness
53
-    N_l0c: 5e-1 # centrosome link rest length
54
-    N_r_init: 2.0e-1 # initial nucleus radius
55
-flags:
56
-    confine:              true
57
-    adjust_drag:          false
58
-    polymerize:           true
59
-    dryrun:               true
60
-    plot:                 true
61
-    pretty:               false
62
-    continuous:           true
63
-    innerloop:            false
64
-    weighted_confinement: true
65
-    write_animation:      true
66
-    landscape_plot:       false
67
-    plot_drag:            true
68
-    circular_wall:        false
69
-    cortex:               true
70
-    centrosome:           true
71
-    nucleus:              true
72
-    follow_cam:           false
73
-    debug:                false
74
-    force_cortex_area:    true
75
-    write_metrics:        true
76
-    explicit_mt_force:    false
77
-
78
-config:
79
-    output_prefix: test_variamu60
80
-
81
-    plot_period: 20
82
-
83
-    recompute_nucleus_each: 5
84
-
85
-    metrics:
86
-        start_iteration: 1
87
-        periods: 2
88
-
89
-    load_state:
90
-        do_load: false
91
-        filename: "initial_conditions/torino_width_1_x.csv"
92
-        do_resample: false
93
-        do_recenter: false
94
-
95
-        init_centro: false
96
-        filename_centro: "initial_conditions/torino_width_1_centro.csv"
97
-
98
-        # shift_nucleus: 0.25 # 3/4 of a period = 3/4 * 2π / f.ω0
99
-
100
-   

+ 86
- 37
src/CellSim.jl View File

@@ -277,6 +277,16 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
277 277
 
278 278
 
279 279
     x = copy(x_init)
280
+    x_candidate = copy(x_init)
281
+
282
+    solver_params = CSC.SolverParams(
283
+            P.δt, # current time step
284
+            1e-3, # maximum time step
285
+            1e-6, # minimum time step
286
+            1.4,  # stepping factor
287
+            1e-4,  # step up error
288
+            1e-2  # step down error
289
+    )
280 290
 
281 291
     Cortex.init_FD_matrices(P)
282 292
     coords, coords_s = Cortex.new_PointCoords(x, P)
@@ -337,6 +347,9 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
337 347
         centro_vr = missing
338 348
     end
339 349
 
350
+    centro_x_candidate = copy(coords.centro_x)
351
+    centro_angle_candidate = copy(coords.centro_angle)
352
+
340 353
     resi, resi_J = Cortex.wrap_residuals(coords, coords_s, potentials, P, F, plotables)
341 354
     if F.innerloop
342 355
         resi_solver = NLsolve.DifferentiableSparseMultivariateFunction(resi, resi_J)
@@ -365,6 +378,7 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
365 378
 
366 379
     k = 0
367 380
     prev_height = 0.0
381
+    current_time = 0
368 382
 
369 383
     metrics = Dict{String, Float64}()
370 384
     post_init_periods = get(config["metrics"], "post_init_periods", 0)
@@ -372,7 +386,6 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
372 386
     post_init_target_max_y = Inf
373 387
     metrics_started = false
374 388
 
375
-
376 389
     stepping = F.debug
377 390
     breakpoint = -1
378 391
     if !stepping
@@ -467,48 +480,84 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
467 480
         end
468 481
 
469 482
 
470
-        try
483
+        # try
471 484
         # inner loop
472 485
         if k > 1
473 486
             Cortex.update_coords(coords, P, x)
474 487
         end
475 488
 
476
-        if F.nucleus
477
-            fill!(potentials.N_W, 0.0)
478
-            fill!(potentials.N_∇W, 0.0)
479
-            fill!(potentials.C_∇W, 0.0)
480
-            fill!(potentials.CS_∇W, 0.0)
481
-            Nucleus.compute_contact_force(potentials, coords, nucleus_coords, P, F)
482
-            if F.centrosome
483
-                Nucleus.compute_centronuclear_force(potentials, coords, nucleus_coords, P, F)
489
+        decreasing_step = false
490
+
491
+        while true
492
+            # adaptive
493
+
494
+            if F.nucleus
495
+                fill!(potentials.N_W, 0.0)
496
+                fill!(potentials.N_∇W, 0.0)
497
+                fill!(potentials.C_∇W, 0.0)
498
+                fill!(potentials.CS_∇W, 0.0)
499
+                Nucleus.compute_contact_force(potentials, coords, nucleus_coords, P, F)
500
+                if F.centrosome
501
+                    Nucleus.compute_centronuclear_force(potentials, coords, nucleus_coords, P, F)
502
+                end
503
+                recompute = (k % config["recompute_nucleus_each"] == 0)
504
+                Nucleus.update_coords(old_nucleus_coords, nucleus_coords, potentials, P, F, temparrays, recompute)
484 505
             end
485
-            recompute = (k % config["recompute_nucleus_each"] == 0)
486
-            Nucleus.update_coords(old_nucleus_coords, nucleus_coords, potentials, P, F, temparrays, recompute)
487
-        end
488 506
 
489
-        if F.cortex
490
-            # cortex evolution
491
-            resi(vec(x), r_x)
492
-            resi_J(vec(x), Jr_x)
493
-        end
507
+            if F.cortex
508
+                # cortex evolution
509
+                resi(vec(x), r_x)
510
+                resi_J(vec(x), Jr_x)
511
+            end
494 512
 
495
-        if F.centrosome
496
-            (centro_A, centro_id_comp, centro_b_ce, centro_b_ce_rhs, centro_b_co_rhs) = Centrosome.assemble_system(P, F, coords, centro_bufs, centro_vr, centro_qw, centro_pc, plotables, potentials)
497
-            # centrosome evolution
498
-            M = ([[Jr_x+centro_id_comp centro_b_ce'];[centro_b_ce centro_A]])
513
+            if F.centrosome
514
+                (centro_A, centro_id_comp, centro_b_ce, centro_b_ce_rhs, centro_b_co_rhs) = Centrosome.assemble_system(P, F, coords, centro_bufs, centro_vr, centro_qw, centro_pc, plotables, potentials)
515
+                # centrosome evolution
516
+                M = ([[Jr_x+centro_id_comp centro_b_ce'];[centro_b_ce centro_A]])
499 517
 
500
-            rhs = [-r_x+P.δt*centro_b_co_rhs; P.δt*centro_b_ce_rhs]
518
+                rhs = [-r_x+P.δt*centro_b_co_rhs; P.δt*centro_b_ce_rhs]
501 519
 
502
-            δx[:] = M\rhs
520
+                δx[:] = M\rhs
503 521
 
504
-            x .+= reshape(δx[1:2P.N], P.N, 2)
522
+                x_candidate .= x .+ reshape(δx[1:2P.N], P.N, 2)
505 523
 
506
-            coords.centro_x .+= δx[(2P.N+1):(2P.N+2)]
507
-            coords.centro_angle .+= δx[2P.N+3]
508
-            Centrosome.compute_vr(P, coords, centro_bufs, centro_vr)
509
-        else
510
-            δx[:] = -Jr_x\r_x
511
-            x .+= + reshape(δx, P.N, 2)
524
+                centro_x_candidate .= coords.centro_x .+ δx[(2P.N+1):(2P.N+2)]
525
+                centro_angle_candidate .= coords.centro_angle .+ δx[2P.N+3]
526
+            else
527
+                δx[:] = -Jr_x\r_x
528
+                x_candidate .= x .+ reshape(δx, P.N, 2)
529
+            end
530
+
531
+            max_displacement = maximum(abs.(x - x_candidate))
532
+            @show max_displacement
533
+
534
+            if max_displacement < solver_params.step_up_error # too small
535
+                if decreasing_step
536
+                    break
537
+                end
538
+                P.δt = min(solver_params.max_δt, P.δt * solver_params.stepping_factor)
539
+                println()
540
+                println("Time step ↑ $(P.δt)")
541
+                println()
542
+                decreasing_step = false
543
+            elseif max_displacement > solver_params.step_down_error # too large
544
+                P.δt = P.δt / solver_params.stepping_factor
545
+                decreasing_step = true
546
+                println()
547
+                println("Time step ↓ $(P.δt)")
548
+                println()
549
+                if P.δt < solver_params.min_δt
550
+                    throw("Time step became too small")
551
+                end
552
+            else # ok
553
+                x .= x_candidate
554
+                if F.centrosome
555
+                    coords.centro_x .= centro_x_candidate
556
+                    coords.centro_angle .= centro_angle_candidate
557
+                    Centrosome.compute_vr(P, coords, centro_bufs, centro_vr)
558
+                end
559
+                break
560
+            end
512 561
         end
513 562
 
514 563
         if metrics_started
@@ -521,12 +570,12 @@ function launch(P::CSC.Params, F::CSC.Flags, config)
521 570
         if F.nucleus
522 571
             Nucleus.copy(old_nucleus_coords, nucleus_coords)
523 572
         end
524
-        catch e
525
-            println()
526
-            println("ERROR at iteration ", k, ":")
527
-            println(e)
528
-            break
529
-        end
573
+        # catch e
574
+            # println()
575
+            # println("ERROR at iteration ", k, ":")
576
+            # println(e)
577
+            # break
578
+        # end
530 579
 
531 580
         # plot
532 581
         if (F.plot && ((k % plot_period == 0) || stepping))

+ 9
- 0
src/CellSimCommon.jl View File

@@ -31,6 +31,15 @@ function init_circ_idx(N::Int64)
31 31
              circshift(1:N, -1), circshift(1:N, -2))
32 32
 end
33 33
 
34
+mutable struct SolverParams
35
+    δt::Float64
36
+    max_δt::Float64
37
+    min_δt::Float64
38
+    stepping_factor::Float64
39
+    step_up_error::Float64
40
+    step_down_error::Float64
41
+end
42
+
34 43
 mutable struct Params
35 44
     M::Int64      # max. number of iterations
36 45
 

Loading…
Cancel
Save