← Back to branch summary

~mikael-mortensen/cbcpdesys/trunk

« back to all changes in this revision

Viewing changes to cbc/cfd/oasis/demo_doubly_periodic.py

  • Committer: Mikael Mortensen
  • Date: 2013-06-18 06:27:55 UTC
  • Revision ID: mikael.mortensen@gmail.com-20130618062755-rmwx447ozeyvq9ht
Adding method for projection on nonmatching meshes that works in parllel

Show diffs side-by-side

added added

removed removed

Lines of Context:
cbc/cfd/oasis/demo_doubly_periodic.py
103
103
Lx = 2.*pi
104
104
Ly = 2.
105
105
Lz = pi
106
 
Nx = 257
107
 
Ny = 193
108
 
Nz = 193
 
106
Nx = 25
 
107
Ny = 19
 
108
Nz = 19
109
109
 
110
110
mesh = BoxMesh(0., -Ly/2., -Lz/2., Lx, Ly/2., Lz/2., Nx, Ny, Nz)
111
111
# Create stretched mesh in y-direction
152
152
max_iter = 1                   # Pressure velocity iterations on given timestep
153
153
iters_on_first_timestep = 2    # Pressure velocity iterations on first timestep
154
154
max_error = 1e-6
155
 
check = 1000                     # print out info and save solution every check timestep 
 
155
check = 1                     # print out info and save solution every check timestep 
156
156
save_restart_file = 5000        # Saves two previous timesteps needed for a clean restart
157
157
    
158
158
# Specify body force
176
176
    folder = path.join(folder, str(max(map(eval, previous)) + 1))
177
177
 
178
178
MPI.barrier()
179
 
vtkfolder = path.join(folder, "VTK")
 
179
h5folder = path.join(folder, "HDF5")
180
180
statsfolder = path.join(folder, "Stats")
181
181
if MPI.process_number() == 0:
182
182
    try:
183
 
        makedirs(vtkfolder)
 
183
        makedirs(h5folder)
184
184
    except:
185
185
        pass
186
186
    try:
187
187
        makedirs(statsfolder)
188
188
    except:
189
189
        pass
190
 
u_stats_file = File(path.join(statsfolder, "umean.xml.gz"))
191
 
k_stats_file = File(path.join(statsfolder, "kmean.xml.gz"))
 
190
#u_stats_file = File(path.join(statsfolder, "umean.xml.gz"))
 
191
#k_stats_file = File(path.join(statsfolder, "kmean.xml.gz"))
192
192
 
193
193
#### Set a folder that contains xml.gz files of the solution. 
194
194
restart_folder = None        
253
253
    u0x = project(u0[0], V)
254
254
    u1x = project(u0[1], V)
255
255
    u2x = project(u0[2], V)
256
 
    y = interpolate(Expression("0.1335*((1+x[1])*(1-x[1]))"), V)
257
 
    q_['u0'].vector()[:] = y.vector()[:] 
 
256
    y = interpolate(Expression("x[1] > 0 ? 1-x[1] : 1+x[1]"), V)
 
257
    #y.vector()[y.vector().array() == 0] = 1e-12
 
258
    uu = project(1.25*(utau/0.41*ln(conditional(y<1e-12, 1.e-12, y)*utau/nu)+5.*utau), V)
 
259
    q_['u0'].vector()[:] = uu.vector()[:] 
258
260
    q_['u0'].vector().axpy(1.0, u0x.vector())
259
261
    q_['u1'].vector()[:] = u1x.vector()[:]
260
262
    q_['u2'].vector()[:] = u2x.vector()[:]
490
492
        if MPI.process_number()==0:
491
493
            info_green('Time = {0:2.4e}, timestep = {1:6d}, End time = {2:2.4e}'.format(t, tstep, T)) 
492
494
        newfolder = path.join(folder, 'timestep='+str(tstep))        
493
 
        #plot(u_[0], rescale=True)
 
495
        plot(u_[0], rescale=True)
494
496
        u1 = assemble(dot(u_, normal)*ds(1), exterior_facet_domains=facets)
495
497
        
496
498
        if MPI.process_number() == 0:
506
508
            for ui in u_components:
507
509
                newfile_1 = File(path.join(newfolder, ui + '_1.xml.gz'))
508
510
                newfile_1 << q_1[ui]
509
 
        
 
511
 
510
512
        stats.tovtk(0, filename=statsfolder+"/dump_mean_{}.vtk".format(tstep))
511
 
        voluviz(q_['u0'])
512
 
        voluviz.toh5_lowmem(0, tstep, filename=vtkfolder+"/snapshot_u0_{}.h5".format(tstep))
513
 
        voluviz.probes.clear()
514
 
        voluviz(q_['u1'])
515
 
        voluviz.toh5_lowmem(0, tstep, filename=vtkfolder+"/snapshot_u1_{}.h5".format(tstep))
516
 
        voluviz.probes.clear()
517
 
        voluviz(q_['u2'])
518
 
        voluviz.toh5_lowmem(0, tstep, filename=vtkfolder+"/snapshot_u2_{}.h5".format(tstep))
519
 
        voluviz.probes.clear()
 
513
        voluviz(q_['u0'])
 
514
        voluviz.toh5(0, tstep, filename=h5folder+"/snapshot_u0_{}.h5".format(tstep))
 
515
        voluviz.probes.clear()
 
516
        voluviz(q_['u1'])
 
517
        voluviz.toh5(0, tstep, filename=h5folder+"/snapshot_u1_{}.h5".format(tstep))
 
518
        voluviz.probes.clear()
 
519
        voluviz(q_['u2'])
 
520
        voluviz.toh5(0, tstep, filename=h5folder+"/snapshot_u2_{}.h5".format(tstep))
 
521
        voluviz.probes.clear()
520
522
        t1 = time.time()
521
523
        if MPI.process_number()==0:
522
524
            ff = open(newfolder+"/timeprstep_{}.txt".format(tottime/check), "w")