<jupyter_start><jupyter_text>Initial condition<jupyter_code># 2D honeycomb mesh
n_x = 5
n_y = 5
coords = utils.generate_honeycomb_coordinates(n_x, n_y)

# make cell_list for the sheet
sheet = [cl.Cell(i, [x,y], -6.0, True, lambda t: 6 + t) for i, (x, y) in enumerate(coords)]
# delete cells to make it circular
del sheet[24]
del sheet[23]
del sheet[21]
del sheet[20]
del sheet[4]
del sheet[0]
# plot to check what happend
utils.plot_2d_population(sheet)
#prepare consistent initial data
solver_scipy = cbmos.CBModel(ff.PiecewisePolynomial(), scpi.solve_ivp, dim)
t_data_init = [0, 0.0001]
_, initial_sheet = solver_scipy.simulate(sheet, t_data_init, {'muA': 0.21*9.1, 'muR': 9.1, 'rA': rA, 'rR': 1.0/(1.0-np.sqrt(0.21)/3.0), 'n': 1.0, 'p': 1.0}, {}, seed=seed)[-1]


# plot to check what happend
utils.plot_2d_population(initial_sheet)<jupyter_output><empty_output><jupyter_text>Convergence study  Using parameters fitted to relaxation time<jupyter_code>params_cubic = {"mu": 5.70, "s": s, "rA": rA}
muR = 9.1
ratio = 0.21
params_poly = {'muA': ratio*muR, 'muR': muR, 'rA': rA, 'rR': 1.0/(1.0-np.sqrt(ratio)/3.0), 'n': 1.0, 'p': 1.0}
mu_gls=1.95
params_gls = {'mu': mu_gls, 'a':-2*np.log(0.002/mu_gls)}
params = {'cubic': params_cubic, 'pw. quad.': params_poly, 'GLS': params_gls}
dt_ref = 0.0005
N_ref = int(1/dt_ref*tf)+1
t_data_ref = np.arange(0, tf, dt_ref)

ref_sol_dicts = {}
for solver in solver_names:
    print(solver)
    solvers = solver_dicts[solver]
    ref_traj = {}
    for force in force_names: 
        print('>'+force)
        t_data_ref, history = solvers[force].simulate(initial_sheet, t_data_ref, params[force], {'dt': dt_ref}, seed=seed)
        ref_traj[force] = {
            cell.ID: np.array([cell_list[i].position for cell_list in history])
            for i, cell in enumerate(history[0])
        }
    ref_sol_dicts[solver] = ref_traj 
print('Done.')
dt_values = [0.001*1.25**n for n in range(0, 22)]
sol_dicts = {}
sol = {'cubic': [], 'pw. quad.': [], 'GLS': []}
for dt in dt_values:
    N = int(1/dt*tf) + 1
    print([ N, dt])
    #t_data = np.linspace(0,1,N)
    t_data = np.arange(0,tf,dt)
    
    for force in force_names:
        t_data, history = solver_dicts['EF'][force].simulate(initial_sheet, t_data, params[force], {'dt': dt}, seed=seed)
        length = min(len(t_data), len(history))
        
        traj = {
            cell.ID: np.array([cell_list[i].position for cell_list in history])
            for i, cell in enumerate(history[0])
        }
                
        errorx = 0
        errory = 0
        for ID, tr in traj.items():
            interx = np.interp(t_data_ref[:], t_data[:length], tr[:length, 0])
            intery = np.interp(t_data_ref[:], t_data[:length], tr[:length, 1])
                        
            #splinex = CubicSpline(t_data[1:length], tr[:length, 0])
            #spliney = CubicSpline(t_data[1:length], tr[:length, 1])
            #interx = splinex(t_data_ref[1:])
            #intery = spliney(t_data_ref[1:])
            ref = ref_sol_dicts['EF'][force][ID]
            refx = ref[:,0]
            refy = ref[:,1]
            errorx = errorx + np.linalg.norm(interx-refx)/np.linalg.norm(refx)
            errory = errory + np.linalg.norm(intery-refy)/np.linalg.norm(refy)
        error = np.array([errorx, errory])/len(history[0])
        sol[force].append(error)
sol_dicts['EF'] = sol

current_solver = 'midpoint'
sol = {'cubic': [], 'pw. quad.': [], 'GLS': []}
for dt in dt_values:
    N = int(1/dt*tf) + 1
    #print([ N, dt])
    #t_data = np.linspace(0,1,N)
    t_data = np.arange(0,tf,dt)
    
    for force in force_names:
        t_data, history = solver_dicts[current_solver][force].simulate(initial_sheet, t_data, params[force], {'dt': dt}, seed=seed)
        length = min(len(t_data), len(history))
        
        traj = {
            cell.ID: np.array([cell_list[i].position for cell_list in history])
            for i, cell in enumerate(history[0])
        }
        
        errorx = 0
        errory = 0
        for ID, tr in traj.items():
            splinex = CubicSpline(t_data[:length], tr[:length, 0])
            spliney = CubicSpline(t_data[:length], tr[:length, 1])
            interx = splinex(t_data_ref[:])
            intery = spliney(t_data_ref[:])
            ref = ref_sol_dicts[current_solver][force][ID]
            refx = ref[:,0]
            refy = ref[:,1]
            errorx = errorx + np.linalg.norm(interx-refx)/np.linalg.norm(refx)
            errory = errory + np.linalg.norm(intery-refy)/np.linalg.norm(refy)
        error = np.array([errorx, errory])/len(history[0])
        sol[force].append(error)
sol_dicts[current_solver] = sol

current_solver = 'AB'
sol = {'cubic': [], 'pw. quad.': [], 'GLS': []}
for dt in dt_values:
    N = int(1/dt*tf) + 1
    #print([ N, dt])
    #t_data = np.linspace(0,1,N)
    t_data = np.arange(0,tf,dt)
    
    for force in force_names:
        t_data, history = solver_dicts[current_solver][force].simulate(initial_sheet, t_data, params[force], {'dt': dt}, seed=seed)
        length = min(len(t_data), len(history))
        
        traj = {
            cell.ID: np.array([cell_list[i].position for cell_list in history])
            for i, cell in enumerate(history[0])
        }
        
        errorx = 0
        errory = 0
        for ID, tr in traj.items():
            splinex = CubicSpline(t_data[:length], tr[:length, 0])
            spliney = CubicSpline(t_data[:length], tr[:length, 1])
            interx = splinex(t_data_ref[:])
            intery = spliney(t_data_ref[:])
            ref = ref_sol_dicts[current_solver][force][ID]
            refx = ref[:,0]
            refy = ref[:,1]
            errorx = errorx + np.linalg.norm(interx-refx)/np.linalg.norm(refx)
            errory = errory + np.linalg.norm(intery-refy)/np.linalg.norm(refy)
        error = np.array([errorx, errory])/len(history[0])
        sol[force].append(error)
sol_dicts[current_solver] = sol

params = {'legend.fontsize': 'xx-large',
          'figure.figsize': (6.75, 5),
          'lines.linewidth': 3.0,
         'axes.labelsize': 'xx-large',
         'axes.titlesize':'xx-large',
         'xtick.labelsize':'xx-large',
         'ytick.labelsize':'xx-large',
          'legend.fontsize': 'xx-large',
         'font.size': 12,
          'font.family': 'serif',
          "mathtext.fontset": "dejavuserif",
         'axes.titlepad': 12,
        'axes.labelpad': 12}
plt.rcParams.update(params)
# single figure

fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(19.5, 5), sharey=True, gridspec_kw={'wspace': 0.1})

# ax1
for force in force_names: 
    
    ax1.loglog(dt_values, np.sum(np.array(sol_dicts['EF'][force]), axis=1), label=force, color=colors[force], linestyle=linestyles[force])
    #plt.loglog(dt_values, np.array(sol_dicts['EF'][force])[:,1], label=force+' y')
ax1.loglog(dt_values[1:-1], np.array(dt_values[1:-1])*0.2, ':', label='$f(\Delta t)= \Delta t$', color='grey')
ax1.legend(loc='lower right', borderaxespad=0.)
#plt.legend()
ax1.set_title(r'$\bf{(g)}$')
ax1.set(ylabel='$\epsilon_{rel}$')
ax1.set_xlim([7*1e-4, 1.5*1e-1])
ax1.set_ylim([5*1e-10, 5*1e-0])
ax1.set_xticklabels([])
#ax1.set_ylim([5*1e-8, 5*1e-0])
ax1.text(0.0001, 0.00000001, 'Monolayer (38 cells)', fontsize=22, rotation='vertical')

#ax2
for force in force_names: 
    ax2.loglog(dt_values, np.sum(np.array(sol_dicts['midpoint'][force]), axis=1), label=force, color=colors[force], linestyle=linestyles[force])
    #plt.loglog(dt_values, np.array(sol_dicts[current_solver][force])[:,1], label=force+' y')
ax2.loglog(dt_values[1:-1], np.array(dt_values[1:-1])**2*0.4, ':', label='$f(\Delta t)= \Delta t^2$', color='grey')
ax2.legend(borderaxespad=0.)
#ax2.set(xlabel='$\Delta t$')
ax2.set_xlim([7*1e-4, 1.5*1e-1])
ax2.set_title(r'$\bf{(h)}$')
ax2.set_xticklabels([])

#ax2.text(0.0025, 150, 'midpoint method', fontsize=22)


#ax3
for force in force_names: 
    ax3.loglog(dt_values, np.sum(np.array(sol_dicts['AB'][force]), axis=1), label=force, color=colors[force], linestyle=linestyles[force])
    #plt.loglog(dt_values, np.array(sol_dicts[current_solver][force])[:,1], label=force+' y')
ax3.loglog(dt_values[1:-1], np.array(dt_values[1:-1])**2*0.4, ':', label='$f(\Delta t)= \Delta t^2$', color='grey')
ax3.legend(borderaxespad=0.)
#ax3.set(xlabel='$\Delta t$')
ax3.set_xlim([7*1e-4, 1.5*1e-1])
ax3.set_title(r'$\bf{(i)}$')
ax3.set_xticklabels([])

#ax3.text(0.001, 150, 'Adams-Bashforth method', fontsize=22)


plt.savefig('Fig13c_combined.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18453077463140652, 'fraction_numerical': 0.034051018020126375, 'mean_word_length': 3.2416873449131516, 'pattern_counts': {'":': 4, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '400178', 'license': 'MIT', 'n_tokens_mistral': 3557, 'n_tokens_neox': 3416, 'n_words': 604, 'path': 'examples/exp-convergence_study_38cells.ipynb', 'repo_name': 'Ambistic/CellTissue', 'size': '8546'}