Coverage for pySDC/projects/AllenCahn_Bayreuth/run_temp_forcing_verification.py: 99%
109 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +0000
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +0000
1from argparse import ArgumentParser
2import json
3import numpy as np
4from mpi4py import MPI
5from mpi4py_fft import newDistArray
7from pySDC.helpers.stats_helper import get_sorted
9from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
10from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
11from pySDC.implementations.problem_classes.AllenCahn_Temp_MPIFFT import allencahn_temp_imex
12from pySDC.implementations.transfer_classes.TransferMesh_MPIFFT import fft_to_fft
15def run_simulation(name='', spectral=None, nprocs_time=None, nprocs_space=None, dt=None, cwd='.'):
16 """
17 A test program to do PFASST runs for the AC equation with temperature-based forcing
19 (slightly inefficient, but will run for a few seconds only)
21 Args:
22 name (str): name of the run, will be used to distinguish different setups
23 spectral (bool): run in real or spectral space
24 nprocs_time (int): number of processors in time
25 nprocs_space (int): number of processors in space (None if serial)
26 dt (float): time-step size
27 cwd (str): current working directory
28 """
30 # set MPI communicator
31 comm = MPI.COMM_WORLD
33 world_rank = comm.Get_rank()
34 world_size = comm.Get_size()
36 # split world communicator to create space-communicators
37 if nprocs_space is not None:
38 color = int(world_rank / nprocs_space)
39 else:
40 color = int(world_rank / 1)
41 space_comm = comm.Split(color=color)
42 space_rank = space_comm.Get_rank()
43 space_size = space_comm.Get_size()
45 assert world_size == space_size, 'This script cannot run parallel-in-time with MPI, only spatial parallelism'
47 # initialize level parameters
48 level_params = dict()
49 level_params['restol'] = 1e-12
50 level_params['dt'] = dt
51 level_params['nsweeps'] = [1]
53 # initialize sweeper parameters
54 sweeper_params = dict()
55 sweeper_params['quad_type'] = 'RADAU-RIGHT'
56 sweeper_params['num_nodes'] = [3]
57 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
58 sweeper_params['initial_guess'] = 'spread'
60 # initialize problem parameters
61 problem_params = dict()
62 problem_params['L'] = 1.0
63 problem_params['nvars'] = [(128, 128), (32, 32)]
64 problem_params['eps'] = [0.04]
65 problem_params['radius'] = 0.25
66 problem_params['TM'] = 1.0
67 problem_params['D'] = 0.1
68 problem_params['dw'] = [21.0]
69 problem_params['comm'] = space_comm
70 problem_params['init_type'] = 'circle'
71 problem_params['spectral'] = spectral
73 # initialize step parameters
74 step_params = dict()
75 step_params['maxiter'] = 50
77 # initialize controller parameters
78 controller_params = dict()
79 controller_params['logger_level'] = 30 if space_rank == 0 else 99 # set level depending on rank
80 controller_params['predict_type'] = 'pfasst_burnin'
82 # fill description dictionary for easy step instantiation
83 description = dict()
84 description['problem_params'] = problem_params # pass problem parameters
85 description['sweeper_class'] = imex_1st_order
86 description['sweeper_params'] = sweeper_params # pass sweeper parameters
87 description['level_params'] = level_params # pass level parameters
88 description['step_params'] = step_params # pass step parameters
89 description['space_transfer_class'] = fft_to_fft
90 description['problem_class'] = allencahn_temp_imex
92 # set time parameters
93 t0 = 0.0
94 Tend = 1 * 0.001
96 if space_rank == 0:
97 out = f'---------> Running {name} with spectral={spectral} and {space_size} process(es) in space...'
98 print(out)
100 # instantiate controller
101 controller = controller_nonMPI(num_procs=nprocs_time, controller_params=controller_params, description=description)
103 # get initial values on finest level
104 P = controller.MS[0].levels[0].prob
105 uinit = P.u_exact(t0)
107 # call main function to get things done...
108 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
110 if space_rank == 0:
111 # convert filtered statistics of iterations count, sorted by time
112 iter_counts = get_sorted(stats, type='niter', sortby='time')
113 niters = np.mean(np.array([item[1] for item in iter_counts]))
114 out = f'Mean number of iterations: {niters:.4f}'
115 print(out)
117 # get setup time
118 timing = get_sorted(stats, type='timing_setup', sortby='time')
119 out = f'Setup time: {timing[0][1]:.4f} sec.'
120 print(out)
122 # get running time
123 timing = get_sorted(stats, type='timing_run', sortby='time')
124 out = f'Time to solution: {timing[0][1]:.4f} sec.'
125 print(out)
127 refname = f'{cwd}/data/AC-reference-tempforce_00001000'
128 with open(f'{refname}.json', 'r') as fp:
129 obj = json.load(fp)
131 array = np.fromfile(f'{refname}.dat', dtype=obj['datatype'])
132 array = array.reshape(obj['shape'], order='C')
134 if spectral:
135 ureal = newDistArray(P.fft, False)
136 ureal = P.fft.backward(uend[..., 0], ureal)
137 Treal = newDistArray(P.fft, False)
138 Treal = P.fft.backward(uend[..., 1], Treal)
139 err = max(np.amax(abs(ureal - array[..., 0])), np.amax(abs(Treal - array[..., 1])))
140 else:
141 err = abs(array - uend)
143 out = '...Done <---------\n'
144 print(out)
146 return err
149def main(nprocs_space=None, cwd='.'):
150 """
151 Little helper routine to run the whole thing
153 Args:
154 nprocs_space (int): number of processors in space (None if serial)
155 cwd (str): current working directory
156 """
157 name = 'AC-test-tempforce'
159 nsteps = [2**i for i in range(4)]
161 errors = [1]
162 orders = []
163 for n in nsteps:
164 err = run_simulation(name=name, spectral=False, nprocs_time=n, nprocs_space=nprocs_space, dt=1e-03 / n, cwd=cwd)
165 errors.append(err)
166 orders.append(np.log(errors[-1] / errors[-2]) / np.log(0.5))
167 print(f'Error: {errors[-1]:6.4e}')
168 print(f'Order of accuracy: {orders[-1]:4.2f}\n')
170 assert errors[2 + 1] < 8e-10, f'Errors are too high, got {errors[2 + 1]}'
171 assert np.isclose(orders[3], 5.3, rtol=2e-02), f'Order of accuracy is not within tolerance, got {orders[3]}'
173 print()
175 errors = [1]
176 orders = []
177 for n in nsteps:
178 err = run_simulation(name=name, spectral=True, nprocs_time=n, nprocs_space=nprocs_space, dt=1e-03 / n, cwd=cwd)
179 errors.append(err)
180 orders.append(np.log(errors[-1] / errors[-2]) / np.log(0.5))
181 print(f'Error: {errors[-1]:6.4e}')
182 print(f'Order of accuracy: {orders[-1]:4.2f}\n')
184 assert errors[2 + 1] < 8e-10, f'Errors are too high, got {errors[2 + 1]}'
185 assert np.isclose(
186 orders[1].view(np.ndarray), 4.6, rtol=7e-02
187 ), f'Order of accuracy is not within tolerance, got {orders[1]}'
190if __name__ == "__main__":
191 # Add parser to get number of processors in space (have to do this here to enable automatic testing)
192 parser = ArgumentParser()
193 parser.add_argument("-n", "--nprocs_space", help='Specifies the number of processors in space', type=int)
194 args = parser.parse_args()
196 main(nprocs_space=args.nprocs_space)