Skip to content

Commit 9cb8ecd

Browse files
committed
Added input from @MichaelFlec. This closes #114
1 parent e100ee5 commit 9cb8ecd

File tree

6 files changed

+157
-11
lines changed

6 files changed

+157
-11
lines changed

.github/workflows/ci_pipeline.yml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ jobs:
1717

1818
- name: Add packages
1919
run: |
20+
sudo apt-get update
2021
sudo apt-get install texlive-latex-recommended texlive-fonts-recommended texlive-latex-extra cm-super dvipng
2122
2223
- name: Cache conda

pySDC/implementations/problem_classes/AllenCahn_Temp_MPIFFT.py

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ def eval_f(self, u, t):
122122
f.expl[..., 0] = self.fft.forward(tmpf)
123123

124124
f.impl[..., 1] = -self.params.D * self.K2 * u[..., 1]
125-
f.expl[..., 1] = -f.impl[..., 0] - f.expl[..., 0]
125+
f.expl[..., 1] = f.impl[..., 0] + f.expl[..., 0]
126126

127127
else:
128128

@@ -138,7 +138,7 @@ def eval_f(self, u, t):
138138
u_hat = self.fft.forward(u[..., 1])
139139
lap_u_hat = -self.params.D * self.K2 * u_hat
140140
f.impl[..., 1] = self.fft.backward(lap_u_hat, f.impl[..., 1])
141-
f.expl[..., 1] = -f.impl[..., 0] - f.expl[..., 0]
141+
f.expl[..., 1] = f.impl[..., 0] + f.expl[..., 0]
142142

143143
return f
144144

@@ -225,10 +225,10 @@ def circle_rand():
225225
def sine():
226226
tmp_me = newDistArray(self.fft, self.params.spectral)
227227
if self.params.spectral:
228-
tmp = np.sin(2 * np.pi * self.X[0]) * np.sin(2 * np.pi * self.X[1])
228+
tmp = 0.5 * (2.0 + 0.0 * np.sin(2 * np.pi * self.X[0]) * np.sin(2 * np.pi * self.X[1]))
229229
tmp_me[:] = self.fft.forward(tmp)
230230
else:
231-
tmp_me[:] = np.sin(2 * np.pi * self.X[0]) * np.sin(2 * np.pi * self.X[1])
231+
tmp_me[:] = 0.5 * (2.0 + 0.0 * np.sin(2 * np.pi * self.X[0]) * np.sin(2 * np.pi * self.X[1]))
232232
return tmp_me
233233

234234
assert t == 0, 'ERROR: u_exact only valid for t=0'
Binary file not shown.

pySDC/projects/AllenCahn_Bayreuth/run_temp_forcing_benchmark.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -61,8 +61,8 @@ def run_simulation(name=None, nprocs_space=None):
6161
problem_params['eps'] = [0.04]
6262
problem_params['radius'] = 0.25
6363
problem_params['TM'] = 1.0
64-
problem_params['D'] = 0.1
65-
problem_params['dw'] = [21.0]
64+
problem_params['D'] = 1.0
65+
problem_params['dw'] = [300.0]
6666
problem_params['comm'] = space_comm
6767
problem_params['name'] = name
6868
problem_params['init_type'] = 'circle_rand'
@@ -90,7 +90,7 @@ def run_simulation(name=None, nprocs_space=None):
9090

9191
# set time parameters
9292
t0 = 0.0
93-
Tend = 32 * 0.001
93+
Tend = 100 * 0.001
9494

9595
if space_rank == 0 and time_rank == 0:
9696
out = f'---------> Running {name} with {time_size} process(es) in time and {space_size} process(es) in space...'
Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
from argparse import ArgumentParser
2+
import numpy as np
3+
from mpi4py import MPI
4+
5+
from pySDC.helpers.stats_helper import filter_stats, sort_stats
6+
from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
7+
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
8+
from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
9+
from pySDC.implementations.problem_classes.AllenCahn_Temp_MPIFFT import allencahn_temp_imex
10+
11+
from pySDC.projects.AllenCahn_Bayreuth.AllenCahn_dump import dump
12+
13+
14+
def run_simulation(name='', spectral=None, nprocs_space=None):
15+
"""
16+
A test program to create reference data for the AC equation with temporal forcing
17+
Args:
18+
name (str): name of the run, will be used to distinguish different setups
19+
spectral (bool): run in real or spectral space
20+
nprocs_space (int): number of processors in space (None if serial)
21+
"""
22+
23+
# set MPI communicator
24+
comm = MPI.COMM_WORLD
25+
26+
world_rank = comm.Get_rank()
27+
world_size = comm.Get_size()
28+
29+
# split world communicator to create space-communicators
30+
if nprocs_space is not None:
31+
color = int(world_rank / nprocs_space)
32+
else:
33+
color = int(world_rank / 1)
34+
space_comm = comm.Split(color=color)
35+
space_rank = space_comm.Get_rank()
36+
space_size = space_comm.Get_size()
37+
38+
assert world_size == space_size, 'This script cannot run parallel-in-time with MPI, only spatial parallelism'
39+
40+
# initialize level parameters
41+
level_params = dict()
42+
level_params['restol'] = 1E-12
43+
level_params['dt'] = 1E-03
44+
level_params['nsweeps'] = [1]
45+
46+
# initialize sweeper parameters
47+
sweeper_params = dict()
48+
sweeper_params['collocation_class'] = CollGaussRadau_Right
49+
sweeper_params['num_nodes'] = [7]
50+
sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
51+
sweeper_params['initial_guess'] = 'spread'
52+
53+
# initialize problem parameters
54+
problem_params = dict()
55+
problem_params['L'] = 1.0
56+
problem_params['nvars'] = [(128, 128)]
57+
problem_params['eps'] = [0.03]
58+
problem_params['radius'] = 0.35682
59+
problem_params['TM'] = 1.0
60+
problem_params['D'] = 10.0
61+
problem_params['dw'] = [1.0]
62+
problem_params['comm'] = space_comm
63+
problem_params['name'] = name
64+
problem_params['init_type'] = 'circle'
65+
problem_params['spectral'] = spectral
66+
67+
# initialize step parameters
68+
step_params = dict()
69+
step_params['maxiter'] = 50
70+
71+
# initialize controller parameters
72+
controller_params = dict()
73+
controller_params['logger_level'] = 20 if space_rank == 0 else 99 # set level depending on rank
74+
controller_params['hook_class'] = dump
75+
76+
# fill description dictionary for easy step instantiation
77+
description = dict()
78+
description['problem_params'] = problem_params # pass problem parameters
79+
description['sweeper_class'] = imex_1st_order
80+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
81+
description['level_params'] = level_params # pass level parameters
82+
description['step_params'] = step_params # pass step parameters
83+
description['problem_class'] = allencahn_temp_imex
84+
85+
# set time parameters
86+
t0 = 0.0
87+
Tend = 100 * 0.001
88+
89+
if space_rank == 0:
90+
out = f'---------> Running {name} with spectral={spectral} and {space_size} process(es) in space...'
91+
print(out)
92+
93+
# instantiate controller
94+
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
95+
96+
# get initial values on finest level
97+
P = controller.MS[0].levels[0].prob
98+
uinit = P.u_exact(t0)
99+
100+
# call main function to get things done...
101+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
102+
103+
if space_rank == 0:
104+
105+
print()
106+
107+
# convert filtered statistics of iterations count, sorted by time
108+
iter_counts = sort_stats(filter_stats(stats, type='niter'), sortby='time')
109+
niters = np.mean(np.array([item[1] for item in iter_counts]))
110+
out = f'Mean number of iterations: {niters:.4f}'
111+
print(out)
112+
113+
# get setup time
114+
timing = sort_stats(filter_stats(stats, type='timing_setup'), sortby='time')
115+
out = f'Setup time: {timing[0][1]:.4f} sec.'
116+
print(out)
117+
118+
# get running time
119+
timing = sort_stats(filter_stats(stats, type='timing_run'), sortby='time')
120+
out = f'Time to solution: {timing[0][1]:.4f} sec.'
121+
print(out)
122+
123+
out = '...Done <---------\n'
124+
print(out)
125+
126+
127+
def main(nprocs_space=None):
128+
"""
129+
Little helper routine to run the whole thing
130+
Args:
131+
nprocs_space (int): number of processors in space (None if serial)
132+
"""
133+
name = 'AC-realistic-tempforce'
134+
run_simulation(name=name, spectral=False, nprocs_space=nprocs_space)
135+
# run_simulation(name=name, spectral=True, nprocs_space=nprocs_space)
136+
137+
138+
if __name__ == "__main__":
139+
140+
# Add parser to get number of processors in space (have to do this here to enable automatic testing)
141+
parser = ArgumentParser()
142+
parser.add_argument("-n", "--nprocs_space", help='Specifies the number of processors in space', type=int)
143+
args = parser.parse_args()
144+
145+
main(nprocs_space=args.nprocs_space)

pySDC/projects/AllenCahn_Bayreuth/run_temp_forcing_verification.py

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -169,8 +169,8 @@ def main(nprocs_space=None, cwd='.'):
169169
print(f'Error: {errors[-1]:6.4e}')
170170
print(f'Order of accuracy: {orders[-1]:4.2f}\n')
171171

172-
assert errors[2 + 1] < 1.4E-09, f'Errors are too high, got {errors[2 + 1]}'
173-
assert np.isclose(orders[3], 5, rtol=2E-02), f'Order of accuracy is not within tolerance, got {orders[3]}'
172+
assert errors[2 + 1] < 8E-10, f'Errors are too high, got {errors[2 + 1]}'
173+
assert np.isclose(orders[3], 5.3, rtol=2E-02), f'Order of accuracy is not within tolerance, got {orders[3]}'
174174

175175
print()
176176

@@ -183,8 +183,8 @@ def main(nprocs_space=None, cwd='.'):
183183
print(f'Error: {errors[-1]:6.4e}')
184184
print(f'Order of accuracy: {orders[-1]:4.2f}\n')
185185

186-
assert errors[2 + 1] < 1.4E-09, f'Errors are too high, got {errors[2 + 1]}'
187-
assert np.isclose(orders[1], 5, rtol=7E-02), f'Order of accuracy is not within tolerance, got {orders[1]}'
186+
assert errors[2 + 1] < 8E-10, f'Errors are too high, got {errors[2 + 1]}'
187+
assert np.isclose(orders[1], 4.6, rtol=7E-02), f'Order of accuracy is not within tolerance, got {orders[1]}'
188188

189189

190190
if __name__ == "__main__":

0 commit comments

Comments
 (0)