convergence test for SharpClaw 3D using acoustics_3d_heterogeneous from PyClaw examples.

convergence_test.py

def verify_classic_sharp(claw1,claw2,p,base_name):
    import os
    import numpy as np    

    pfinal1  = claw1.frames[claw1.num_output_times].state.get_q_global()       
    pfinal2  = claw2.frames[claw2.num_output_times].state.get_q_global()
    grid = claw1.solution.state.grid

    plot_claws(grid.x.centers,pfinal1[0,:,0,0],pfinal2[0,:,0,0],base_name+'_x_'+str(2**p))
    plot_claws(grid.y.centers,pfinal1[0,0,:,0],pfinal2[0,0,:,0],base_name+'_y_'+str(2**p))
    plot_claws(grid.z.centers,pfinal1[0,0,0,:],pfinal2[0,0,0,:],base_name+'_z_'+str(2**p))

    pfinal1 = pfinal1[0,:,:,:].reshape(-1)
    pfinal2 = pfinal2[0,:,:,:].reshape(-1)

    final_difference = np.prod(grid.delta)*np.linalg.norm(pfinal1-pfinal2,ord=1)
    np.save(base_name+'_diff_res'+str(p),final_difference)
    print 'difference at same order: ',final_difference

    return final_difference

def plot_claws(grid,claw1,claw2,fig_name):
    from matplotlib import pylab as plt

    plt.close('all')
    plt.figure()
    plt.plot(grid,claw1,'b-',label='sharp')
    plt.plot(grid,claw2,'r-',label='classic')
    plt.title(fig_name)
    plt.legend(frameon=False)
    plt.savefig('./'+fig_name+'.png',dpi=240)
    plt.close() 


def verify_sharp_vs_classic(claw1,claw2,base_name):
    import os
    import numpy as np

    pfinal1  = claw1.frames[claw1.num_output_times].state.get_q_global()       
    pfinal2  = claw2.frames[claw2.num_output_times].state.get_q_global()
    grid = claw1.solution.state.grid    

    pfinal2 = pfinal2[0,::2,::2,::2] + pfinal2[0,1::2,::2,::2] + pfinal2[0,::2,::2,::2] + pfinal2[0,1::2,::2,::2] + \
            pfinal2[0,::2,::2,1::2] + pfinal2[0,1::2,::2,1::2] + pfinal2[0,::2,::2,1::2] + pfinal2[0,1::2,::2,1::2]

    pfinal2 = pfinal2/8.0

    pfinal2 = pfinal2[:,:,:].reshape(-1)
    pfinal1 = pfinal1[0,:,:,:].reshape(-1)

    final_difference = np.prod(grid.delta)*np.linalg.norm(pfinal1-pfinal2,ord=1)
    np.save(base_name+'_diff_res_vs_classic'+str(p),final_difference)
    print 'sharp vs classic: ',final_difference
    return final_difference

def verify_self_convergence(claw_old,claw_new,base_name):
    import os
    import numpy as np

    pfinal1  = claw_old.frames[claw_old.num_output_times].state.get_q_global()
    pfinal2  = claw_new.frames[claw_new.num_output_times].state.get_q_global()       
    grid = claw_new.solution.state.grid

    pfinal2 = pfinal2[0,::2,::2,::2] + pfinal2[0,1::2,::2,::2] + pfinal2[0,::2,::2,::2] + pfinal2[0,1::2,::2,::2] + \
            pfinal2[0,::2,::2,1::2] + pfinal2[0,1::2,::2,1::2] + pfinal2[0,::2,::2,1::2] + pfinal2[0,1::2,::2,1::2]

    pfinal2 = pfinal2/8.0

    pfinal2 = pfinal2[:,:,:].reshape(-1)
    pfinal1 = pfinal1[0,:,:,:].reshape(-1)

    self_difference = np.prod(grid.delta)*np.linalg.norm(pfinal1-pfinal2,ord=1)
    np.save(base_name+'_diff_res_self'+str(p),self_difference)
    print 'self difference: ',self_difference
    return self_difference


import acoustics_3d_interface as acoustics
import numpy as np
import matplotlib
matplotlib.use('Agg')

base_name = 'sol2diff'
test_results = np.empty([2,3,7])
test_sharp   = acoustics.setup(solver_type='sharpclaw',test='heterogeneous',disable_output=False,hmx=8,hmy=8,hmz=8,use_petsc=0)
test_classic = acoustics.setup(solver_type='classic',test='heterogeneous',disable_output=False,hmx=16,hmy=16,hmz=16,use_petsc=0)
test_sharp.run()
test_classic.run()

for p in xrange(4, 8, 1):

    test_sharp_old = test_sharp
    test_classic_old = test_classic
    ns = 2**p
    nc = 2**(p+1)
    test_sharp = acoustics.setup(solver_type='sharpclaw',test='heterogeneous',disable_output=False,hmx=ns,hmy=ns,hmz=ns,use_petsc=0)
    test_sharp.run()
    test_classic = acoustics.setup(solver_type='classic',test='heterogeneous',disable_output=False,hmx=nc,hmy=nc,hmz=nc,use_petsc=0)
    test_classic.run()

    test_results[:,0,p-5] = ns
    test_results[0,1,p-5] = verify_classic_sharp(test_sharp,test_classic_old,p,base_name)
    test_results[1,1,p-5] = verify_sharp_vs_classic(test_sharp,test_classic,base_name)
    test_results[1,2,p-5] = verify_self_convergence(test_sharp_old,test_sharp,base_name)



np.save(base_name,test_results)