Skip to content
Snippets Groups Projects
Commit 14015815 authored by Franck Pérignon's avatar Franck Pérignon
Browse files

Remove useless files

parent 303235d5
No related branches found
No related tags found
No related merge requests found
......@@ -21,7 +21,7 @@ from parmepy.problem.navier_stokes import NSProblem
from parmepy.operator.monitors.printer import Printer
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
from dataTG import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
WITH_PROJ, OUTPUT_FREQ, FILENAME, simu
OUTPUT_FREQ, FILENAME, simu
from parmepy.constants import CONSERVATIVE
......@@ -96,7 +96,6 @@ poisson = Poisson(velo, vorti,
vorti: nbElem})
# projection=WITH_PROJ)
# Bridges between the different topologies in order to
# redistribute data.
# 1 -Advection to stretching
......@@ -110,7 +109,7 @@ distrStrPoiss = Redistribute([vorti, velo], stretch, poisson)
## Define the problem to solve
pb = NSProblem(operators=[advec, distrAdvStr, stretch,
distrStrPoiss,
diffusion, poisson], # , distrPoissEnergy],
diffusion, poisson],
simulation=simu, dumpFreq=-1)
## Setting solver to Problem (only operators for computational tasks)
......@@ -128,23 +127,16 @@ energy = Energy_enstrophy(velo, vorti,
## frequency=500,
## prefix='./res/TG_',
## ext='.vtk')
vorti.setTopoInit(topofft)
velo.setTopoInit(topofft)
pb.addMonitors([energy])
pb.setUp()
def run():
## Solve problem
pb.solve()
def printV(step):
for vd in velo.discreteFields.values():
print 'velo ', step, vd.norm()
for vd in vorti.discreteFields.values():
print 'vorti ', step, vd.norm()
## Solve problem
from parmepy.mpi import MPI
print "Start computation ..."
time = MPI.Wtime()
......
......@@ -15,9 +15,9 @@ WITH_PROJ = True
# post treatment ...
OUTPUT_FREQ = 1
FILENAME = './resp2_noproj/energy.dat'
FILENAME = './TG/energy2.dat'
# simulation parameters
from parmepy.problem.simulation import Simulation
simu = Simulation(tinit=0.0, tend=10., timeStep=0.01,
simu = Simulation(tinit=0.0, tend=1., timeStep=0.002,
iterMax=1000000)
import scitools.filetable as ft
import matplotlib.pyplot as plt
file1 = open("resp1/energy.dat")
ener = ft.read(file1)
time = ener[:, 0]
energy = ener[:, 1]
enstrophy = ener[:, 2]
#ratio = ener[:,3]
#dedt = ener[:,4]
#nuS = ener[:,5]
#nuES = ener[:,6]
size = energy.shape[0]
import numpy as np
em2 = np.zeros(size - 2)
em1 = np.zeros(size - 2)
em2[:] = energy[:-2]
em1[:] = energy[1:-1]
dt = np.zeros(size - 1)
dt[:] = time[1:] - time[:-1]
en = energy[2:]
dedt_calc = - (3. * en[:] - 4. * em1[:] + em2[:]) / (2. * dt[1:])
## file1 = open("resp1/energy.dat")
## ener = ft.read(file1)
## time = ener[:, 0]
## energy = ener[:, 1]
## enstrophy = ener[:, 2]
## #ratio = ener[:,3]
## #dedt = ener[:,4]
## #nuS = ener[:,5]
## #nuES = ener[:,6]
VISCOSITY = 1. / 1600.
nuEff = dedt_calc / enstrophy[2:]
ratio_calc = nuEff * 1600
nuS_calc = VISCOSITY * enstrophy[2:]
nuES_calc = nuEff * enstrophy[2:]
## nuS_calc = VISCOSITY * enstrophy[2:]
# multiproc #
file2 = open("resp2_noproj/energy.dat")
file2 = open("resp2_RK2/energy.dat")
ener2 = ft.read(file2)
time2 = ener2[:, 0]
energy2 = ener2[:, 1]
enstrophy2 = ener2[:, 2]
nuS_calc2 = enstrophy2[:] * VISCOSITY
size = energy2.shape[0]
em22 = np.zeros(size - 2)
em12 = np.zeros(size - 2)
em22[:] = energy2[:-2]
em12[:] = energy2[1:-1]
dt2 = np.zeros(size - 1)
dt2[:] = time2[1:] - time2[:-1]
en2 = energy2[2:]
nuS_calc2 = VISCOSITY * enstrophy2[2:]
plt.ioff()
plt.figure(1)
print time2[2:].shape
print nuS_calc2.shape
plt.plot(time[2:], nuS_calc)
plt.plot(time2[2:], nuS_calc2, 'o--')
#plt.plot(time, nuS_calc)
plt.plot(time2, enstrophy2, 'o--')
plt.show()
## size = energy.shape[0]
## import numpy as np
## em2 = np.zeros(size - 2)
## em1 = np.zeros(size - 2)
## em2[:] = energy[:-2]
## em1[:] = energy[1:-1]
## dt = np.zeros(size - 1)
## dt[:] = time[1:] - time[:-1]
## en = energy[2:]
## dedt_calc = - (3. * en[:] - 4. * em1[:] + em2[:]) / (2. * dt[1:])
## nuEff = dedt_calc / enstrophy[2:]
## ratio_calc = nuEff * 1600
## nuES_calc = nuEff * enstrophy[2:]
## size = energy2.shape[0]
## em22 = np.zeros(size - 2)
## em12 = np.zeros(size - 2)
## em22[:] = energy2[:-2]
## em12[:] = energy2[1:-1]
## dt2 = np.zeros(size - 1)
## dt2[:] = time2[1:] - time2[:-1]
## en2 = energy2[2:]
# -*- coding: utf-8 -*-
"""
@file fct2op.py
Convert differential operator for stretching into a function
to Discrete time integration
"""
from parmepy.numerics.differentialOperator import DifferentialOperator
import numpy as np
class Fct2Op(object):
"""
Fct2Op
make the link between differentialOperator and
the integration time scheme (Euler,RK,...)
"""
def __init__(self, field1, field2, method='', choice='', topology=None):
"""
Constructor.
store field2 to work with it in the apply function
@param field1 : vorticity field
@param field2 : velocity field or grad(U)
depending on 'choice' parameter
@param choice : can be 'gradV' for GradU.w or 'divWU' for div(wu)
"""
self.field1 = field1
self.field2 = field2
self.method = method
self.choice = choice
self.topology = topology
def apply(self, t, y):
"""
Apply operator.
Fct2Op.apply is used like a def function
"""
# self.synchro.apply()
if self.choice == 'gradV': # field2 = gradU
temp0 = y[0][...] * 0.
temp1 = y[0][...] * 0.
temp2 = y[0][...] * 0.
ind0a = self.topology.mesh.local_start[0]
ind0b = self.topology.mesh.local_end[0] + 1
ind1a = self.topology.mesh.local_start[1]
ind1b = self.topology.mesh.local_end[1] + 1
ind2a = self.topology.mesh.local_start[2]
ind2b = self.topology.mesh.local_end[2] + 1
ind0 = np.arange(ind0a, ind0b)
ind1 = np.arange(ind1a, ind1b)
ind2 = np.arange(ind2a, ind2b)
## Perform Matrix/Scalar product on local mesh
temp0[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
self.field2[3][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[4][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[5][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
self.field2[6][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[7][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
self.field2[8][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
# result = np.concatenate((
# np.array([temp0]),np.array([temp1]),np.array([temp2])))
return temp0, temp1, temp2
if self.choice == 'divWU': # field2 = velocity field
result = DifferentialOperator(
y, self.field2, self.method, self.choice).__call__()
return result[0], result[1], result[2]
if self.choice == 'hessianU': # field2 = velocity field
gradU, maxgersh = DifferentialOperator(
self.field1, self.field2,
self.method, choice='gradV').__call__()
GradUline1 = np.concatenate(([gradU[0]], [gradU[1]], [gradU[2]]))
GradUline2 = np.concatenate(([gradU[3]], [gradU[4]], [gradU[5]]))
GradUline3 = np.concatenate(([gradU[6]], [gradU[7]], [gradU[8]]))
result1, maxgersh = DifferentialOperator(
self.field1, GradUline1,
self.method, choice='gradV').__call__()
result2, maxgersh = DifferentialOperator(
self.field1, GradUline2,
self.method, choice='gradV').__call__()
result3, maxgersh = DifferentialOperator(
self.field1, GradUline3,
self.method, choice='gradV').__call__()
return result1, result2, result3
else:
print 'choice', choice, 'in Fct2Op is not defined'
def __str__(self):
s = "fct2Op(DiscreteOperator). " + DiscreteOperator.__str__(self)
return s
if __name__ == "__main__":
print __doc__
print "- Provided class : Fct2Op"
print Fct2Op.__doc__
......@@ -4,7 +4,6 @@
Discrete MultiPhase Rot Grad P
"""
from discrete import DiscreteOperator
from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.constants import np, debug
from parmepy.tools.timers import timed_function
......@@ -14,7 +13,7 @@ class DensityVisco_d(DiscreteOperator):
"""
@debug
def __init__(self, density, viscosity,
def __init__(self, density, viscosity,
method=None, densityVal=None, viscoVal=None):
"""
Constructor.
......
......@@ -4,8 +4,6 @@
Discrete MultiPhase Rot Grad P
"""
from discrete import DiscreteOperator
from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.numerics.fct2op import Fct2Op
from parmepy.constants import np, debug
from parmepy.tools.timers import timed_function
......
......@@ -3,10 +3,8 @@
@file compute_forces.py
Compute forces
"""
from parmepy.constants import debug, np, PARMES_REAL, ORDER, \
PARMES_INTEGER, PI, XDIR, YDIR, ZDIR
from parmepy.mpi.topology import Cartesian
from parmepy.numerics.fct2op import Fct2Op
from parmepy.constants import debug, np, ORDER, \
PARMES_INTEGER, PI, XDIR, YDIR, ZDIR
from parmepy.operator.updateGhosts import UpdateGhosts
from parmepy.numerics.differential_operations import DivStressTensor3D
from parmepy.numerics.finite_differences import FD_C_4
......
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment