-
Notifications
You must be signed in to change notification settings - Fork 2
/
ns
351 lines (301 loc) · 13.2 KB
/
ns
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
#!/usr/bin/env python
__author__ = "Anders Logg <[email protected]>"
__date__ = "2008-04-11"
__copyright__ = "Copyright (C) 2008-2010 " + __author__
__license__ = "GNU GPL version 3 or any later version"
import sys, time
from dolfin import set_log_active, parameters, list_timings, interactive, MPI, warning
from problems import Problem, problems
from solvers import Solver, solvers
# List of mesh sizes
mesh_sizes = [5,10,20,40,80,100]# 11, 16, 23, 32, 45, 64, 91, 128, 181, 256, 362]
#if len(sys.argv) > 1 and sys.argv[1] in ['channel', 'g2ref']:
# mesh_sizes = [7, 11, 15, 23, 31, 45, 63, 91, 127, 181, 255, 362]
master = MPI.process_number() == 0
import os
#if master:
# print
# print "PATH variable is:"
# print os.environ["PATH"]
# print
# Default options
OPTIONS = {
#simulations parameters
"problem_name" :sys.argv[1],
"solver_name" :sys.argv[2],
"segregated" :True,
"mesh_name" :"lshape_2",
"timesteps" :1000,
"cycles" :3,
"period" :951,
"viscosity" :0.0035,
"uOrder" :2,
"pOrder" :1,
"tOrder" :None,
#Boundary conditions information
"bc_in_ids" :None,
"bc_out_ids" :None,
"bc_in_profile" :"womersley", #parabolic for steady, womesley for pulsatile
"bc_in_data" :["FC_MCA"],
"bc_in_Qmean" :None,
"bc_out_type" :"zero_pressure",
#Compute various quantities
"compute_lplus" :True, #still needs implementations
#Saving/loading solutions
"save_frequency" :10000000,
"save_previous_ts" :True,
"restart" :False,
"restart_time" :1,
"restart_path" :"results/lshape_piso_lshape_2_ts1000_cycles3_uOrder1",
#Computing parameters
"processors" :8,
"node_type" :"debug",
"wct_hrs" :0,
"wct_mins" :2,
"time_out" :2850,
#case name
"case_name" :None,
"no_of_restart" :0,
"check_point" :None,
"check_point_frequency" :1,
##################################
"save_systole": False,
"max_steps": None,
"save_check_points": True,
"restart_ns":False,
"refinement_level": 0,
"dt_division": 0,
"save_solution": False,
"save_number": None, # max number of saves in total, at fixed frequency
"check_mem_usage": False,
"check_frequency": 10,
"save_solution_at_t=T": False,
"save_xml": False,
"plot_solution": False,
"plot_functional": False,
"compute_stress": False,
"compute_divergence": False,
"debug": True,
"timer": False,
"krylov_solver_absolute_tolerance": 1e-25,
"krylov_solver_relative_tolerance": 1e-12,
"krylov_solver_monitor_convergence": False,
"file_name": None, #Filename of the results file
"casename": None,
"postfix": None,
"dt": None,
"T": None,
# Used for restart_nsing a solution at a given time
"restart_ns_time": None,
"restart_ns_timestep": None,
"restart_ns_casename": None ,
"current_cycle" : None,
# Parameters for specific problems or solvers. If these become used
# in more places, please move them out of this section, and maybe
# change the default to None (i.e., move the default to the problem)
# if there is no sensible common default.
# ipcs solvers
"solver.u_tent": "gmres,hypre_euclid",
"solver.p_neumann": "gmres,hypre_amg",
"solver.p_dirichlet": "gmres,ml_amg",
"solver.p": None, # overrides neumann/dirichlet if given
"solver.u_corr": "bicgstab,hypre_euclid",
# grpc solver
# Parameters specific for Challenge
"test_case": 1,
"plot_probes": False,
"store_probes": False,
# Parameters specific for csf_flow_tapered
"waveform": "pulse",
"f1": 1.0,
#parameters for periodic boundary condition
"constrained_domain": None
}
OPTIONS["case_name"]=str(OPTIONS["problem_name"])+"_"+str(OPTIONS["solver_name"])+"_"+str(OPTIONS["mesh_name"])+"_ts"+str(OPTIONS["timesteps"])+"_cycles"+str(OPTIONS["cycles"])+"_uOrder"+str(OPTIONS["uOrder"])
#OPTIONS["restart_ns_path"]="./results/"+OPTIONS["case_name"]
def save_results(case_dir, problem, solver, num_dofs, cputime, wct, functional, dt_division, error):
"Save results to file."
# Print summary
if master:
print ""
print "Problem |", problem
print "Solver |", solver
print "Unknowns |", num_dofs
print "CPU time |", cputime
print "WCT time |", wct
print "Overhead |", wct - cputime
print "Functional |", functional
print "Error |", error
try:
if error != functional:
errpct = 100*error/functional
print "Error (%) |", errpct
except:
pass
# Print DOLFIN summary
set_log_active(True)
list_timings()
# Append to file
filename = os.path.join(case_dir, "results.log")
try:
with open(filename, "a") as file:
file.write("%s, %s, %s, %d, %.15g, %.15g, %.15g, %s, %s\n" %
(time.asctime(), problem, solver, num_dofs, cputime, wct,
functional, str(dt_division) , str(error)))
except Exception as e:
warning("Failed to write final results to log file: %s"%str(e))
def usage():
"Print usage"
if master:
print "Usage: ns problem solver\n\nAvailable problems:\n"
print "\n".join(" " + p for p in problems)
print "\nAvailable solvers:\n"
print "\n".join(" " + s for s in solvers)
def main(args):
"Parse command-line arguments and run solver"
# Check arguments
if not len(args) >= 2:
usage()
return 2
# Get problem and solver
problem_name, solver_name = args[:2]
# Get options
options = OPTIONS.copy()
for arg in args[2:]:
try:
key, value = arg.split("=")
#Try if there is a list
try:
list_value=value.split("[")[1].split("]")[0].split(",")
try: options[key]=[eval(i) for i in list_value]
except: options[key]=[i for i in list_value]
#Try if there is not a list
except:
try: options[key] = eval(value)
except: options[key] = str(value)
except:
if master: print "Warning: Unhandled command-line argument", arg
########### ADDED BY OWAIS ############################################################################
#define t_order, segregated and non-segredate, and save_previous time step
if options["tOrder"]==None and options["solver_name"]=="ipcs_ab_cn":
options["tOrder"]=2
if master: print "No tOrder defined; using ",str(options["tOrder"])," as default"
if options["tOrder"]==None and options["solver_name"]=="piso" :
options["tOrder"]=1
if master: print "No tOrder defined; using ",str(options["tOrder"])," as default"
if options["segregated"]==None and options["solver_name"]=="ipcs_ab_cn":
options["segregated"]=True
if master: print "No Segregated option defined; using ",options["segregated"]," as default"
if options["segregated"]==None and options["solver_name"]=="piso" :
options["segregated"]=False
if master: print "No Segregated option defined; using ",options["segregated"]," as default"
if options["save_previous_ts"]==None and options["solver_name"]=="ipcs_ab_cn":
options["save_previous_ts"]=True
if master: print "No save_previous_ts value defined; using ",options["save_previous_ts"]," as default"
if (options["save_previous_ts"]==None or options["save_previous_ts"]==True)and options["solver_name"]=="piso" :
options["save_previous_ts"]=False
if master: print "No save_previous_ts value defined; using ",options["save_previous_ts"]," as default"
#Obtain in a out ids from /.data/ids if None
try:
ids_infile=open("./data/ids")
ids_infile.readline()
for line in ids_infile:
line=line.split("\n")[0].split(" ")
if options["mesh_name"]==line[0]:
bc_in_ids=[]; bc_out_ids=[]; bc_in_Qmean=[]
if options["bc_in_ids"] ==None:
for i in line[1].split(","): bc_in_ids.append(int(i))
options["bc_in_ids"] = bc_in_ids
if master: print "No inlet ids defined; reading from ./data/ids", bc_in_ids
if options["bc_out_ids"] ==None:
for i in line[2].split(","): bc_out_ids.append(int(i))
options["bc_out_ids"]= bc_out_ids
if master: print "No outlet ids defined; reading from ./data/ids", bc_out_ids
if options["bc_in_Qmean"] ==None:
for i in line[3].split(","): bc_in_Qmean.append(float(i))
options["bc_in_Qmean"]= bc_in_Qmean
if master: print "No Qmean defined; reading from ./data/ids", bc_in_Qmean
except:
print "Inlet/Outets not reading from ./data/ids"
#######################################################################################################################
for key, value in options.iteritems():
if key.startswith("solver.") and isinstance(value, str):
options[key] = value.split(',')
# Set global DOLFIN parameters
parameters["form_compiler"]["cpp_optimize"] = True
parameters["krylov_solver"]["absolute_tolerance"] = options["krylov_solver_absolute_tolerance"]
parameters["krylov_solver"]["relative_tolerance"] = options["krylov_solver_relative_tolerance"]
parameters["krylov_solver"]["monitor_convergence"] = options["krylov_solver_monitor_convergence"]
# Set debug level
set_log_active(options["debug"])
# Set refinement level
options["N"] = mesh_sizes[options["refinement_level"]]
# Set refinement level
dt_division = str(options["dt_division"])
# Create problem and solver
if master: print "Making problem..."
problem = Problem(problem_name, options)
if master: print "Making solver..."
solver = Solver(solver_name, options)
if master:
print "Problem: " + str(problem)
print "Solver: " + str(solver)
# Lazy attempt at keeping compatibility with old way
if options["casename"] is None:
options["casename"] = solver.prefix(problem) + \
"_refinement_level_"+ str(options["refinement_level"])
if options["postfix"] is None:
options["casename_full"] = options["casename"]
else:
options["casename_full"] = options["casename"]+"_"+options["postfix"]
# Make directory for result files
if options["case_name"] is None:
casedir = os.path.join("results", options["casename_full"])
else:
casedir = os.path.join("results",options["case_name"])
options["casedir"] = casedir
if master:
if not os.path.exists(casedir):
os.mkdir(casedir)
os.mkdir(casedir+"/check_point")
print "Made case directory '%s'" % casedir
else:
print "NB! Directory already exists: '%s'." % casedir
# Solve problem with solver
wct = time.time()
if options["restart_ns_timestep"] is not None:
from solvers.restart_nsinfo import RestartInfo
restart_ns = RestartInfo(options["restart_ns_casename"] or options["casename_full"],
timestep=options["restart_ns_timestep"])
if master:
print "Restarting from timestep %d (time %.3g)" % (restart_ns.timestep, restart_ns.t)
u, p = solver.solve(problem, restart_ns)
elif options["restart_ns_time"] is not None:
from solvers.restart_nsinfo import RestartInfo
restart_ns = RestartInfo(options["restart_ns_casename"] or options["casename_full"],
t=options["restart_ns_time"])
if master:
print "Restarting from timestep %d (time %.3g)" % (restart_ns.timestep, restart_ns.t)
u, p = solver.solve(problem, restart_ns)
else:
u, p = solver.solve(problem)
# Compute elapsed time
wct = time.time() - wct
# Compute number of degrees of freedom
if options['segregated']:
num_dofs = sum(_u.vector().size() for _u in u) + p.vector().size()
else:
num_dofs = u.vector().size() + p.vector().size()
# Get functional value and error.
functional, error = solver.eval()
# Save results
save_results(casedir, problem, solver, num_dofs, solver.cputime(), wct, functional, dt_division, error)
if options['plot_solution']:
interactive()
return 0
if __name__ == "__main__":
args = sys.argv[1:]
if master: print "Running main with args", args
result = main(args)
if master: print "Result from main is", result
sys.exit(result)