diff --git a/fftma_module/gen/FFTMA_TEST.py b/fftma_module/gen/FFTMA_TEST.py index a0823c4..d1c5d16 100644 --- a/fftma_module/gen/FFTMA_TEST.py +++ b/fftma_module/gen/FFTMA_TEST.py @@ -3,34 +3,32 @@ from refine import refine as ref import numpy as np from random import randint as rdi -N=128 +N = 128 for i in range(1): - nx, ny, nz = N,N,N - dx, dy, dz = 1.0, 1.0, 1.0 - seed= 1548762 #rdi(10000,99999) - var=1 - vario=2 - alpha=1 - lcx=2 - lcy=4 - lcz=16 - ap1x=1 - ap1y=0 - ap1z=0 - ap2x=0 - ap2y=1 - ap2z=0 - - v1 = (var, vario, alpha, lcx, lcy, lcz, ap1x, ap1y, ap1z, ap2x, ap2y, ap2z) - variograms = [v1] - - mean=15.3245987 - variance=3.5682389 - typ=1 - - k=gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ) - - np.save("out"+str(i)+".npy",ref(k,4,4,4)) - - + nx, ny, nz = N, N, N + dx, dy, dz = 1.0, 1.0, 1.0 + seed = 1548762 # rdi(10000,99999) + var = 1 + vario = 2 + alpha = 1 + lcx = 2 + lcy = 4 + lcz = 16 + ap1x = 1 + ap1y = 0 + ap1z = 0 + ap2x = 0 + ap2y = 1 + ap2z = 0 + + v1 = (var, vario, alpha, lcx, lcy, lcz, ap1x, ap1y, ap1z, ap2x, ap2y, ap2z) + variograms = [v1] + + mean = 15.3245987 + variance = 3.5682389 + typ = 1 + + k = gen(nx, ny, nz, dx, dy, dz, seed, variograms, mean, variance, typ) + + np.save("out" + str(i) + ".npy", ref(k, 4, 4, 4)) diff --git a/fftma_module/gen/setup.py b/fftma_module/gen/setup.py index 57b742c..5855a50 100644 --- a/fftma_module/gen/setup.py +++ b/fftma_module/gen/setup.py @@ -1,10 +1,74 @@ from distutils.core import setup, Extension -module_FFTMA = Extension('FFTMA', include_dirs = ['./include'],sources=["moduleFFTMA.c","./lib_src/Py_getvalues.c","./lib_src/Py_kgeneration.c","./lib_src/genlib.c","./lib_src/random.c","./lib_src/simpio.c","./lib_src/strlib.c","./lib_src/symtab.c","./lib_src/scanadt.c","./lib_src/stack.c","./lib_src/gammf.c","./lib_src/fftma.c","./lib_src/addstat.c","./lib_src/axes.c","./lib_src/cgrid.c","./lib_src/covariance.c","./lib_src/fourt.c","./lib_src/length.c","./lib_src/maxfactor.c","./lib_src/test_fact.c","./lib_src/cov_value.c","./lib_src/generate.c","./lib_src/gasdev.c","./lib_src/ran2.c","./lib_src/stable.c","./lib_src/gaussian.c","./lib_src/power.c","./lib_src/cubic.c","./lib_src/spherical.c","./lib_src/nugget.c","./lib_src/exponential.c","./lib_src/cardsin.c","./lib_src/nor2log.c","./lib_src/kgeneration.c","./lib_src/kgeneration2.c","./lib_src/fftma2.c","./lib_src/prebuild_gwn.c","./lib_src/build_real.c","./lib_src/addstat2.c","./lib_src/clean_real.c","./lib_src/pgeneration.c","./lib_src/pgeneration2.c","./lib_src/FFTPressure.c","./lib_src/FFTtest.c","./lib_src/build_pressure.c","./lib_src/build_velocity.c","./lib_src/total_pressure.c","./lib_src/total_velocity.c","./lib_src/clean_real2.c","./lib_src/waveVectorCompute3D.c","./lib_src/mat_vec.c","./lib_src/derivReal.c","./lib_src/inputdata.c","./lib_src/inputfiledata.c","./lib_src/debuginput.c","./lib_src/readdata.c","./lib_src/readfile_bin.c","./lib_src/writefile.c","./lib_src/writefile_bin.c","./lib_src/testmemory.c","./lib_src/testopenfile.c","./lib_src/readdata3.c"]) - - - +module_FFTMA = Extension( + "FFTMA", + include_dirs=["./include"], + sources=[ + "moduleFFTMA.c", + "./lib_src/Py_getvalues.c", + "./lib_src/Py_kgeneration.c", + "./lib_src/genlib.c", + "./lib_src/random.c", + "./lib_src/simpio.c", + "./lib_src/strlib.c", + "./lib_src/symtab.c", + "./lib_src/scanadt.c", + "./lib_src/stack.c", + "./lib_src/gammf.c", + "./lib_src/fftma.c", + "./lib_src/addstat.c", + "./lib_src/axes.c", + "./lib_src/cgrid.c", + "./lib_src/covariance.c", + "./lib_src/fourt.c", + "./lib_src/length.c", + "./lib_src/maxfactor.c", + "./lib_src/test_fact.c", + "./lib_src/cov_value.c", + "./lib_src/generate.c", + "./lib_src/gasdev.c", + "./lib_src/ran2.c", + "./lib_src/stable.c", + "./lib_src/gaussian.c", + "./lib_src/power.c", + "./lib_src/cubic.c", + "./lib_src/spherical.c", + "./lib_src/nugget.c", + "./lib_src/exponential.c", + "./lib_src/cardsin.c", + "./lib_src/nor2log.c", + "./lib_src/kgeneration.c", + "./lib_src/kgeneration2.c", + "./lib_src/fftma2.c", + "./lib_src/prebuild_gwn.c", + "./lib_src/build_real.c", + "./lib_src/addstat2.c", + "./lib_src/clean_real.c", + "./lib_src/pgeneration.c", + "./lib_src/pgeneration2.c", + "./lib_src/FFTPressure.c", + "./lib_src/FFTtest.c", + "./lib_src/build_pressure.c", + "./lib_src/build_velocity.c", + "./lib_src/total_pressure.c", + "./lib_src/total_velocity.c", + "./lib_src/clean_real2.c", + "./lib_src/waveVectorCompute3D.c", + "./lib_src/mat_vec.c", + "./lib_src/derivReal.c", + "./lib_src/inputdata.c", + "./lib_src/inputfiledata.c", + "./lib_src/debuginput.c", + "./lib_src/readdata.c", + "./lib_src/readfile_bin.c", + "./lib_src/writefile.c", + "./lib_src/writefile_bin.c", + "./lib_src/testmemory.c", + "./lib_src/testopenfile.c", + "./lib_src/readdata3.c", + ], +) setup(ext_modules=[module_FFTMA]) diff --git a/fftma_module/ref/setup.py b/fftma_module/ref/setup.py index 13fe47c..0393449 100644 --- a/fftma_module/ref/setup.py +++ b/fftma_module/ref/setup.py @@ -1,7 +1,7 @@ from distutils.core import setup, Extension -module = Extension('refine', sources=['FINALrefine.c']) +module = Extension("refine", sources=["FINALrefine.c"]) setup(ext_modules=[module]) diff --git a/fftma_module/ref/test.py b/fftma_module/ref/test.py index 48ab0e9..3827d09 100644 --- a/fftma_module/ref/test.py +++ b/fftma_module/ref/test.py @@ -2,15 +2,11 @@ from time import time import numpy as np import refine -size=420 -a=np.arange(size**3).astype('f8').reshape((size,size,size)) -ti=time() -b=refine.refine(a,2,2,2) -tf=time() -dt=tf-ti - -print a -print b -print dt +size = 420 +a = np.arange(size ** 3).astype("f8").reshape((size, size, size)) +ti = time() +b = refine.refine(a, 2, 2, 2) +tf = time() +dt = tf - ti raw_input("") diff --git a/fftma_module/testRes.py b/fftma_module/testRes.py index 502cae1..226da4e 100644 --- a/fftma_module/testRes.py +++ b/fftma_module/testRes.py @@ -2,63 +2,55 @@ import numpy as np import sys from refine import refine as ref -def get_p(pn, pdir, pprefix): - - p=np.load(pdir+pprefix+"0"+'.npy') - for i in range(1,pn): - p=np.concatenate((p,np.load(pdir+pprefix+str(i)+'.npy')),axis=0) +def get_p(pn, pdir, pprefix): - return p + p = np.load(pdir + pprefix + "0" + ".npy") + for i in range(1, pn): + p = np.concatenate((p, np.load(pdir + pprefix + str(i) + ".npy")), axis=0) + return p def get_k(pn, kdir, kprefix): - k=(np.load(kdir+kprefix+'0'+'.npy'))[1:-1,:,:] - for i in range(1,pn): - k=np.concatenate((k,(np.load(kdir+kprefix+str(i)+'.npy'))[1:-1,:,:]),axis=0) - return ref(k,2,2,2) - - - + k = (np.load(kdir + kprefix + "0" + ".npy"))[1:-1, :, :] + for i in range(1, pn): + k = np.concatenate( + (k, (np.load(kdir + kprefix + str(i) + ".npy"))[1:-1, :, :]), axis=0 + ) + return ref(k, 2, 2, 2) +def kef(P, K, i, j, k, pbc): + # tx=2*K[:,:,i]*K[:,:,i+1]/(K[:,:,i]+K[:,:,i+1]) + # ty=2*K[:,j,:]*K[:,j+1,:]/(K[:,j,:]+K[:,j+1,:]) + tz = 2 * K[k, :, :] * K[k + 1, :, :] / (K[k, :, :] + K[k + 1, :, :]) -def kef(P,K,i,j,k,pbc): - #tx=2*K[:,:,i]*K[:,:,i+1]/(K[:,:,i]+K[:,:,i+1]) - #ty=2*K[:,j,:]*K[:,j+1,:]/(K[:,j,:]+K[:,j+1,:]) - tz=2*K[k,:,:]*K[k+1,:,:]/(K[k,:,:]+K[k+1,:,:]) + # qx=tx*(P[:,:,i+1]-P[:,:,i]) + # qy=ty*(P[:,j+1,:]-P[:,j,:]) + qz = -tz * (P[k + 1, :, :] - P[k, :, :]) - #qx=tx*(P[:,:,i+1]-P[:,:,i]) - #qy=ty*(P[:,j+1,:]-P[:,j,:]) - qz=-tz*(P[k+1,:,:]-P[k,:,:]) - - kz=qz.sum()*(K.shape[0]+1)/(pbc*K.shape[1]*K.shape[2]) - return kz + kz = qz.sum() * (K.shape[0] + 1) / (pbc * K.shape[1] * K.shape[2]) + return kz def test(pn, kdir, pdir, kprefix, pprefix): - K=get_k(pn, kdir, kprefix) - print(K.shape) - P=get_p(pn, pdir, pprefix) - print(P) - print(P.shape) - print(kef(P,K,1,1,1,1000)) - return + K = get_k(pn, kdir, kprefix) + print(K.shape) + P = get_p(pn, pdir, pprefix) + print(P) + print(P.shape) + print(kef(P, K, 1, 1, 1, 1000)) + return -pn=4 -kdir="./test/" -pdir="./test/" -kprefix="k" -pprefix="P" +pn = 4 +kdir = "./test/" +pdir = "./test/" +kprefix = "k" +pprefix = "P" test(pn, kdir, pdir, kprefix, pprefix) - - - - - diff --git a/mpirunner.py b/mpirunner.py index 504a772..e209797 100755 --- a/mpirunner.py +++ b/mpirunner.py @@ -1,6 +1,7 @@ import numpy as np from mpi4py import MPI -#from tools.realization import realization + +# from tools.realization import realization from tools.generation.config import DotheLoop, get_config import os import sys @@ -8,59 +9,67 @@ from tools.Prealization import realization from utilities.conditional_decorator import * from memory_profiler import profile -CONFIG_FILE_PATH = 'config.ini' if 'CONFIG_FILE_PATH' not in os.environ else os.environ['CONFIG_FILE_PATH'] -IS_TEST = False if 'TEST' not in os.environ else True +CONFIG_FILE_PATH = ( + "config.ini" + if "CONFIG_FILE_PATH" not in os.environ + else os.environ["CONFIG_FILE_PATH"] +) +IS_TEST = False if "TEST" not in os.environ else True + def main(): - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - pn = comm.Get_size() + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + pn = comm.Get_size() + + if pn == 1: + sequential() + return + if rank == 0: + manager() + else: + worker() + return - if pn==1: - sequential() - return - if rank==0: - manager() - else: - worker() - return @conditional_decorator(profile, IS_TEST) def sequential(): - comm = MPI.COMM_WORLD - conffile = CONFIG_FILE_PATH - parser,iterables = get_config(conffile) - njobs = DotheLoop(-1,parser,iterables) - start_job=0 - for job in range(start_job,njobs): - realization(job) - return + comm = MPI.COMM_WORLD + conffile = CONFIG_FILE_PATH + parser, iterables = get_config(conffile) + njobs = DotheLoop(-1, parser, iterables) + start_job = 0 + for job in range(start_job, njobs): + realization(job) + return + def manager(): - comm = MPI.COMM_WORLD - conffile = CONFIG_FILE_PATH - parser,iterables = get_config(conffile) - njobs = DotheLoop(-1,parser,iterables) - start_job=0 - for job in range(start_job,njobs): - dest=comm.recv(source=MPI.ANY_SOURCE) - comm.send(job,dest=dest) - for i in range(comm.Get_size()-1): - dest=comm.recv(source=MPI.ANY_SOURCE) - comm.send(-1,dest=dest) + comm = MPI.COMM_WORLD + conffile = CONFIG_FILE_PATH + parser, iterables = get_config(conffile) + njobs = DotheLoop(-1, parser, iterables) + start_job = 0 + for job in range(start_job, njobs): + dest = comm.recv(source=MPI.ANY_SOURCE) + comm.send(job, dest=dest) + for i in range(comm.Get_size() - 1): + dest = comm.recv(source=MPI.ANY_SOURCE) + comm.send(-1, dest=dest) + + return - return @conditional_decorator(profile, IS_TEST) def worker(): - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - job=1 - while job!=-1: - comm.send(rank,dest=0) - job = comm.recv(source=0) - realization(job) - return + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + job = 1 + while job != -1: + comm.send(rank, dest=0) + job = comm.recv(source=0) + realization(job) + return main() diff --git a/tests/integration/test.py b/tests/integration/test.py index 997a989..c03c211 100644 --- a/tests/integration/test.py +++ b/tests/integration/test.py @@ -4,6 +4,7 @@ import numpy as np import unittest from numpy.lib.function_base import diff + def find_relative_errors(path_original, path): binary_original = np.load(path_original) binary = np.load(path) @@ -18,22 +19,31 @@ def find_relative_errors(path_original, path): for y in range(len(diffs)): for z in range(len(diffs)): if type(diffs[x][y][z]) != type([]): - relative_error = 0 if binary_original[x][y][z] == 0 else diffs[x][y][z] / binary_original[x][y][z] + relative_error = ( + 0 + if binary_original[x][y][z] == 0 + else diffs[x][y][z] / binary_original[x][y][z] + ) relative_errors.append(abs(relative_error)) else: for w in range(len(diffs)): - relative_error = 0 if binary_original[x][y][z][w] == 0 else diffs[x][y][z][w] / binary_original[x][y][z][w] + relative_error = ( + 0 + if binary_original[x][y][z][w] == 0 + else diffs[x][y][z][w] / binary_original[x][y][z][w] + ) relative_errors.append(abs(relative_error)) return relative_errors -BINARIES = ['Cmap', 'D', 'P', 'V', 'k'] -class TestIntegration(unittest.TestCase): +BINARIES = ["Cmap", "D", "P", "V", "k"] + +class TestIntegration(unittest.TestCase): @classmethod def setUpClass(cls): - os.chdir('../..') + os.chdir("../..") config_file = os.path.abspath("./tests/integration/conf_test.ini") os.system(f"CONFIG_FILE_PATH={config_file} mpirun -np 1 python3 mpirunner.py") @@ -43,16 +53,21 @@ class TestIntegration(unittest.TestCase): for i in range(90): for binary in BINARIES: - path = './tests/integration/tmp_output/{}/{}.npy'.format(i, binary) - path_original = './test_loop/{}/{}.npy'.format(i, binary) + path = "./tests/integration/tmp_output/{}/{}.npy".format(i, binary) + path_original = "./test_loop/{}/{}.npy".format(i, binary) relative_errors = find_relative_errors(path_original, path) binary_results[binary].append(relative_errors) cls.binary_stats = {} for binary in binary_results: - binary_results[binary] = [item for sublist in binary_results[binary] for item in sublist] + binary_results[binary] = [ + item for sublist in binary_results[binary] for item in sublist + ] if len(binary_results[binary]) != 0: - cls.binary_stats[binary] = {"max": max(binary_results[binary]), "avg": sum(binary_results[binary]) / len(binary_results[binary])} + cls.binary_stats[binary] = { + "max": max(binary_results[binary]), + "avg": sum(binary_results[binary]) / len(binary_results[binary]), + } @classmethod def tearDownClass(cls): @@ -84,5 +99,5 @@ class TestIntegration(unittest.TestCase): self.assertLess(V_stats["avg"], 0.05) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/tests/performance/connectivity.py b/tests/performance/connectivity.py index c71634c..a1f84a8 100644 --- a/tests/performance/connectivity.py +++ b/tests/performance/connectivity.py @@ -1,7 +1,7 @@ import os from benchmarker import Benchmarker -os.chdir('../..') +os.chdir("../..") config_gen_file_64 = os.path.abspath("./tests/performance/conf_gen_64.ini") config_conn_file_64 = os.path.abspath("./tests/performance/conf_conn_64.ini") @@ -13,7 +13,7 @@ index_1 = 0 index_8 = 0 -''' +""" Esta etapa tarda mucho tiempo y no es muy independiente de la generación de medios. Si se generan medios con los parámetros dados: [Iterables] @@ -26,32 +26,39 @@ Se generan 90 medios: 15 (p[2]) * 2 (seeds[1]) * 3 (len(connectivity)) * 1 (len( Pero si se toman esos medios generados y se aplica solo la etapa de conectividad Se calcula la conectividad sobre 6 medios: 2 (seeds[1]) * 3 (len(connectivity)) * 1 (len(variances))* 1 (len(lc)) Solucion: marcar en la etapa de generacion binary = yes -> esta bien esto? -''' +""" with Benchmarker() as bench: for i in range(len(CONN_CONFIG_FILES)): - size = 2**(6+i) + size = 2 ** (6 + i) @bench(f"Connectivity 1 core with size {size}") def _(bm): global index_1 - os.system(f"CONFIG_FILE_PATH={GEN_CONFIG_FILES[index_1]} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={GEN_CONFIG_FILES[index_1]} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py" + ) with bm: - os.system(f"CONFIG_FILE_PATH={CONN_CONFIG_FILES[index_1]} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={CONN_CONFIG_FILES[index_1]} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py" + ) ## teardown os.system("rm -rf ./tests/performance/tmp_gen_output") - index_1 +=1 + index_1 += 1 @bench(f"Connectivity 8 core with size {size}") def _(bm): global index_8 - os.system(f"CONFIG_FILE_PATH={GEN_CONFIG_FILES[index_8]} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={GEN_CONFIG_FILES[index_8]} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py" + ) with bm: - os.system(f"CONFIG_FILE_PATH={CONN_CONFIG_FILES[index_8]} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={CONN_CONFIG_FILES[index_8]} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py" + ) ## teardown os.system("rm -rf ./tests/performance/tmp_gen_output") - index_8 +=1 - + index_8 += 1 diff --git a/tests/performance/generation.py b/tests/performance/generation.py index 837467d..c5e214d 100644 --- a/tests/performance/generation.py +++ b/tests/performance/generation.py @@ -1,7 +1,7 @@ import os from benchmarker import Benchmarker -os.chdir('../..') +os.chdir("../..") config_file_64 = os.path.abspath("./tests/performance/conf_gen_64.ini") @@ -18,26 +18,30 @@ index_8 = 0 with Benchmarker() as bench: for i in range(len(CONFIG_FILES)): - size = 2**(6+i) + size = 2 ** (6 + i) @bench(f"generation 1 core {size} tamaño") def _(bm): global index_1 config_file = CONFIG_FILES[index_1] with bm: - os.system(f"CONFIG_FILE_PATH={config_file} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={config_file} TEST=True mpirun -oversubscribe -np 1 python3 mpirunner.py" + ) ## teardown os.system("rm -rf ./tests/performance/tmp_gen_output") - index_1 +=1 + index_1 += 1 @bench(f"generation 8 core {size} tamaño") def _(bm): global index_8 config_file = CONFIG_FILES[index_8] with bm: - os.system(f"CONFIG_FILE_PATH={config_file} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py") + os.system( + f"CONFIG_FILE_PATH={config_file} TEST=True mpirun -oversubscribe -np 8 python3 mpirunner.py" + ) ## teardown os.system("rm -rf ./tests/performance/tmp_gen_output") - index_8 +=1 \ No newline at end of file + index_8 += 1 diff --git a/tools/Prealization.py b/tools/Prealization.py index 21a566f..16a3aae 100755 --- a/tools/Prealization.py +++ b/tools/Prealization.py @@ -1,5 +1,5 @@ import numpy as np -from tools.generation.config import DotheLoop, get_config +from tools.generation.config import DotheLoop, get_config import os import sys from mpi4py import MPI @@ -10,73 +10,75 @@ from tools.solver.comp_Kperm_scale import comp_kperm_sub from tools.solver.Ndar import PetscP from tools.generation.fftma_gen import fftmaGenerator -CONFIG_FILE_PATH = 'config.ini' if 'CONFIG_FILE_PATH' not in os.environ else os.environ['CONFIG_FILE_PATH'] - -def realization(job): - - if job==-1: - return - - - conffile=CONFIG_FILE_PATH - parser, iterables = get_config(conffile) - start_job=int(parser.get('General',"startJob")) - - if job1: - icomm=MPI.COMM_SELF.Spawn(sys.executable, args=['./tools/solver/Ndar.py',datadir,str(ref),'0',Rtol,'1'], maxprocs=n_p) - icomm.Disconnect() - else: - PetscP(datadir,ref,'0',True,float(Rtol),0) - - - compkperm=parser.get('K-Postprocess',"kperm") - if compkperm!='no': - #print('start kperm') - comp_kperm_sub(parser,datadir,nr) - #print('finished job ' +str(job)) - - postP=parser.get('K-Postprocess',"postprocess") - if postP!='no': - comp_postKeff(parser,datadir,nr) - - - return - -def create_dir(datadir,job): - - try: - os.makedirs(datadir) - except: - print('Warning: Unable to create dir job: '+str(job)) - - return - +CONFIG_FILE_PATH = ( + "config.ini" + if "CONFIG_FILE_PATH" not in os.environ + else os.environ["CONFIG_FILE_PATH"] +) +def realization(job): + if job == -1: + return + + conffile = CONFIG_FILE_PATH + parser, iterables = get_config(conffile) + start_job = int(parser.get("General", "startJob")) + + if job < start_job: + return + + rdir = "./" + parser.get("General", "simDir") + "/" + datadir = rdir + str(job) + "/" + create_dir(datadir, job) + if job == 0: + copyfile(conffile, rdir + "config.ini") + + genera = parser.get("Generation", "genera") + if genera != "no": + fftmaGenerator(datadir, job, CONFIG_FILE_PATH) + # os.system('CONFIG_FILE_PATH=' + CONFIG_FILE_PATH + ' python3 ./tools/generation/fftma_gen.py ' + datadir +' ' + str(job)) + + nr = DotheLoop(job, parser, iterables)[3] - iterables["seeds"][0] + Cconec = parser.get("Connectivity", "conec") + if Cconec != "no": + comp_connec(parser, datadir, nr) + + n_p = int(parser.get("Solver", "num_of_cores")) + ref = int(parser.get("Solver", "ref")) + solv = parser.get("Solver", "solve") + Rtol = parser.get("Solver", "rtol") + + if solv != "no": + if n_p > 1: + icomm = MPI.COMM_SELF.Spawn( + sys.executable, + args=["./tools/solver/Ndar.py", datadir, str(ref), "0", Rtol, "1"], + maxprocs=n_p, + ) + icomm.Disconnect() + else: + PetscP(datadir, ref, "0", True, float(Rtol), 0) + + compkperm = parser.get("K-Postprocess", "kperm") + if compkperm != "no": + # print('start kperm') + comp_kperm_sub(parser, datadir, nr) + # print('finished job ' +str(job)) + + postP = parser.get("K-Postprocess", "postprocess") + if postP != "no": + comp_postKeff(parser, datadir, nr) + + return + + +def create_dir(datadir, job): + + try: + os.makedirs(datadir) + except: + print("Warning: Unable to create dir job: " + str(job)) + + return diff --git a/tools/connec/JoinCmaps.py b/tools/connec/JoinCmaps.py index af20e5e..0c17fff 100755 --- a/tools/connec/JoinCmaps.py +++ b/tools/connec/JoinCmaps.py @@ -1,136 +1,142 @@ import numpy as np -def joinCmapX(cmap1,cmap2): - nclus1 = np.max(cmap1) - cmap2=np.where(cmap2!=0,cmap2+nclus1,0) - - old_nclus=0 - new_nclus=1 - - while new_nclus!= old_nclus: - - old_nclus=new_nclus - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap2=np.where(cmap2==cmap2[0,i,j],cmap1[-1,i,j],cmap2) - - - - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap1=np.where(cmap1==cmap1[-1,i,j],cmap2[0,i,j],cmap1) - - cmap=np.append(cmap1,cmap2,axis=0) - y = np.bincount(cmap.reshape(-1).astype(int)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - new_nclus=cf.shape[0] #cantidad de clusters - #print(new_nclus) - - return cmap - - -def joinCmapY(cmap1,cmap2): - - nclus1 = np.max(cmap1) - cmap2=np.where(cmap2!=0,cmap2+nclus1,0) - - old_nclus=0 - new_nclus=1 - - while new_nclus!= old_nclus: - - old_nclus=new_nclus - for i in range(cmap1.shape[0]): - for j in range(cmap1.shape[2]): - if cmap1[i,-1,j] != 0 and cmap2[i,0,j] !=0: - if cmap1[i,-1,j] != cmap2[i,0,j]: - cmap2=np.where(cmap2==cmap2[i,0,j],cmap1[i,-1,j],cmap2) - - - - for i in range(cmap1.shape[0]): - for j in range(cmap1.shape[2]): - if cmap1[i,-1,j] != 0 and cmap2[i,0,j] !=0: - if cmap1[i,-1,j] != cmap2[i,0,j]: - cmap1=np.where(cmap1==cmap1[i,-1,j],cmap2[i,0,j],cmap1) - - cmap=np.append(cmap1,cmap2,axis=1) - y = np.bincount(cmap.reshape(-1).astype(int)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - new_nclus=cf.shape[0] #cantidad de clusters - #print(new_nclus) - - return cmap - - -def joinCmapZ(cmap1,cmap2): - - nclus1 = np.max(cmap1) - cmap2=np.where(cmap2!=0,cmap2+nclus1,0) - - old_nclus=0 - new_nclus=1 - - while new_nclus!= old_nclus: - - old_nclus=new_nclus - for i in range(cmap1.shape[0]): - for j in range(cmap1.shape[1]): - if cmap1[i,j,-1] != 0 and cmap2[i,j,0] !=0: - if cmap1[i,j,-1] != cmap2[i,j,0]: - cmap2=np.where(cmap2==cmap2[i,j,0],cmap1[i,j,-1],cmap2) - - - - for i in range(cmap1.shape[0]): - for j in range(cmap1.shape[1]): - if cmap1[i,j,-1] != 0 and cmap2[i,j,0] !=0: - if cmap1[i,j,-1] != cmap2[i,j,0]: - cmap1=np.where(cmap1==cmap1[i,j,-1],cmap2[i,j,0],cmap1) - - cmap=np.append(cmap1,cmap2,axis=2) - y = np.bincount(cmap.reshape(-1).astype(int)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - new_nclus=cf.shape[0] #cantidad de clusters - #print(new_nclus) - - return cmap - -def joinBox(vec,join_y,join_z): - - Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2] - nx = Nx//2 - ny,nz=Ny, Nz - if join_y: - ny=Ny//2 - if join_z: - nz=Nz//2 - - vec[:,:ny,:nz] = joinCmapX(vec[:nx,:ny,:nz],vec[nx:,:ny,:nz]) - if not join_z and not join_y: - return vec - - if join_y: - vec[:,ny:,:nz] = joinCmapX(vec[:nx,ny:,:nz],vec[nx:,ny:,:nz]) - if join_z: - vec[:,:ny,nz:] = joinCmapX(vec[:nx,:ny,nz:],vec[nx:,:ny,nz:]) - if join_z and join_y: - vec[:,ny:,nz:] = joinCmapX(vec[:nx,ny:,nz:],vec[nx:,ny:,nz:]) - if join_y: - vec[:,:,:nz] = joinCmapY(vec[:,:ny,:nz],vec[:,ny:,:nz]) - if join_z: - if join_y: - vec[:,:,nz:] = joinCmapY(vec[:,:ny,nz:],vec[:,ny:,nz:]) - vec[:,:,:] = joinCmapZ(vec[:,:,:nz],vec[:,:,nz:]) - - return vec +def joinCmapX(cmap1, cmap2): + nclus1 = np.max(cmap1) + cmap2 = np.where(cmap2 != 0, cmap2 + nclus1, 0) + + old_nclus = 0 + new_nclus = 1 + + while new_nclus != old_nclus: + old_nclus = new_nclus + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap2 = np.where( + cmap2 == cmap2[0, i, j], cmap1[-1, i, j], cmap2 + ) + + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap1 = np.where( + cmap1 == cmap1[-1, i, j], cmap2[0, i, j], cmap1 + ) + + cmap = np.append(cmap1, cmap2, axis=0) + y = np.bincount(cmap.reshape(-1).astype(int)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + new_nclus = cf.shape[0] # cantidad de clusters + # print(new_nclus) + + return cmap + + +def joinCmapY(cmap1, cmap2): + + nclus1 = np.max(cmap1) + cmap2 = np.where(cmap2 != 0, cmap2 + nclus1, 0) + + old_nclus = 0 + new_nclus = 1 + + while new_nclus != old_nclus: + + old_nclus = new_nclus + for i in range(cmap1.shape[0]): + for j in range(cmap1.shape[2]): + if cmap1[i, -1, j] != 0 and cmap2[i, 0, j] != 0: + if cmap1[i, -1, j] != cmap2[i, 0, j]: + cmap2 = np.where( + cmap2 == cmap2[i, 0, j], cmap1[i, -1, j], cmap2 + ) + + for i in range(cmap1.shape[0]): + for j in range(cmap1.shape[2]): + if cmap1[i, -1, j] != 0 and cmap2[i, 0, j] != 0: + if cmap1[i, -1, j] != cmap2[i, 0, j]: + cmap1 = np.where( + cmap1 == cmap1[i, -1, j], cmap2[i, 0, j], cmap1 + ) + + cmap = np.append(cmap1, cmap2, axis=1) + y = np.bincount(cmap.reshape(-1).astype(int)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + new_nclus = cf.shape[0] # cantidad de clusters + # print(new_nclus) + + return cmap + + +def joinCmapZ(cmap1, cmap2): + + nclus1 = np.max(cmap1) + cmap2 = np.where(cmap2 != 0, cmap2 + nclus1, 0) + + old_nclus = 0 + new_nclus = 1 + + while new_nclus != old_nclus: + + old_nclus = new_nclus + for i in range(cmap1.shape[0]): + for j in range(cmap1.shape[1]): + if cmap1[i, j, -1] != 0 and cmap2[i, j, 0] != 0: + if cmap1[i, j, -1] != cmap2[i, j, 0]: + cmap2 = np.where( + cmap2 == cmap2[i, j, 0], cmap1[i, j, -1], cmap2 + ) + + for i in range(cmap1.shape[0]): + for j in range(cmap1.shape[1]): + if cmap1[i, j, -1] != 0 and cmap2[i, j, 0] != 0: + if cmap1[i, j, -1] != cmap2[i, j, 0]: + cmap1 = np.where( + cmap1 == cmap1[i, j, -1], cmap2[i, j, 0], cmap1 + ) + + cmap = np.append(cmap1, cmap2, axis=2) + y = np.bincount(cmap.reshape(-1).astype(int)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + new_nclus = cf.shape[0] # cantidad de clusters + # print(new_nclus) + + return cmap + + +def joinBox(vec, join_y, join_z): + + Nx, Ny, Nz = vec.shape[0], vec.shape[1], vec.shape[2] + nx = Nx // 2 + ny, nz = Ny, Nz + if join_y: + ny = Ny // 2 + if join_z: + nz = Nz // 2 + + vec[:, :ny, :nz] = joinCmapX(vec[:nx, :ny, :nz], vec[nx:, :ny, :nz]) + if not join_z and not join_y: + return vec + + if join_y: + vec[:, ny:, :nz] = joinCmapX(vec[:nx, ny:, :nz], vec[nx:, ny:, :nz]) + if join_z: + vec[:, :ny, nz:] = joinCmapX(vec[:nx, :ny, nz:], vec[nx:, :ny, nz:]) + if join_z and join_y: + vec[:, ny:, nz:] = joinCmapX(vec[:nx, ny:, nz:], vec[nx:, ny:, nz:]) + if join_y: + vec[:, :, :nz] = joinCmapY(vec[:, :ny, :nz], vec[:, ny:, :nz]) + if join_z: + if join_y: + vec[:, :, nz:] = joinCmapY(vec[:, :ny, nz:], vec[:, ny:, nz:]) + vec[:, :, :] = joinCmapZ(vec[:, :, :nz], vec[:, :, nz:]) + + return vec diff --git a/tools/connec/PostConec (copy).py b/tools/connec/PostConec (copy).py index b390256..be4e818 100755 --- a/tools/connec/PostConec (copy).py +++ b/tools/connec/PostConec (copy).py @@ -6,272 +6,301 @@ import os import collections - -def ConnecInd(cmap,scales,datadir): - - datadir=datadir+'ConnectivityMetrics/' - try: - os.makedirs(datadir) - except: - nada=0 - - for scale in scales: - res=dict() - res=doforsubS_computeCmap(res,cmap,scale,postConec) - np.save(datadir+str(scale)+'.npy',res) - - return - -def doforsubS_computeCmap(res,cmap,l,funpost): - - - L=cmap.shape[0] - Nx, Ny,Nz=cmap.shape[0],cmap.shape[1],cmap.shape[2] - - nblx=Nx//l #for each dimension - nbly=Ny//l - if cmap.shape[2]==1: - lz=1 - nblz=1 - else: - lz=l - nblz=Nz//l - - keys=funpost(np.array([]),res,0,0) - - for key in keys: - res[key]=np.zeros((nblx,nbly,nblz)) - - for i in range(nblx): - for j in range(nbly): - for k in range(nblz): - res=funpost(cmap[i*l:(i+1)*l,j*l:(j+1)*l,k*l:(k+1)*lz],res,(i,j,k),1) - - return res - -def postConec(cmap,results,ind,flag): - - - if flag==0: - keys=[] - keys+=['PPHA'] - keys+=['VOLALE'] - keys+=['ZNCC'] - keys+=['GAMMA'] - keys+=['spanning', 'npz', 'npy', 'npx'] - keys+=['Plen','S','P'] - return keys - - - dim=3 - if cmap.shape[2]==1: - cmap=cmap[:,:,0] - dim=2 - - y = np.bincount(cmap.reshape(-1)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - - if cf[0,0]==0: - cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - - if cf.shape[0]>0: - - - spanning, pclusZ, pclusY, pclusX =get_perco(cmap,dim) - plen=Plen(spanning,cmap,cf,dim) - nper=np.sum(cf[:,1]) #num de celdas permeables - nclus=cf.shape[0] #cantidad de clusters - - results['PPHA'][ind]=nper/np.size(cmap) #ppha - results['VOLALE'][ind]=np.max(cf[:,1])/nper #volale #corregido va entre [0,p] - results['ZNCC'][ind]=nclus #zncc - results['GAMMA'][ind]=np.sum(cf[:,1]**2)/np.size(cmap)/nper #gamma, recordar zintcc =gamma*p - results['spanning'][ind],results['npz'][ind], results['npy'][ind], results['npx'][ind]=spanning, len(pclusZ), len(pclusY), len(pclusX) - results['Plen'][ind],results['S'][ind],results['P'][ind] = plen[0],plen[1],plen[2] - - - if cf.shape[0]==0: - for key in keys: - results[key][ind]=0 - return results - - -#ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA - - -def get_pos2D(cmap,cdis): - - Ns=cdis.shape[0] - pos=dict() - i=0 - for cnum in cdis[:,0]: - pos[cnum]=np.zeros((cdis[i,1]+1,2)) #+1 porque uso de flag - i+=1 - - for i in range(cmap.shape[0]): - for j in range(cmap.shape[1]): - if cmap[i,j] != 0: - flag=int(pos[cmap[i,j]][0,0])+1 - pos[cmap[i,j]][0,0]=flag - pos[cmap[i,j]][flag,0]=i - pos[cmap[i,j]][flag,1]=j - return pos - - -def get_pos3D(cmap,cdis): - - Ns=cdis.shape[0] - pos=dict() - i=0 - for cnum in cdis[:,0]: - pos[cnum]=np.zeros((cdis[i,1]+1,3)) - i+=1 - for i in range(cmap.shape[0]): - for j in range(cmap.shape[1]): - for k in range(cmap.shape[2]): - - if cmap[i,j,k] != 0: - flag=int(pos[cmap[i,j,k]][0,0])+1 - pos[cmap[i,j,k]][0,0]=flag - pos[cmap[i,j,k]][flag,0]=i - pos[cmap[i,j,k]][flag,1]=j - pos[cmap[i,j,k]][flag,2]=k - - - return pos - -def Plen(spannng,cmap,cdis,dim): - - if dim==2: - return P_len2D(spannng,cmap,cdis) - if dim==3: - return P_len3D(spannng,cmap,cdis) - return [] - -def P_len2D(spanning,cmap,cdis): - - pos = get_pos2D(cmap,cdis) - #print(summary['NpcY'],summary['NpcX'],summary['PPHA']) - - den=0 - num=0 - - nperm=np.sum(cdis[:,1]) - if spanning > 0: - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - else: - P=0 - - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - - - -def P_len3D(spanning,cmap,cdis): - - - pos = get_pos3D(cmap,cdis) - #print(summary['NpcY'],summary['NpcX'],summary['PPHA']) - - den=0 - num=0 - - nperm=np.sum(cdis[:,1]) - if spanning > 0: - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - else: - P=0 - - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy, mposz = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]), np.mean(pos[cnum][1:,2]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2+(pos[cnum][1:,2]-mposz)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - - - -def get_perco(cmap,dim): - - if dim==2: - - pclusY=[] #list of the percolating clusters - for i in range(cmap.shape[0]): - if cmap[i,0] != 0: - if cmap[i,0] not in pclusY: - if cmap[i,0] in cmap[:,-1]: - pclusY+=[cmap[i,0]] - - - pclusZ=[] #list of the percolating clusters Z direction, this one is the main flow in Ndar.py, the fixed dimension is the direction used to see if pecolates - for i in range(cmap.shape[1]): - if cmap[0,i] != 0: - if cmap[0,i] not in pclusZ: - if cmap[0,i] in cmap[-1,:]: #viendo sin en la primer cara esta el mismo cluster que en la ultima - pclusZ+=[cmap[0,i]] - - - pclusX=[] - spanning=0 - if len(pclusZ)==1 and pclusZ==pclusY: - spanning=1 - - - if dim==3: - - - pclusX=[] #list of the percolating clusters - for i in range(cmap.shape[0]): # Z - for j in range(cmap.shape[1]): #X - if cmap[i,j,0] != 0: - if cmap[i,j,0] not in pclusX: - if cmap[i,j,0] in cmap[:,:,-1]: - pclusX+=[cmap[i,j,0]] - - pclusY=[] #list of the percolating clusters - for i in range(cmap.shape[0]): # Z - for k in range(cmap.shape[2]): #X - if cmap[i,0,k] != 0: - if cmap[i,0,k] not in pclusY: - if cmap[i,0,k] in cmap[:,-1,:]: - pclusY+=[cmap[i,0,k]] - - pclusZ=[] #list of the percolating clusters - for k in range(cmap.shape[2]): #x - for j in range(cmap.shape[1]): #y - if cmap[0,j,k] != 0: - if cmap[0,j,k] not in pclusZ: - if cmap[0,j,k] in cmap[-1,:,:]: - pclusZ+=[cmap[0,j,k]] #this is the one - - spanning=0 - if len(pclusZ)==1 and pclusZ==pclusY and pclusZ==pclusX: - spanning=1 - - - return spanning, pclusZ, pclusY, pclusX +def ConnecInd(cmap, scales, datadir): + + datadir = datadir + "ConnectivityMetrics/" + try: + os.makedirs(datadir) + except: + nada = 0 + + for scale in scales: + res = dict() + res = doforsubS_computeCmap(res, cmap, scale, postConec) + np.save(datadir + str(scale) + ".npy", res) + + return + + +def doforsubS_computeCmap(res, cmap, l, funpost): + + L = cmap.shape[0] + Nx, Ny, Nz = cmap.shape[0], cmap.shape[1], cmap.shape[2] + + nblx = Nx // l # for each dimension + nbly = Ny // l + if cmap.shape[2] == 1: + lz = 1 + nblz = 1 + else: + lz = l + nblz = Nz // l + + keys = funpost(np.array([]), res, 0, 0) + + for key in keys: + res[key] = np.zeros((nblx, nbly, nblz)) + + for i in range(nblx): + for j in range(nbly): + for k in range(nblz): + res = funpost( + cmap[ + i * l : (i + 1) * l, j * l : (j + 1) * l, k * l : (k + 1) * lz + ], + res, + (i, j, k), + 1, + ) + + return res + + +def postConec(cmap, results, ind, flag): + + if flag == 0: + keys = [] + keys += ["PPHA"] + keys += ["VOLALE"] + keys += ["ZNCC"] + keys += ["GAMMA"] + keys += ["spanning", "npz", "npy", "npx"] + keys += ["Plen", "S", "P"] + return keys + + dim = 3 + if cmap.shape[2] == 1: + cmap = cmap[:, :, 0] + dim = 2 + + y = np.bincount(cmap.reshape(-1)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + + if cf[0, 0] == 0: + cf = cf[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + + if cf.shape[0] > 0: + + spanning, pclusZ, pclusY, pclusX = get_perco(cmap, dim) + plen = Plen(spanning, cmap, cf, dim) + nper = np.sum(cf[:, 1]) # num de celdas permeables + nclus = cf.shape[0] # cantidad de clusters + + results["PPHA"][ind] = nper / np.size(cmap) # ppha + results["VOLALE"][ind] = ( + np.max(cf[:, 1]) / nper + ) # volale #corregido va entre [0,p] + results["ZNCC"][ind] = nclus # zncc + results["GAMMA"][ind] = ( + np.sum(cf[:, 1] ** 2) / np.size(cmap) / nper + ) # gamma, recordar zintcc =gamma*p + ( + results["spanning"][ind], + results["npz"][ind], + results["npy"][ind], + results["npx"][ind], + ) = (spanning, len(pclusZ), len(pclusY), len(pclusX)) + results["Plen"][ind], results["S"][ind], results["P"][ind] = ( + plen[0], + plen[1], + plen[2], + ) + + if cf.shape[0] == 0: + for key in keys: + results[key][ind] = 0 + return results + + +# ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA + + +def get_pos2D(cmap, cdis): + + Ns = cdis.shape[0] + pos = dict() + i = 0 + for cnum in cdis[:, 0]: + pos[cnum] = np.zeros((cdis[i, 1] + 1, 2)) # +1 porque uso de flag + i += 1 + + for i in range(cmap.shape[0]): + for j in range(cmap.shape[1]): + if cmap[i, j] != 0: + flag = int(pos[cmap[i, j]][0, 0]) + 1 + pos[cmap[i, j]][0, 0] = flag + pos[cmap[i, j]][flag, 0] = i + pos[cmap[i, j]][flag, 1] = j + return pos + + +def get_pos3D(cmap, cdis): + + Ns = cdis.shape[0] + pos = dict() + i = 0 + for cnum in cdis[:, 0]: + pos[cnum] = np.zeros((cdis[i, 1] + 1, 3)) + i += 1 + for i in range(cmap.shape[0]): + for j in range(cmap.shape[1]): + for k in range(cmap.shape[2]): + + if cmap[i, j, k] != 0: + flag = int(pos[cmap[i, j, k]][0, 0]) + 1 + pos[cmap[i, j, k]][0, 0] = flag + pos[cmap[i, j, k]][flag, 0] = i + pos[cmap[i, j, k]][flag, 1] = j + pos[cmap[i, j, k]][flag, 2] = k + + return pos + + +def Plen(spannng, cmap, cdis, dim): + + if dim == 2: + return P_len2D(spannng, cmap, cdis) + if dim == 3: + return P_len3D(spannng, cmap, cdis) + return [] + + +def P_len2D(spanning, cmap, cdis): + + pos = get_pos2D(cmap, cdis) + # print(summary['NpcY'],summary['NpcX'],summary['PPHA']) + + den = 0 + num = 0 + + nperm = np.sum(cdis[:, 1]) + if spanning > 0: + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + else: + P = 0 + + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy = np.mean(pos[cnum][1:, 0]), np.mean( + pos[cnum][1:, 1] + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + (pos[cnum][1:, 1] - mposy) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] + + +def P_len3D(spanning, cmap, cdis): + + pos = get_pos3D(cmap, cdis) + # print(summary['NpcY'],summary['NpcX'],summary['PPHA']) + + den = 0 + num = 0 + + nperm = np.sum(cdis[:, 1]) + if spanning > 0: + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + else: + P = 0 + + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy, mposz = ( + np.mean(pos[cnum][1:, 0]), + np.mean(pos[cnum][1:, 1]), + np.mean(pos[cnum][1:, 2]), + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + + (pos[cnum][1:, 1] - mposy) ** 2 + + (pos[cnum][1:, 2] - mposz) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] + + +def get_perco(cmap, dim): + + if dim == 2: + + pclusY = [] # list of the percolating clusters + for i in range(cmap.shape[0]): + if cmap[i, 0] != 0: + if cmap[i, 0] not in pclusY: + if cmap[i, 0] in cmap[:, -1]: + pclusY += [cmap[i, 0]] + + pclusZ = ( + [] + ) # list of the percolating clusters Z direction, this one is the main flow in Ndar.py, the fixed dimension is the direction used to see if pecolates + for i in range(cmap.shape[1]): + if cmap[0, i] != 0: + if cmap[0, i] not in pclusZ: + if ( + cmap[0, i] in cmap[-1, :] + ): # viendo sin en la primer cara esta el mismo cluster que en la ultima + pclusZ += [cmap[0, i]] + + pclusX = [] + spanning = 0 + if len(pclusZ) == 1 and pclusZ == pclusY: + spanning = 1 + + if dim == 3: + + pclusX = [] # list of the percolating clusters + for i in range(cmap.shape[0]): # Z + for j in range(cmap.shape[1]): # X + if cmap[i, j, 0] != 0: + if cmap[i, j, 0] not in pclusX: + if cmap[i, j, 0] in cmap[:, :, -1]: + pclusX += [cmap[i, j, 0]] + + pclusY = [] # list of the percolating clusters + for i in range(cmap.shape[0]): # Z + for k in range(cmap.shape[2]): # X + if cmap[i, 0, k] != 0: + if cmap[i, 0, k] not in pclusY: + if cmap[i, 0, k] in cmap[:, -1, :]: + pclusY += [cmap[i, 0, k]] + + pclusZ = [] # list of the percolating clusters + for k in range(cmap.shape[2]): # x + for j in range(cmap.shape[1]): # y + if cmap[0, j, k] != 0: + if cmap[0, j, k] not in pclusZ: + if cmap[0, j, k] in cmap[-1, :, :]: + pclusZ += [cmap[0, j, k]] # this is the one + + spanning = 0 + if len(pclusZ) == 1 and pclusZ == pclusY and pclusZ == pclusX: + spanning = 1 + + return spanning, pclusZ, pclusY, pclusX diff --git a/tools/connec/PostConec.py b/tools/connec/PostConec.py index 02ea2fa..e5e31a1 100755 --- a/tools/connec/PostConec.py +++ b/tools/connec/PostConec.py @@ -5,229 +5,271 @@ import os import collections - -def ConnecInd(cmap,scales,datadir): - - datadir=datadir+'ConnectivityMetrics/' - try: - os.makedirs(datadir) - except: - nada=0 - - for scale in scales: - res=dict() - res=doforsubS_computeCmap(res,cmap,scale,postConec) - np.save(datadir+str(scale)+'.npy',res) - - return - -def doforsubS_computeCmap(res,cmap,l,funpost): - - - L=cmap.shape[0] - Nx, Ny,Nz=cmap.shape[0],cmap.shape[1],cmap.shape[2] - - nblx=Nx//l #for each dimension - - ly=l - nbly=Ny//l - - lz=l - nblz=Nz//l - - if nbly==0: #si l> Ny - nbly=1 - ly=Ny - - if nblz==0: - lz=1 - nblz=1 - - - keys=funpost(np.array([]),res,0,0) - - for key in keys: - res[key]=np.zeros((nblx,nbly,nblz)) - - for i in range(nblx): - for j in range(nbly): - for k in range(nblz): - res=funpost(cmap[i*l:(i+1)*l,j*ly:(j+1)*ly,k*l:(k+1)*lz],res,(i,j,k),1) - - return res - -def postConec(cmap,results,ind,flag): - - keys=[] - keys+=['PPHA'] - keys+=['VOLALE'] - keys+=['ZNCC'] - keys+=['GAMMA'] - keys+=['spanning', 'npz', 'npy', 'npx'] - keys+=['Plen','S','P'] - keys+=['PlenX','SX','PX'] - if flag==0: - - return keys - - - y = np.bincount(cmap.reshape(-1)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - - if cf[0,0]==0: - cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - - if cf.shape[0]>0: - spanning, pclusX, pclusY, pclusZ =get_perco(cmap) - plen=Plen(spanning,cmap,cf) - #print(pclusX,spanning) - if len(pclusX) >0 and spanning ==0: - plenX=PlenX(pclusX,cmap,cf) - else: - plenX=plen - nper=np.sum(cf[:,1]) #num de celdas permeables - nclus=cf.shape[0] #cantidad de clusters - results['PPHA'][ind]=nper/np.size(cmap) #ppha - results['VOLALE'][ind]=np.max(cf[:,1])/nper #volale #corregido va entre [0,p] - results['ZNCC'][ind]=nclus #zncc - results['GAMMA'][ind]=np.sum(cf[:,1]**2)/nper**2 #gamma, recordar zintcc =gamma*nper - results['spanning'][ind],results['npz'][ind], results['npy'][ind], results['npx'][ind]=spanning, len(pclusZ), len(pclusY), len(pclusX) - results['Plen'][ind],results['S'][ind],results['P'][ind] = plen[0],plen[1],plen[2] - results['PlenX'][ind],results['SX'][ind],results['PX'][ind] = plenX[0],plenX[1],plenX[2] - if cf.shape[0]==0: - for key in keys: - results[key][ind]=0 - return results - -def get_pos(cmap,cdis): - - Ns=cdis.shape[0] - pos=dict() - i=0 - for cnum in cdis[:,0]: - pos[cnum]=np.zeros((cdis[i,1]+1,3)) - i+=1 - for i in range(cmap.shape[0]): - for j in range(cmap.shape[1]): - for k in range(cmap.shape[2]): - - if cmap[i,j,k] != 0: - flag=int(pos[cmap[i,j,k]][0,0])+1 - pos[cmap[i,j,k]][0,0]=flag - pos[cmap[i,j,k]][flag,0]=i - pos[cmap[i,j,k]][flag,1]=j - pos[cmap[i,j,k]][flag,2]=k - - return pos - -def Plen(spanning,cmap,cdis): - - pos = get_pos(cmap,cdis) - #print(summary['NpcY'],summary['NpcX'],summary['PPHA']) - - den=0 - num=0 - - nperm=np.sum(cdis[:,1]) - if spanning > 0: - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - else: - P=0 - - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy, mposz = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]), np.mean(pos[cnum][1:,2]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2+(pos[cnum][1:,2]-mposz)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - - -def PlenX(pclusX,cmap,cdis): - - #guarda que solo se entra en esta funcion si no es spanning pero hay al menos 1 cluster percolante en X - - for cluster in pclusX[1:]: - cmap=np.where(cmap==cluster,pclusX[0],cmap) - - - y = np.bincount(cmap.reshape(-1)) - ii = np.nonzero(y)[0] - cdis=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - - if cdis[0,0]==0: - cdis=cdis[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - - - pos = get_pos(cmap,cdis) - nperm=np.sum(cdis[:,1]) - - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - den=0 - num=0 - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy, mposz = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]), np.mean(pos[cnum][1:,2]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2+(pos[cnum][1:,2]-mposz)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - +def ConnecInd(cmap, scales, datadir): + + datadir = datadir + "ConnectivityMetrics/" + try: + os.makedirs(datadir) + except: + nada = 0 + + for scale in scales: + res = dict() + res = doforsubS_computeCmap(res, cmap, scale, postConec) + np.save(datadir + str(scale) + ".npy", res) + + return + + +def doforsubS_computeCmap(res, cmap, l, funpost): + + L = cmap.shape[0] + Nx, Ny, Nz = cmap.shape[0], cmap.shape[1], cmap.shape[2] + + nblx = Nx // l # for each dimension + + ly = l + nbly = Ny // l + + lz = l + nblz = Nz // l + + if nbly == 0: # si l> Ny + nbly = 1 + ly = Ny + + if nblz == 0: + lz = 1 + nblz = 1 + + keys = funpost(np.array([]), res, 0, 0) + + for key in keys: + res[key] = np.zeros((nblx, nbly, nblz)) + + for i in range(nblx): + for j in range(nbly): + for k in range(nblz): + res = funpost( + cmap[ + i * l : (i + 1) * l, j * ly : (j + 1) * ly, k * l : (k + 1) * lz + ], + res, + (i, j, k), + 1, + ) + + return res + + +def postConec(cmap, results, ind, flag): + + keys = [] + keys += ["PPHA"] + keys += ["VOLALE"] + keys += ["ZNCC"] + keys += ["GAMMA"] + keys += ["spanning", "npz", "npy", "npx"] + keys += ["Plen", "S", "P"] + keys += ["PlenX", "SX", "PX"] + if flag == 0: + + return keys + + y = np.bincount(cmap.reshape(-1)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + + if cf[0, 0] == 0: + cf = cf[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + + if cf.shape[0] > 0: + spanning, pclusX, pclusY, pclusZ = get_perco(cmap) + plen = Plen(spanning, cmap, cf) + # print(pclusX,spanning) + if len(pclusX) > 0 and spanning == 0: + plenX = PlenX(pclusX, cmap, cf) + else: + plenX = plen + nper = np.sum(cf[:, 1]) # num de celdas permeables + nclus = cf.shape[0] # cantidad de clusters + results["PPHA"][ind] = nper / np.size(cmap) # ppha + results["VOLALE"][ind] = ( + np.max(cf[:, 1]) / nper + ) # volale #corregido va entre [0,p] + results["ZNCC"][ind] = nclus # zncc + results["GAMMA"][ind] = ( + np.sum(cf[:, 1] ** 2) / nper ** 2 + ) # gamma, recordar zintcc =gamma*nper + ( + results["spanning"][ind], + results["npz"][ind], + results["npy"][ind], + results["npx"][ind], + ) = (spanning, len(pclusZ), len(pclusY), len(pclusX)) + results["Plen"][ind], results["S"][ind], results["P"][ind] = ( + plen[0], + plen[1], + plen[2], + ) + results["PlenX"][ind], results["SX"][ind], results["PX"][ind] = ( + plenX[0], + plenX[1], + plenX[2], + ) + if cf.shape[0] == 0: + for key in keys: + results[key][ind] = 0 + return results + + +def get_pos(cmap, cdis): + + Ns = cdis.shape[0] + pos = dict() + i = 0 + for cnum in cdis[:, 0]: + pos[cnum] = np.zeros((cdis[i, 1] + 1, 3)) + i += 1 + for i in range(cmap.shape[0]): + for j in range(cmap.shape[1]): + for k in range(cmap.shape[2]): + + if cmap[i, j, k] != 0: + flag = int(pos[cmap[i, j, k]][0, 0]) + 1 + pos[cmap[i, j, k]][0, 0] = flag + pos[cmap[i, j, k]][flag, 0] = i + pos[cmap[i, j, k]][flag, 1] = j + pos[cmap[i, j, k]][flag, 2] = k + + return pos + + +def Plen(spanning, cmap, cdis): + + pos = get_pos(cmap, cdis) + # print(summary['NpcY'],summary['NpcX'],summary['PPHA']) + + den = 0 + num = 0 + + nperm = np.sum(cdis[:, 1]) + if spanning > 0: + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + else: + P = 0 + + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy, mposz = ( + np.mean(pos[cnum][1:, 0]), + np.mean(pos[cnum][1:, 1]), + np.mean(pos[cnum][1:, 2]), + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + + (pos[cnum][1:, 1] - mposy) ** 2 + + (pos[cnum][1:, 2] - mposz) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] + + +def PlenX(pclusX, cmap, cdis): + + # guarda que solo se entra en esta funcion si no es spanning pero hay al menos 1 cluster percolante en X + + for cluster in pclusX[1:]: + cmap = np.where(cmap == cluster, pclusX[0], cmap) + + y = np.bincount(cmap.reshape(-1)) + ii = np.nonzero(y)[0] + cdis = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + + if cdis[0, 0] == 0: + cdis = cdis[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + + pos = get_pos(cmap, cdis) + nperm = np.sum(cdis[:, 1]) + + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + den = 0 + num = 0 + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy, mposz = ( + np.mean(pos[cnum][1:, 0]), + np.mean(pos[cnum][1:, 1]), + np.mean(pos[cnum][1:, 2]), + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + + (pos[cnum][1:, 1] - mposy) ** 2 + + (pos[cnum][1:, 2] - mposz) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] def get_perco(cmap): - - pclusX=[] #list of the percolating clusters - for k in range(cmap.shape[2]): #x - for j in range(cmap.shape[1]): #y - if cmap[0,j,k] != 0: - if cmap[0,j,k] not in pclusX: - if cmap[0,j,k] in cmap[-1,:,:]: - pclusX+=[cmap[0,j,k]] #this is the one - - pclusY=[] #list of the percolating clusters - for i in range(cmap.shape[0]): # Z - for k in range(cmap.shape[2]): #X - if cmap[i,0,k] != 0: - if cmap[i,0,k] not in pclusY: - if cmap[i,0,k] in cmap[:,-1,:]: - pclusY+=[cmap[i,0,k]] - - pclusZ=[] #list of the percolating clusters - if cmap.shape[2]>1: - for i in range(cmap.shape[0]): # Z - for j in range(cmap.shape[1]): #X - if cmap[i,j,0] != 0: - if cmap[i,j,0] not in pclusZ: - if cmap[i,j,0] in cmap[:,:,-1]: - pclusZ+=[cmap[i,j,0]] - - spanning=0 - if len(pclusZ)==1 and pclusZ==pclusY and pclusZ==pclusX: - spanning=1 - else: - spanning=0 - if len(pclusX)==1 and pclusY==pclusX: - spanning=1 - - return spanning, pclusX, pclusY, pclusZ + pclusX = [] # list of the percolating clusters + for k in range(cmap.shape[2]): # x + for j in range(cmap.shape[1]): # y + if cmap[0, j, k] != 0: + if cmap[0, j, k] not in pclusX: + if cmap[0, j, k] in cmap[-1, :, :]: + pclusX += [cmap[0, j, k]] # this is the one + + pclusY = [] # list of the percolating clusters + for i in range(cmap.shape[0]): # Z + for k in range(cmap.shape[2]): # X + if cmap[i, 0, k] != 0: + if cmap[i, 0, k] not in pclusY: + if cmap[i, 0, k] in cmap[:, -1, :]: + pclusY += [cmap[i, 0, k]] + + pclusZ = [] # list of the percolating clusters + if cmap.shape[2] > 1: + for i in range(cmap.shape[0]): # Z + for j in range(cmap.shape[1]): # X + if cmap[i, j, 0] != 0: + if cmap[i, j, 0] not in pclusZ: + if cmap[i, j, 0] in cmap[:, :, -1]: + pclusZ += [cmap[i, j, 0]] + + spanning = 0 + if len(pclusZ) == 1 and pclusZ == pclusY and pclusZ == pclusX: + spanning = 1 + else: + spanning = 0 + if len(pclusX) == 1 and pclusY == pclusX: + spanning = 1 + + return spanning, pclusX, pclusY, pclusZ diff --git a/tools/connec/PostConec_v2.py b/tools/connec/PostConec_v2.py index e4cd7aa..90962d7 100755 --- a/tools/connec/PostConec_v2.py +++ b/tools/connec/PostConec_v2.py @@ -6,363 +6,395 @@ import os import collections - def main(): - - #scales=[4,6,8,16,24,32] - #numofseeds=np.array([10,10,10,48,100,200]) - #startseed=1 - - scales=[2,4,8,12,16,20,26,32] - numofseeds=np.array([1,2,12,16,20,25,30,50]) - - startseed=1 - dim=3 - - numofseeds=numofseeds+startseed - - mapa=np.loadtxt(('vecconec.txt')).astype(int) - - if dim==2: - LL=int(np.sqrt(mapa.shape[0])) - mapa=mapa.reshape(LL,LL) - - if dim==3: - LL=int(np.cbrt(mapa.shape[0])) - mapa=mapa.reshape(LL,LL,LL) - res, names=doforsubS_computeCmap(mapa,scales,postConec, compCon,dim,[],numofseeds) - - with open('keysCon.txt', 'w') as f: - for item in names: - f.write("%s\n" % item) - - np.save('ConResScales.npy',res) - - - return - -def doforsubS_computeCmap(mapa,scales,funpost, funcompCmap,dim,args,numofseeds): - - L=mapa.shape[0] - res=dict() - names=[] - - - with open('Kfield.don') as f: - seed = int(f.readline()) - - for iscale in range(len(scales)): - - l=scales[iscale] - - if numofseeds[iscale] > seed: #guarda aca - - nblocks=L//l #for each dimension - - if dim==2: - for i in range(nblocks): - for j in range(nblocks): - cmapa=funcompCmap(mapa[i*l:(i+1)*l,j*l:(j+1)*l],dim) - dats,names=funpost(cmapa,dim,args) - if i== 0 and j==0: - for icon in range(len(names)): - res[l,names[icon]]=[] - for icon in range(len(names)): - res[l,names[icon]]+=[dats[icon]] - - - if dim==3: - for i in range(nblocks): - for j in range(nblocks): - for k in range(nblocks): - cmapa=funcompCmap(mapa[i*l:(i+1)*l,j*l:(j+1)*l,k*l:(k+1)*l],dim) - dats, names=funpost(cmapa,dim,args) - if i== 0 and j==0 and k==0: - for icon in range(len(names)): - res[l,names[icon]]=[] - for icon in range(len(names)): - res[l,names[icon]]+=[dats[icon]] - - - - return res, names - - -def ConConfig(L,dim): - - params=[] - if dim==2: - params=['1','4','imap.txt',str(L)+' '+str(L),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec2d' - - if dim==3: - params=['1','6','imap.txt',str(L)+' '+str(L)+' ' +str(L),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec3d' - return params, execCon - - - -def compCon(mapa,dim): - - exeDir='./' - L=mapa.shape[0] - params,execCon=ConConfig(L,dim) - - with open(exeDir+'coninput.txt', 'w') as f: - for item in params: - f.write("%s\n" % item) - np.savetxt(exeDir+params[2],mapa.reshape(-1)) - - #wiam=os.getcwd() - #os.chdir(exeDir) - os.system('cp ../../../tools/conec3d ./') - os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir+ - - cmapa=np.loadtxt(params[-2]).reshape(mapa.shape).astype(int) #exeDir+ - #os.chdir(wiam) - return cmapa - - - - - -def postConec(cmap,dim,args): - - names=['PPHA','VOLALE','ZNCC','zintcc','spaninning','npz','npy','npx',] - - - L=cmap.shape[0] - results=[] - names=[] - - - y = np.bincount(cmap.reshape(-1)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - - if cf[0,0]==0: - cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - - if cf.shape[0]>0: - # headers=['N','p','Csize','CLenX','CLenY','CmaxVol','MaxLenX','MaxLenY','NpcX','NpcY'] - - nper=np.sum(cf[:,1]) #num de celdas permeables - nclus=cf.shape[0] #cantidad de clusters - - - #ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA - results+=[nper/np.size(cmap)] #ppha - results+=[np.max(cf[:,1])/nper] #volale #corregido va entre [0,p] - results+=[nclus] #zncc - results+=[np.sum(cf[:,1]**2)/np.size(cmap)/nper] #gamma, recordar zintcc =gamma*p - - spanning, pclusZ, pclusY, pclusX =get_perco(cmap,dim) - results+=[spanning, len(pclusZ), len(pclusY), len(pclusX)] - - - results+=Plen(spanning,cmap,cf,dim) - - - names+=['PPHA'] - names+=['VOLALE'] - names+=['ZNCC'] - names+=['ZINTCC'] - names+=['spanning', 'npz', 'npy', 'npx'] - names+=['Plen','S','P'] - - if cf.shape[0]==0: - for i in range(len(names)): - results+=[0] - return results, names - - -#ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA - - -def get_pos2D(cmap,cdis): - - Ns=cdis.shape[0] - pos=dict() - i=0 - for cnum in cdis[:,0]: - pos[cnum]=np.zeros((cdis[i,1]+1,2)) #+1 porque uso de flag - i+=1 - - for i in range(cmap.shape[0]): - for j in range(cmap.shape[1]): - if cmap[i,j] != 0: - flag=int(pos[cmap[i,j]][0,0])+1 - pos[cmap[i,j]][0,0]=flag - pos[cmap[i,j]][flag,0]=i - pos[cmap[i,j]][flag,1]=j - - - - return pos - - -def get_pos3D(cmap,cdis): - - Ns=cdis.shape[0] - pos=dict() - i=0 - for cnum in cdis[:,0]: - pos[cnum]=np.zeros((cdis[i,1]+1,3)) - i+=1 - for i in range(cmap.shape[0]): - for j in range(cmap.shape[1]): - for k in range(cmap.shape[2]): - - if cmap[i,j,k] != 0: - flag=int(pos[cmap[i,j,k]][0,0])+1 - pos[cmap[i,j,k]][0,0]=flag - pos[cmap[i,j,k]][flag,0]=i - pos[cmap[i,j,k]][flag,1]=j - pos[cmap[i,j,k]][flag,2]=k - - - return pos - -def Plen(spannng,cmap,cdis,dim): - - if dim==2: - return P_len2D(spannng,cmap,cdis) - if dim==3: - return P_len3D(spannng,cmap,cdis) - return [] - -def P_len2D(spanning,cmap,cdis): - - - pos = get_pos2D(cmap,cdis) - #print(summary['NpcY'],summary['NpcX'],summary['PPHA']) - - den=0 - num=0 - - nperm=np.sum(cdis[:,1]) - if spanning > 0: - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - else: - P=0 - - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - - - -def P_len3D(spanning,cmap,cdis): - - - pos = get_pos3D(cmap,cdis) - #print(summary['NpcY'],summary['NpcX'],summary['PPHA']) - - den=0 - num=0 - - nperm=np.sum(cdis[:,1]) - if spanning > 0: - amax=np.argmax(cdis[:,1]) - P=cdis[amax,1]/nperm - cdis=np.delete(cdis,amax,axis=0) - - else: - P=0 - - i=0 - if cdis.shape[0]> 0: - S=np.sum(cdis[:,1])/(cdis.shape[0]) - for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k- - mposx, mposy, mposz = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]), np.mean(pos[cnum][1:,2]) #el 1: de sacar el flag - Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2+(pos[cnum][1:,2]-mposz)**2) #Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik - num += cdis[i,1]**2 * Rs - den+=cdis[i,1]**2 - i+=1 - return [np.sqrt(num/den), S, P] - else: - return [0,0,P] - - - - -def get_perco(cmap,dim): - - if dim==2: - - pclusY=[] #list of the percolating clusters - for i in range(cmap.shape[0]): - if cmap[i,0] != 0: - if cmap[i,0] not in pclusY: - if cmap[i,0] in cmap[:,-1]: - pclusY+=[cmap[i,0]] - - - pclusZ=[] #list of the percolating clusters Z direction, this one is the main flow in Ndar.py, the fixed dimension is the direction used to see if pecolates - for i in range(cmap.shape[1]): - if cmap[0,i] != 0: - if cmap[0,i] not in pclusZ: - if cmap[0,i] in cmap[-1,:]: #viendo sin en la primer cara esta el mismo cluster que en la ultima - pclusZ+=[cmap[0,i]] - - - pclusX=[] - spanning=0 - if len(pclusZ)==1 and pclusZ==pclusY: - spanning=1 - - - if dim==3: - - - pclusX=[] #list of the percolating clusters - for i in range(cmap.shape[0]): # Z - for j in range(cmap.shape[1]): #X - if cmap[i,j,0] != 0: - if cmap[i,j,0] not in pclusX: - if cmap[i,j,0] in cmap[:,:,-1]: - pclusX+=[cmap[i,j,0]] - - pclusY=[] #list of the percolating clusters - for i in range(cmap.shape[0]): # Z - for k in range(cmap.shape[2]): #X - if cmap[i,0,k] != 0: - if cmap[i,0,k] not in pclusY: - if cmap[i,0,k] in cmap[:,-1,:]: - pclusY+=[cmap[i,0,k]] - - pclusZ=[] #list of the percolating clusters - for k in range(cmap.shape[2]): #x - for j in range(cmap.shape[1]): #y - if cmap[0,j,k] != 0: - if cmap[0,j,k] not in pclusZ: - if cmap[0,j,k] in cmap[-1,:,:]: - pclusZ+=[cmap[0,j,k]] #this is the one - - spanning=0 - if len(pclusZ)==1 and pclusZ==pclusY and pclusZ==pclusX: - spanning=1 - - - return spanning, pclusZ, pclusY, pclusX + # scales=[4,6,8,16,24,32] + # numofseeds=np.array([10,10,10,48,100,200]) + # startseed=1 + + scales = [2, 4, 8, 12, 16, 20, 26, 32] + numofseeds = np.array([1, 2, 12, 16, 20, 25, 30, 50]) + + startseed = 1 + dim = 3 + + numofseeds = numofseeds + startseed + + mapa = np.loadtxt(("vecconec.txt")).astype(int) + + if dim == 2: + LL = int(np.sqrt(mapa.shape[0])) + mapa = mapa.reshape(LL, LL) + + if dim == 3: + LL = int(np.cbrt(mapa.shape[0])) + mapa = mapa.reshape(LL, LL, LL) + res, names = doforsubS_computeCmap( + mapa, scales, postConec, compCon, dim, [], numofseeds + ) + + with open("keysCon.txt", "w") as f: + for item in names: + f.write("%s\n" % item) + + np.save("ConResScales.npy", res) + + return + + +def doforsubS_computeCmap(mapa, scales, funpost, funcompCmap, dim, args, numofseeds): + + L = mapa.shape[0] + res = dict() + names = [] + + with open("Kfield.don") as f: + seed = int(f.readline()) + + for iscale in range(len(scales)): + + l = scales[iscale] + + if numofseeds[iscale] > seed: # guarda aca + + nblocks = L // l # for each dimension + + if dim == 2: + for i in range(nblocks): + for j in range(nblocks): + cmapa = funcompCmap( + mapa[i * l : (i + 1) * l, j * l : (j + 1) * l], dim + ) + dats, names = funpost(cmapa, dim, args) + if i == 0 and j == 0: + for icon in range(len(names)): + res[l, names[icon]] = [] + for icon in range(len(names)): + res[l, names[icon]] += [dats[icon]] + + if dim == 3: + for i in range(nblocks): + for j in range(nblocks): + for k in range(nblocks): + cmapa = funcompCmap( + mapa[ + i * l : (i + 1) * l, + j * l : (j + 1) * l, + k * l : (k + 1) * l, + ], + dim, + ) + dats, names = funpost(cmapa, dim, args) + if i == 0 and j == 0 and k == 0: + for icon in range(len(names)): + res[l, names[icon]] = [] + for icon in range(len(names)): + res[l, names[icon]] += [dats[icon]] + + return res, names + + +def ConConfig(L, dim): + + params = [] + if dim == 2: + params = [ + "1", + "4", + "imap.txt", + str(L) + " " + str(L), + "1.0 1.0", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec2d" + + if dim == 3: + params = [ + "1", + "6", + "imap.txt", + str(L) + " " + str(L) + " " + str(L), + "1.0 1.0 1.0", + "30", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec3d" + return params, execCon + + +def compCon(mapa, dim): + + exeDir = "./" + L = mapa.shape[0] + params, execCon = ConConfig(L, dim) + + with open(exeDir + "coninput.txt", "w") as f: + for item in params: + f.write("%s\n" % item) + np.savetxt(exeDir + params[2], mapa.reshape(-1)) + + # wiam=os.getcwd() + # os.chdir(exeDir) + os.system("cp ../../../tools/conec3d ./") + os.system(" ./" + execCon + ">/dev/null") #'cd ' +exeDir+ + + cmapa = np.loadtxt(params[-2]).reshape(mapa.shape).astype(int) # exeDir+ + # os.chdir(wiam) + return cmapa + + +def postConec(cmap, dim, args): + + names = [ + "PPHA", + "VOLALE", + "ZNCC", + "zintcc", + "spaninning", + "npz", + "npy", + "npx", + ] + + L = cmap.shape[0] + results = [] + names = [] + + y = np.bincount(cmap.reshape(-1)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + + if cf[0, 0] == 0: + cf = cf[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + + if cf.shape[0] > 0: + # headers=['N','p','Csize','CLenX','CLenY','CmaxVol','MaxLenX','MaxLenY','NpcX','NpcY'] + + nper = np.sum(cf[:, 1]) # num de celdas permeables + nclus = cf.shape[0] # cantidad de clusters + + # ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA + results += [nper / np.size(cmap)] # ppha + results += [np.max(cf[:, 1]) / nper] # volale #corregido va entre [0,p] + results += [nclus] # zncc + results += [ + np.sum(cf[:, 1] ** 2) / np.size(cmap) / nper + ] # gamma, recordar zintcc =gamma*p + + spanning, pclusZ, pclusY, pclusX = get_perco(cmap, dim) + results += [spanning, len(pclusZ), len(pclusY), len(pclusX)] + + results += Plen(spanning, cmap, cf, dim) + + names += ["PPHA"] + names += ["VOLALE"] + names += ["ZNCC"] + names += ["ZINTCC"] + names += ["spanning", "npz", "npy", "npx"] + names += ["Plen", "S", "P"] + + if cf.shape[0] == 0: + for i in range(len(names)): + results += [0] + return results, names + + +# ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA + + +def get_pos2D(cmap, cdis): + + Ns = cdis.shape[0] + pos = dict() + i = 0 + for cnum in cdis[:, 0]: + pos[cnum] = np.zeros((cdis[i, 1] + 1, 2)) # +1 porque uso de flag + i += 1 + + for i in range(cmap.shape[0]): + for j in range(cmap.shape[1]): + if cmap[i, j] != 0: + flag = int(pos[cmap[i, j]][0, 0]) + 1 + pos[cmap[i, j]][0, 0] = flag + pos[cmap[i, j]][flag, 0] = i + pos[cmap[i, j]][flag, 1] = j + + return pos + + +def get_pos3D(cmap, cdis): + + Ns = cdis.shape[0] + pos = dict() + i = 0 + for cnum in cdis[:, 0]: + pos[cnum] = np.zeros((cdis[i, 1] + 1, 3)) + i += 1 + for i in range(cmap.shape[0]): + for j in range(cmap.shape[1]): + for k in range(cmap.shape[2]): + + if cmap[i, j, k] != 0: + flag = int(pos[cmap[i, j, k]][0, 0]) + 1 + pos[cmap[i, j, k]][0, 0] = flag + pos[cmap[i, j, k]][flag, 0] = i + pos[cmap[i, j, k]][flag, 1] = j + pos[cmap[i, j, k]][flag, 2] = k + + return pos + + +def Plen(spannng, cmap, cdis, dim): + + if dim == 2: + return P_len2D(spannng, cmap, cdis) + if dim == 3: + return P_len3D(spannng, cmap, cdis) + return [] + + +def P_len2D(spanning, cmap, cdis): + + pos = get_pos2D(cmap, cdis) + # print(summary['NpcY'],summary['NpcX'],summary['PPHA']) + + den = 0 + num = 0 + + nperm = np.sum(cdis[:, 1]) + if spanning > 0: + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + else: + P = 0 + + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy = np.mean(pos[cnum][1:, 0]), np.mean( + pos[cnum][1:, 1] + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + (pos[cnum][1:, 1] - mposy) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] + + +def P_len3D(spanning, cmap, cdis): + + pos = get_pos3D(cmap, cdis) + # print(summary['NpcY'],summary['NpcX'],summary['PPHA']) + + den = 0 + num = 0 + + nperm = np.sum(cdis[:, 1]) + if spanning > 0: + amax = np.argmax(cdis[:, 1]) + P = cdis[amax, 1] / nperm + cdis = np.delete(cdis, amax, axis=0) + + else: + P = 0 + + i = 0 + if cdis.shape[0] > 0: + S = np.sum(cdis[:, 1]) / (cdis.shape[0]) + for cnum in cdis[ + :, 0 + ]: # los clusters estan numerados a partir de 1, cluster cero es k- + mposx, mposy, mposz = ( + np.mean(pos[cnum][1:, 0]), + np.mean(pos[cnum][1:, 1]), + np.mean(pos[cnum][1:, 2]), + ) # el 1: de sacar el flag + Rs = np.mean( + (pos[cnum][1:, 0] - mposx) ** 2 + + (pos[cnum][1:, 1] - mposy) ** 2 + + (pos[cnum][1:, 2] - mposz) ** 2 + ) # Rs cuadrado ecuacion 12.9 libro Harvey Gould, Jan Tobochnik + num += cdis[i, 1] ** 2 * Rs + den += cdis[i, 1] ** 2 + i += 1 + return [np.sqrt(num / den), S, P] + else: + return [0, 0, P] + + +def get_perco(cmap, dim): + + if dim == 2: + + pclusY = [] # list of the percolating clusters + for i in range(cmap.shape[0]): + if cmap[i, 0] != 0: + if cmap[i, 0] not in pclusY: + if cmap[i, 0] in cmap[:, -1]: + pclusY += [cmap[i, 0]] + + pclusZ = ( + [] + ) # list of the percolating clusters Z direction, this one is the main flow in Ndar.py, the fixed dimension is the direction used to see if pecolates + for i in range(cmap.shape[1]): + if cmap[0, i] != 0: + if cmap[0, i] not in pclusZ: + if ( + cmap[0, i] in cmap[-1, :] + ): # viendo sin en la primer cara esta el mismo cluster que en la ultima + pclusZ += [cmap[0, i]] + + pclusX = [] + spanning = 0 + if len(pclusZ) == 1 and pclusZ == pclusY: + spanning = 1 + + if dim == 3: + + pclusX = [] # list of the percolating clusters + for i in range(cmap.shape[0]): # Z + for j in range(cmap.shape[1]): # X + if cmap[i, j, 0] != 0: + if cmap[i, j, 0] not in pclusX: + if cmap[i, j, 0] in cmap[:, :, -1]: + pclusX += [cmap[i, j, 0]] + + pclusY = [] # list of the percolating clusters + for i in range(cmap.shape[0]): # Z + for k in range(cmap.shape[2]): # X + if cmap[i, 0, k] != 0: + if cmap[i, 0, k] not in pclusY: + if cmap[i, 0, k] in cmap[:, -1, :]: + pclusY += [cmap[i, 0, k]] + + pclusZ = [] # list of the percolating clusters + for k in range(cmap.shape[2]): # x + for j in range(cmap.shape[1]): # y + if cmap[0, j, k] != 0: + if cmap[0, j, k] not in pclusZ: + if cmap[0, j, k] in cmap[-1, :, :]: + pclusZ += [cmap[0, j, k]] # this is the one + + spanning = 0 + if len(pclusZ) == 1 and pclusZ == pclusY and pclusZ == pclusX: + spanning = 1 + + return spanning, pclusZ, pclusY, pclusX main() - - - - diff --git a/tools/connec/back/comp_cmap_scales.py b/tools/connec/back/comp_cmap_scales.py index 6da035b..4cfbf21 100755 --- a/tools/connec/back/comp_cmap_scales.py +++ b/tools/connec/back/comp_cmap_scales.py @@ -2,129 +2,182 @@ import numpy as np import os import time from JoinCmaps import * -#k[x,y,z] - -def div_veccon(kc,kh,nbl,rundir): - - t0=time.time() - kc=np.where(kc==kh,1,0).astype(int) - - tcmaps=time.time() - kc=get_smallCmap(kc,nbl,rundir) - tcmaps=time.time()-tcmaps - #if s_scale/dev/null') #'cd ' +exeDir++'>/dev/null' - vec=np.loadtxt(params[-2]).reshape(vec.shape[0],vec.shape[1],vec.shape[2]).astype(int) - return vec - -def ConConfig(sx,sy,sz,Nz,rundir): - - params=[] - if Nz==1: - params=['1','4','vecconec.txt',str(sx)+' '+str(sy),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec2d' - - else: - params=['1','6','vecconec.txt',str(sx)+' '+str(sy)+' ' +str(sz),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec3d' - - - with open(rundir+'coninput.txt', 'w') as f: - for item in params: - f.write("%s\n" % item) - - return params, execCon - -def join(vec,nbl): - - - Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2] - sx,sy,sz = Nx//nbl,Ny//nbl,Nz//nbl - ex,ey,ez=np.log2(Nx),np.log2(Ny),np.log2(Nz) - - if Nz==1: - sz=1 - nbz=1 - ez=1 - esz=1 - else: - esz=np.log2(sz) - - - esx,esy=np.log2(sx),np.log2(sy) - - - - for bs in range(0,int(ex-esx)): - - nbx,nby = int(2**(ex-esx-bs-1)),int(2**(ey-esy-bs-1)) - if Nz==1: - sz=1 - nbz=1 - else: - nbz=int(2**(ez-esz-bs-1)) - sz=Nz//nbz - sx,sy=Nx//nbx,Ny//nby - - for i in range(nbx): - for j in range(nby): - for k in range(nbz): - a=2 - vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=joinBox(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],True,False) - - return vec - -''' +# k[x,y,z] + + +def div_veccon(kc, kh, nbl, rundir): + + t0 = time.time() + kc = np.where(kc == kh, 1, 0).astype(int) + + tcmaps = time.time() + kc = get_smallCmap(kc, nbl, rundir) + tcmaps = time.time() - tcmaps + # if s_scale/dev/null") #'cd ' +exeDir++'>/dev/null' + vec = ( + np.loadtxt(params[-2]) + .reshape(vec.shape[0], vec.shape[1], vec.shape[2]) + .astype(int) + ) + return vec + + +def ConConfig(sx, sy, sz, Nz, rundir): + + params = [] + if Nz == 1: + params = [ + "1", + "4", + "vecconec.txt", + str(sx) + " " + str(sy), + "1.0 1.0", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec2d" + + else: + params = [ + "1", + "6", + "vecconec.txt", + str(sx) + " " + str(sy) + " " + str(sz), + "1.0 1.0 1.0", + "30", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec3d" + + with open(rundir + "coninput.txt", "w") as f: + for item in params: + f.write("%s\n" % item) + + return params, execCon + + +def join(vec, nbl): + + Nx, Ny, Nz = vec.shape[0], vec.shape[1], vec.shape[2] + sx, sy, sz = Nx // nbl, Ny // nbl, Nz // nbl + ex, ey, ez = np.log2(Nx), np.log2(Ny), np.log2(Nz) + + if Nz == 1: + sz = 1 + nbz = 1 + ez = 1 + esz = 1 + else: + esz = np.log2(sz) + + esx, esy = np.log2(sx), np.log2(sy) + + for bs in range(0, int(ex - esx)): + + nbx, nby = int(2 ** (ex - esx - bs - 1)), int(2 ** (ey - esy - bs - 1)) + if Nz == 1: + sz = 1 + nbz = 1 + else: + nbz = int(2 ** (ez - esz - bs - 1)) + sz = Nz // nbz + sx, sy = Nx // nbx, Ny // nby + + for i in range(nbx): + for j in range(nby): + for k in range(nbz): + a = 2 + vec[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ] = joinBox( + vec[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ], + True, + False, + ) + + return vec + + +""" job=0 k=np.load('../../data/'+str(job)+'/k.npy') div_veccon(k,100,1,'./') div_veccon(k,100,2,'./') div_veccon(k,100,4,'./') -''' +""" for job in range(6): - k=np.load('../../data/'+str(job)+'/k.npy') - print(job) - res=div_veccon(k,100,4,'./') - np.savetxt('../../data/'+str(job)+'/Cmap_res.txt',res) - res=div_veccon(k,100,1,'./') - #div_veccon(k,100,64,'./') - #div_veccon(k,100,128,'./') + k = np.load("../../data/" + str(job) + "/k.npy") + print(job) + res = div_veccon(k, 100, 4, "./") + np.savetxt("../../data/" + str(job) + "/Cmap_res.txt", res) + res = div_veccon(k, 100, 1, "./") + # div_veccon(k,100,64,'./') + # div_veccon(k,100,128,'./') diff --git a/tools/connec/back/get_cmap.py b/tools/connec/back/get_cmap.py index 76ce203..8a49987 100755 --- a/tools/connec/back/get_cmap.py +++ b/tools/connec/back/get_cmap.py @@ -2,116 +2,136 @@ import numpy as np import os import time -def div_veccon(vec,kh,npartes,condir): - vec=np.where(vec==kh,1,0).astype(int) - Nx, Ny,Nz=k.shape[0],k.shape[1],k.shape[2] - rdir='./' - tt=0 - t1=time.time() - nx=Nx//npartes - params,execCon=ConConfig(nx,Ny,Nz) - - with open(condir+'coninput.txt', 'w') as f: - for item in params: - f.write("%s\n" % item) - - wiam=os.getcwd() - os.chdir(condir) - - - i=0 - np.savetxt(condir+params[2],vec[i*nx:(i+1)*nx,:,:].reshape(-1)) - tcon=time.time() - os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir+ - tt=tt+(time.time()-tcon) - cmap=np.loadtxt(params[-2]).reshape(nx,Ny,Nz).astype(int) - - - - for i in range(1,npartes): - np.savetxt(condir+params[2],vec[i*nx:(i+1)*nx,:,:].reshape(-1)) - tcon=time.time() - os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir++'>/dev/null' - tt=tt+(time.time()-tcon) - cmapb=np.loadtxt(params[-2]).reshape(nx,Ny,Nz).astype(int) - cmap=joinCmap(cmap,cmapb) - - if npartes > 1: - np.savetxt(rdir+'cmap.txt',cmap.reshape(-1)) - - Ttotal, frac_solver = time.time()-t1, tt/(time.time()-t1) - - - y = np.bincount(cmap.reshape(-1)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - if cf[0,0]==0: - cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - nclus=cf.shape[0] #cantidad de clusters - nper=np.sum(cf[:,1]) #num de celdas permeables - print(nclus,float(nper)/(vec.size),Ttotal) - return np.array([npartes,nx*Ny*Nz,Ttotal, frac_solver ,nclus,float(nper)/(Nx*Nx)]) - - -def ConConfig(nx,Ny,Nz): - - params=[] - if Nz==1: - params=['1','4','vecconec.txt',str(nx)+' '+str(Ny),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec2d' - - else: - params=['1','6','vecconec.txt',str(nx)+' '+str(Nz)+' ' +str(Nz),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec3d' - return params, execCon - - -def joinCmap(cmap1,cmap2): - - nclus1 = np.max(cmap1) - cmap2=np.where(cmap2!=0,cmap2+nclus1,0) - - old_nclus=0 - new_nclus=1 - - while new_nclus!= old_nclus: - - old_nclus=new_nclus - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap2=np.where(cmap2==cmap2[0,i,j],cmap1[-1,i,j],cmap2) - - - - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap1=np.where(cmap1==cmap1[-1,i,j],cmap2[0,i,j],cmap1) - - cmap=np.append(cmap1,cmap2,axis=0) - y = np.bincount(cmap.reshape(-1).astype(int)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - new_nclus=cf.shape[0] #cantidad de clusters - #print(new_nclus) - - - - return cmap - - -partes=[1,4] +def div_veccon(vec, kh, npartes, condir): + + vec = np.where(vec == kh, 1, 0).astype(int) + Nx, Ny, Nz = k.shape[0], k.shape[1], k.shape[2] + rdir = "./" + tt = 0 + t1 = time.time() + nx = Nx // npartes + params, execCon = ConConfig(nx, Ny, Nz) + + with open(condir + "coninput.txt", "w") as f: + for item in params: + f.write("%s\n" % item) + + wiam = os.getcwd() + os.chdir(condir) + + i = 0 + np.savetxt(condir + params[2], vec[i * nx : (i + 1) * nx, :, :].reshape(-1)) + tcon = time.time() + os.system(" ./" + execCon + ">/dev/null") #'cd ' +exeDir+ + tt = tt + (time.time() - tcon) + cmap = np.loadtxt(params[-2]).reshape(nx, Ny, Nz).astype(int) + + for i in range(1, npartes): + np.savetxt(condir + params[2], vec[i * nx : (i + 1) * nx, :, :].reshape(-1)) + tcon = time.time() + os.system(" ./" + execCon + ">/dev/null") #'cd ' +exeDir++'>/dev/null' + tt = tt + (time.time() - tcon) + cmapb = np.loadtxt(params[-2]).reshape(nx, Ny, Nz).astype(int) + cmap = joinCmap(cmap, cmapb) + + if npartes > 1: + np.savetxt(rdir + "cmap.txt", cmap.reshape(-1)) + + Ttotal, frac_solver = time.time() - t1, tt / (time.time() - t1) + + y = np.bincount(cmap.reshape(-1)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + if cf[0, 0] == 0: + cf = cf[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + nclus = cf.shape[0] # cantidad de clusters + nper = np.sum(cf[:, 1]) # num de celdas permeables + print(nclus, float(nper) / (vec.size), Ttotal) + return np.array( + [npartes, nx * Ny * Nz, Ttotal, frac_solver, nclus, float(nper) / (Nx * Nx)] + ) + + +def ConConfig(nx, Ny, Nz): + + params = [] + if Nz == 1: + params = [ + "1", + "4", + "vecconec.txt", + str(nx) + " " + str(Ny), + "1.0 1.0", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec2d" + + else: + params = [ + "1", + "6", + "vecconec.txt", + str(nx) + " " + str(Nz) + " " + str(Nz), + "1.0 1.0 1.0", + "30", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec3d" + return params, execCon + + +def joinCmap(cmap1, cmap2): + + nclus1 = np.max(cmap1) + cmap2 = np.where(cmap2 != 0, cmap2 + nclus1, 0) + + old_nclus = 0 + new_nclus = 1 + + while new_nclus != old_nclus: + + old_nclus = new_nclus + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap2 = np.where( + cmap2 == cmap2[0, i, j], cmap1[-1, i, j], cmap2 + ) + + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap1 = np.where( + cmap1 == cmap1[-1, i, j], cmap2[0, i, j], cmap1 + ) + + cmap = np.append(cmap1, cmap2, axis=0) + y = np.bincount(cmap.reshape(-1).astype(int)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + new_nclus = cf.shape[0] # cantidad de clusters + # print(new_nclus) + + return cmap + + +partes = [1, 4] for i in range(1): - t00=time.time() - res=np.array([]) - rdir='../../data/'+str(i)+'/' - k=np.load('k643d.npy') - for npar in partes: - res=np.append(res,div_veccon(k,100,npar,'./')) - np.savetxt(rdir+'resTestCon.txt',res.reshape(len(partes),-1)) - #np.savetxt(rdir+'resTestCon.txt',res.reshape(len(partes),-1)) - print(i,time.time()-t00) + t00 = time.time() + res = np.array([]) + rdir = "../../data/" + str(i) + "/" + k = np.load("k643d.npy") + for npar in partes: + res = np.append(res, div_veccon(k, 100, npar, "./")) + np.savetxt(rdir + "resTestCon.txt", res.reshape(len(partes), -1)) + # np.savetxt(rdir+'resTestCon.txt',res.reshape(len(partes),-1)) + print(i, time.time() - t00) diff --git a/tools/connec/back/get_cmapOpt.py b/tools/connec/back/get_cmapOpt.py index 435898c..a312339 100755 --- a/tools/connec/back/get_cmapOpt.py +++ b/tools/connec/back/get_cmapOpt.py @@ -2,110 +2,127 @@ import numpy as np import os import time -def div_veccon(vec,kh,npartes,condir): - vec=np.where(vec==kh,1,0).astype(int) - Nx, Ny,Nz=k.shape[0],k.shape[1],k.shape[2] - rdir='./' - tt=0 - t1=time.time() - nx=Nx//npartes - params,execCon=ConConfig(nx,Ny,Nz) - - with open(condir+'coninput.txt', 'w') as f: - for item in params: - f.write("%s\n" % item) - - wiam=os.getcwd() - os.chdir(condir) - - - i=0 - np.savetxt(condir+params[2],vec[i*nx:(i+1)*nx,:,:].reshape(-1)) - tcon=time.time() - os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir+ - tt=tt+(time.time()-tcon) - cmap=np.loadtxt(params[-2]).reshape(nx,Ny,Nz) - - - - for i in range(1,npartes): - np.savetxt(condir+params[2],vec[i*nx:(i+1)*nx,:,:].reshape(-1)) - tcon=time.time() - os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir++'>/dev/null' - tt=tt+(time.time()-tcon) - cmapb=np.loadtxt(params[-2]).reshape(nx,Ny,Nz) - cmap=joinCmap(cmap,cmapb) - - if npartes > 1: - np.savetxt(rdir+'cmap.txt',cmap.reshape(-1)) - - Ttotal, frac_solver = time.time()-t1, tt/(time.time()-t1) - - - y = np.bincount(cmap.reshape(-1).astype(int)) - ii = np.nonzero(y)[0] - cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia - if cf[0,0]==0: - cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero - nclus=cf.shape[0] #cantidad de clusters - nper=np.sum(cf[:,1]) #num de celdas permeables - - return np.array([npartes,nx*Ny*Nz,Ttotal, frac_solver ,nclus,float(nper)/(Nx*Nx)]) - - -def ConConfig(nx,Ny,Nz): - - params=[] - if Nz==1: - params=['1','4','vecconec.txt',str(nx)+' '+str(Ny),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec2d' - - else: - params=['1','6','vecconec.txt',str(nx)+' '+str(Nz)+' ' +str(Nz),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF'] - execCon='conec3d' - return params, execCon - - -def joinCmap(cmap1,cmap2): - - nclus1 = np.max(cmap1) - cmap2=np.where(cmap2!=0,cmap2+nclus1,0) - - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap2=np.where(cmap2==cmap2[0,i,j],cmap1[-1,i,j],cmap2) - - - - for i in range(cmap1.shape[1]): - for j in range(cmap1.shape[2]): - if cmap1[-1,i,j] != 0 and cmap2[0,i,j] !=0: - if cmap1[-1,i,j] != cmap2[0,i,j]: - cmap1=np.where(cmap1==cmap1[-1,i,j],cmap2[0,i,j],cmap1) - - cmap=np.append(cmap1,cmap2,axis=0) - - - - return cmap - -njobs=2 -partes=[1,4,8,16] +def div_veccon(vec, kh, npartes, condir): + + vec = np.where(vec == kh, 1, 0).astype(int) + Nx, Ny, Nz = k.shape[0], k.shape[1], k.shape[2] + rdir = "./" + tt = 0 + t1 = time.time() + nx = Nx // npartes + params, execCon = ConConfig(nx, Ny, Nz) + + with open(condir + "coninput.txt", "w") as f: + for item in params: + f.write("%s\n" % item) + + wiam = os.getcwd() + os.chdir(condir) + + i = 0 + np.savetxt(condir + params[2], vec[i * nx : (i + 1) * nx, :, :].reshape(-1)) + tcon = time.time() + os.system(" ./" + execCon + ">/dev/null") #'cd ' +exeDir+ + tt = tt + (time.time() - tcon) + cmap = np.loadtxt(params[-2]).reshape(nx, Ny, Nz) + + for i in range(1, npartes): + np.savetxt(condir + params[2], vec[i * nx : (i + 1) * nx, :, :].reshape(-1)) + tcon = time.time() + os.system(" ./" + execCon + ">/dev/null") #'cd ' +exeDir++'>/dev/null' + tt = tt + (time.time() - tcon) + cmapb = np.loadtxt(params[-2]).reshape(nx, Ny, Nz) + cmap = joinCmap(cmap, cmapb) + + if npartes > 1: + np.savetxt(rdir + "cmap.txt", cmap.reshape(-1)) + + Ttotal, frac_solver = time.time() - t1, tt / (time.time() - t1) + + y = np.bincount(cmap.reshape(-1).astype(int)) + ii = np.nonzero(y)[0] + cf = np.vstack((ii, y[ii])).T # numero de cluster, frecuencia + if cf[0, 0] == 0: + cf = cf[ + 1:, : + ] # me quedo solo con la distr de tamanos, elimino info cluster cero + nclus = cf.shape[0] # cantidad de clusters + nper = np.sum(cf[:, 1]) # num de celdas permeables + + return np.array( + [npartes, nx * Ny * Nz, Ttotal, frac_solver, nclus, float(nper) / (Nx * Nx)] + ) + + +def ConConfig(nx, Ny, Nz): + + params = [] + if Nz == 1: + params = [ + "1", + "4", + "vecconec.txt", + str(nx) + " " + str(Ny), + "1.0 1.0", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec2d" + + else: + params = [ + "1", + "6", + "vecconec.txt", + str(nx) + " " + str(Nz) + " " + str(Nz), + "1.0 1.0 1.0", + "30", + "pardol.STA", + "pardol.CCO", + "pardol.COF", + ] + execCon = "conec3d" + return params, execCon + + +def joinCmap(cmap1, cmap2): + + nclus1 = np.max(cmap1) + cmap2 = np.where(cmap2 != 0, cmap2 + nclus1, 0) + + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap2 = np.where(cmap2 == cmap2[0, i, j], cmap1[-1, i, j], cmap2) + + for i in range(cmap1.shape[1]): + for j in range(cmap1.shape[2]): + if cmap1[-1, i, j] != 0 and cmap2[0, i, j] != 0: + if cmap1[-1, i, j] != cmap2[0, i, j]: + cmap1 = np.where(cmap1 == cmap1[-1, i, j], cmap2[0, i, j], cmap1) + + cmap = np.append(cmap1, cmap2, axis=0) + + return cmap + + +njobs = 2 +partes = [1, 4, 8, 16] for i in range(210): - t00=time.time() - res=np.array([]) - rdir='../../data/'+str(i)+'/' - k=np.load(rdir+'k.npy') - for npar in partes: - res=np.append(res,div_veccon(k,100,npar,'./')) - res=res.reshape(len(partes),-1) - try: - rres=np.loadtxt(rdir+'resTestCon.txt') - res=np.append(rres,res,axis=0) - np.savetxt(rdir+'resTestCon.txt',res) - except: - np.savetxt(rdir+'resTestCon.txt',res) - print(i,time.time()-t00) + t00 = time.time() + res = np.array([]) + rdir = "../../data/" + str(i) + "/" + k = np.load(rdir + "k.npy") + for npar in partes: + res = np.append(res, div_veccon(k, 100, npar, "./")) + res = res.reshape(len(partes), -1) + try: + rres = np.loadtxt(rdir + "resTestCon.txt") + res = np.append(rres, res, axis=0) + np.savetxt(rdir + "resTestCon.txt", res) + except: + np.savetxt(rdir + "resTestCon.txt", res) + print(i, time.time() - t00) diff --git a/tools/connec/comp_cmap.py b/tools/connec/comp_cmap.py index 412584e..1737a1b 100755 --- a/tools/connec/comp_cmap.py +++ b/tools/connec/comp_cmap.py @@ -4,182 +4,233 @@ import time from tools.connec.JoinCmaps import * import subprocess from tools.connec.PostConec import ConnecInd -#k[x,y,z] + +# k[x,y,z] import json -def comp_connec(parser,rundir,nr): - - kc=np.load(rundir+'k.npy') - keep_aspect = parser.get('Connectivity','keep_aspect') - kh,sx = float(parser.get('Generation','kh')),int(parser.get('Connectivity','block_size')) - S_min_post = int(parser.get('Connectivity','indicators_MinBlockSize')) - nimax =2** int(parser.get('Connectivity','Max_sample_size')) - - gcon =bool(parser.get('Connectivity','compGconec')) - - if S_min_post ==-1 or S_min_post > kc.shape[0]: - S_min_post=kc.shape[0] #solo calcula indicadores para mayo escala - if S_min_post ==0: - S_min_post=sx #solo calcula indicadores para escalas a partir del optimo - if sx > S_min_post: - sx = get_min_nbl(kc,nimax,nr,S_min_post) #corta en mas artes para tener mediads de conec - - nbl=kc.shape[0]//sx - - - if keep_aspect=='yes': - keep_aspect=True - else: - keep_aspect=False - - t0=time.time() - kc=np.where(kc==kh,1,0).astype(int) - tcmaps=time.time() - kc=get_smallCmap(kc,nbl,rundir,keep_aspect) - tcmaps=time.time()-tcmaps - - kc,PostConTime=join(kc,nbl,keep_aspect,rundir,S_min_post,gcon) - - ttotal=time.time()-t0 - - summary = np.array([nbl,ttotal,tcmaps/ttotal,PostConTime/ttotal]) - np.savetxt(rundir + 'ConnSummary.txt',summary,header='nbl,ttotal,tcmaps/ttotal,PostConTime/ttotal') - np.save(rundir+'Cmap.npy',kc) - return - - -def get_min_nbl(kc,nimax,nr,smin): - - if kc.shape[2]==1: - dim=2.0 - else: - dim=3.0 - - if nr>0: - y=(1/dim)*np.log2(nr*kc.size/(nimax*(smin**dim))) - else: - y=0 - y=int(y) - s=int((2**y) * smin) - if s connec.out 2>&1') #subprocess.call(['./tools/connec/'+execCon],cwd=rundir) #, '>/dev/null' , cwd=rundir - os.chdir(wd) - vec=np.loadtxt(rundir+params[-1]).reshape(vec.shape[0],vec.shape[1],vec.shape[2]).astype(int) - return vec - -def ConConfig(sx,sy,sz,Nz,rundir): - - params=[] - if Nz==1: - params=['1','4','vecconec.txt',str(sx)+' '+str(sy),'1.0 1.0','pardol.CCO'] - execCon='conec2d' - - else: - params=['1','6','vecconec.txt',str(sx)+' '+str(sy)+' ' +str(sz),'1.0 1.0 1.0','pardol.CCO'] - execCon='conec3d' - - - with open(rundir+'coninput.txt', 'w') as f: - for item in params: - f.write("%s\n" % item) - - return params, execCon - -def join(vec,nbl,keep_aspect,datadir,S_min_post,gcon): - - Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2] - sx = Nx//nbl - - if keep_aspect: - sy,sz = Ny//nbl,Nz//nbl - nblx, nbly,nblz = nbl, nbl, nbl - else: - sy,sz = sx,sx - nblx=nbl - nbly, nblz = Ny//sy, Nz//sz - - ex=np.log2(Nx) - esx=np.log2(sx) - join_z=True - join_y=True - if Nz==1: - sz=1 - nblz=1 - - - post_time=0 - sxL=[sx] - for bs in range(0,int(ex-esx)): - - - if vec.shape[0]==vec.shape[1] and sx>=S_min_post: - t0=time.time() - ConnecInd(vec,[sx],datadir) - post_time=time.time()-t0 - sx,sy,sz = 2*sx,2*sy,2*sz - sxL+=[sx] - - if sz > Nz: - sz=Nz - nblz=1 - join_z=False - if sy > Ny: - sy=Ny - nbly=1 - join_y=False - - nblx,nbly,nblz = Nx//sx, Ny//sy, Nz//sz - - for i in range(nblx): - for j in range(nbly): - for k in range(nblz): - vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=joinBox(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],join_y,join_z) - - if vec.shape[0]==vec.shape[1] and sx>=S_min_post: # - t0=time.time() - ConnecInd(vec,[sx],datadir) - post_time=post_time+(time.time()-t0) - if gcon: - ConnecInd(vec,sxL,datadir+'Global') - return vec, post_time +def comp_connec(parser, rundir, nr): + + kc = np.load(rundir + "k.npy") + keep_aspect = parser.get("Connectivity", "keep_aspect") + kh, sx = float(parser.get("Generation", "kh")), int( + parser.get("Connectivity", "block_size") + ) + S_min_post = int(parser.get("Connectivity", "indicators_MinBlockSize")) + nimax = 2 ** int(parser.get("Connectivity", "Max_sample_size")) + + gcon = bool(parser.get("Connectivity", "compGconec")) + + if S_min_post == -1 or S_min_post > kc.shape[0]: + S_min_post = kc.shape[0] # solo calcula indicadores para mayo escala + if S_min_post == 0: + S_min_post = sx # solo calcula indicadores para escalas a partir del optimo + if sx > S_min_post: + sx = get_min_nbl( + kc, nimax, nr, S_min_post + ) # corta en mas artes para tener mediads de conec + + nbl = kc.shape[0] // sx + + if keep_aspect == "yes": + keep_aspect = True + else: + keep_aspect = False + + t0 = time.time() + kc = np.where(kc == kh, 1, 0).astype(int) + tcmaps = time.time() + kc = get_smallCmap(kc, nbl, rundir, keep_aspect) + tcmaps = time.time() - tcmaps + + kc, PostConTime = join(kc, nbl, keep_aspect, rundir, S_min_post, gcon) + + ttotal = time.time() - t0 + + summary = np.array([nbl, ttotal, tcmaps / ttotal, PostConTime / ttotal]) + np.savetxt( + rundir + "ConnSummary.txt", + summary, + header="nbl,ttotal,tcmaps/ttotal,PostConTime/ttotal", + ) + np.save(rundir + "Cmap.npy", kc) + return + + +def get_min_nbl(kc, nimax, nr, smin): + + if kc.shape[2] == 1: + dim = 2.0 + else: + dim = 3.0 + + if nr > 0: + y = (1 / dim) * np.log2(nr * kc.size / (nimax * (smin ** dim))) + else: + y = 0 + y = int(y) + s = int((2 ** y) * smin) + if s < smin: + s = smin + + return s + + +def get_smallCmap(vec, nbl, rundir, keep_aspect): + + Nx, Ny, Nz = vec.shape[0], vec.shape[1], vec.shape[2] + sx = Nx // nbl + if keep_aspect: + sy, sz = Ny // nbl, Nz // nbl + nblx, nbly, nblz = nbl, nbl, nbl + else: + sy, sz = sx, sx + nblx = nbl + nbly, nblz = Ny // sy, Nz // sz + + params, execCon = ConConfig(sx, sy, sz, Nz, rundir) + if Nz == 1: + nblz = 1 + sz = 1 + os.system("cp ./tools/connec/" + execCon + " " + rundir) + for i in range(nblx): + for j in range(nbly): + for k in range(nblz): + vec[ + i * sx : (i + 1) * sx, j * sy : (j + 1) * sy, k * sz : (k + 1) * sz + ] = connec( + vec[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ], + execCon, + params, + rundir, + ) + + try: + temps = ["pardol*", "conec*d", "coninput.txt", "vecconec.txt"] + for temp in temps: + os.system("rm " + rundir + temp) + except: + print("No connectivity temps to delete") + + return vec + + +def connec(vec, execCon, params, rundir): + np.savetxt(rundir + params[2], vec.reshape(-1), fmt="%i") + wd = os.getcwd() + os.chdir(rundir) + os.system( + "nohup ./" + execCon + " > connec.out 2>&1" + ) # subprocess.call(['./tools/connec/'+execCon],cwd=rundir) #, '>/dev/null' , cwd=rundir + os.chdir(wd) + vec = ( + np.loadtxt(rundir + params[-1]) + .reshape(vec.shape[0], vec.shape[1], vec.shape[2]) + .astype(int) + ) + return vec + + +def ConConfig(sx, sy, sz, Nz, rundir): + + params = [] + if Nz == 1: + params = [ + "1", + "4", + "vecconec.txt", + str(sx) + " " + str(sy), + "1.0 1.0", + "pardol.CCO", + ] + execCon = "conec2d" + + else: + params = [ + "1", + "6", + "vecconec.txt", + str(sx) + " " + str(sy) + " " + str(sz), + "1.0 1.0 1.0", + "pardol.CCO", + ] + execCon = "conec3d" + + with open(rundir + "coninput.txt", "w") as f: + for item in params: + f.write("%s\n" % item) + + return params, execCon + + +def join(vec, nbl, keep_aspect, datadir, S_min_post, gcon): + + Nx, Ny, Nz = vec.shape[0], vec.shape[1], vec.shape[2] + sx = Nx // nbl + + if keep_aspect: + sy, sz = Ny // nbl, Nz // nbl + nblx, nbly, nblz = nbl, nbl, nbl + else: + sy, sz = sx, sx + nblx = nbl + nbly, nblz = Ny // sy, Nz // sz + + ex = np.log2(Nx) + esx = np.log2(sx) + join_z = True + join_y = True + if Nz == 1: + sz = 1 + nblz = 1 + + post_time = 0 + sxL = [sx] + for bs in range(0, int(ex - esx)): + + if vec.shape[0] == vec.shape[1] and sx >= S_min_post: + t0 = time.time() + ConnecInd(vec, [sx], datadir) + post_time = time.time() - t0 + sx, sy, sz = 2 * sx, 2 * sy, 2 * sz + sxL += [sx] + + if sz > Nz: + sz = Nz + nblz = 1 + join_z = False + if sy > Ny: + sy = Ny + nbly = 1 + join_y = False + + nblx, nbly, nblz = Nx // sx, Ny // sy, Nz // sz + + for i in range(nblx): + for j in range(nbly): + for k in range(nblz): + vec[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ] = joinBox( + vec[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ], + join_y, + join_z, + ) + + if vec.shape[0] == vec.shape[1] and sx >= S_min_post: # + t0 = time.time() + ConnecInd(vec, [sx], datadir) + post_time = post_time + (time.time() - t0) + if gcon: + ConnecInd(vec, sxL, datadir + "Global") + return vec, post_time diff --git a/tools/generation/config.py b/tools/generation/config.py index cad562f..90ed0d3 100755 --- a/tools/generation/config.py +++ b/tools/generation/config.py @@ -3,75 +3,74 @@ import configparser import json -def get_config(conffile): - - parser = configparser.ConfigParser() - parser.read(conffile) - - cons=json.loads(parser.get('Iterables',"connectivity")) - ps=json.loads(parser.get('Iterables',"p")) - lcs=json.loads(parser.get('Iterables',"lc")) - variances=json.loads(parser.get('Iterables',"variances")) - seeds=json.loads(parser.get('Iterables',"seeds")) - - seeds=np.arange(seeds[0],seeds[1]+seeds[0]) - ps=np.linspace(ps[0],ps[1],ps[2])/100 - - iterables=dict() - iterables['ps'] = ps - iterables['seeds'] = seeds - iterables['lcs'] = lcs - iterables['variances'] = variances - iterables['cons'] = cons - return parser, iterables - - - -def DotheLoop(job,parser,iterables): - - - - - ps = iterables['ps'] - seeds = iterables['seeds'] - lcs = iterables['lcs'] - variances = iterables['variances'] - cons = iterables['cons'] - - if job==-1: - if parser.get('Generation','binary')=='yes': - if 0 not in cons: - njobs=len(ps)*len(cons)*len(seeds)*len(lcs) - else: - njobs=len(ps)*(len(cons)-1)*len(seeds)*len(lcs)+len(ps)*len(seeds) - else: - if 0 not in cons: - njobs=len(variances)*len(cons)*len(seeds)*len(lcs) - else: - njobs=len(variances)*(len(cons)-1)*len(seeds)*len(lcs)+len(variances)*len(seeds) - return njobs - - i=0 - - for con in cons: - if con == 0: - llcs=[0.000001] - else: - llcs=lcs - for lc in llcs: - if parser.get('Generation','binary')=='yes': - for p in ps: - for seed in seeds: - if i==job: - return [con,lc,p,seed] - i+=1 - - else: - for v in variances: - for seed in seeds: - if i==job: - return [con,lc,v,seed] - i+=1 - return [] +def get_config(conffile): + parser = configparser.ConfigParser() + parser.read(conffile) + + cons = json.loads(parser.get("Iterables", "connectivity")) + ps = json.loads(parser.get("Iterables", "p")) + lcs = json.loads(parser.get("Iterables", "lc")) + variances = json.loads(parser.get("Iterables", "variances")) + seeds = json.loads(parser.get("Iterables", "seeds")) + + seeds = np.arange(seeds[0], seeds[1] + seeds[0]) + ps = np.linspace(ps[0], ps[1], ps[2]) / 100 + + iterables = dict() + iterables["ps"] = ps + iterables["seeds"] = seeds + iterables["lcs"] = lcs + iterables["variances"] = variances + iterables["cons"] = cons + return parser, iterables + + +def DotheLoop(job, parser, iterables): + + ps = iterables["ps"] + seeds = iterables["seeds"] + lcs = iterables["lcs"] + variances = iterables["variances"] + cons = iterables["cons"] + + if job == -1: + if parser.get("Generation", "binary") == "yes": + if 0 not in cons: + njobs = len(ps) * len(cons) * len(seeds) * len(lcs) + else: + njobs = len(ps) * (len(cons) - 1) * len(seeds) * len(lcs) + len( + ps + ) * len(seeds) + else: + if 0 not in cons: + njobs = len(variances) * len(cons) * len(seeds) * len(lcs) + else: + njobs = len(variances) * (len(cons) - 1) * len(seeds) * len(lcs) + len( + variances + ) * len(seeds) + return njobs + + i = 0 + + for con in cons: + if con == 0: + llcs = [0.000001] + else: + llcs = lcs + for lc in llcs: + if parser.get("Generation", "binary") == "yes": + for p in ps: + for seed in seeds: + if i == job: + return [con, lc, p, seed] + i += 1 + + else: + for v in variances: + for seed in seeds: + if i == job: + return [con, lc, v, seed] + i += 1 + return [] diff --git a/tools/generation/fftma_gen.py b/tools/generation/fftma_gen.py index 483ebde..ea09dad 100755 --- a/tools/generation/fftma_gen.py +++ b/tools/generation/fftma_gen.py @@ -3,156 +3,238 @@ from scipy.special import erf, erfinv from scipy.optimize import curve_fit from scipy.stats import norm from FFTMA import gen -from tools.generation.config import DotheLoop, get_config +from tools.generation.config import DotheLoop, get_config from tools.generation.variograma import get_lc from scipy.interpolate import interp1d import sys import time import os -#from memory_profiler import profile -def fftmaGenerator(datadir,job,conffile): +# from memory_profiler import profile + + +def fftmaGenerator(datadir, job, conffile): + + t0 = time.time() + parser, iterables = get_config(conffile) + params = DotheLoop(job, parser, iterables) + + binary = parser.get("Generation", "binary") + uselc_bin = parser.get("Generation", "lcBin") + + if binary == "yes": + logn = "no" + con, lc, p, seed = params[0], params[1], params[2], params[3] + variance = 0 + else: + logn = "yes" + con, lc, variance, seed = params[0], params[1], params[2], params[3] + p = 0 + + Nx, Ny, Nz = ( + int(parser.get("Generation", "Nx")), + int(parser.get("Generation", "Ny")), + int(parser.get("Generation", "Nz")), + ) + # N=int(42.666666667*lc) + # Nx,Ny,Nz = N,N,N + # print(N) + + kh, kl, vario = ( + float(parser.get("Generation", "kh")), + float(parser.get("Generation", "kl")), + int(parser.get("Generation", "variogram_type")), + ) + compute_lc = parser.get("Generation", "compute_lc") + generate_K( + Nx, + Ny, + Nz, + con, + lc, + p, + kh, + kl, + seed, + logn, + variance, + vario, + datadir, + compute_lc, + uselc_bin, + ) + + np.savetxt( + datadir + "GenParams.txt", + np.array( + [time.time() - t0, Nx, Ny, Nz, con, lc, p, kh, kl, seed, variance, vario] + ), + header="Runtime,Nx,Ny,Nz,con,lc,p,kh,kl,seed,variance,vario", + ) + + return + + +def obtainLctobin(p, con, vario): + + lc = np.load( + "./tools/generation/lc.npy", allow_pickle=True, encoding="latin1" + ).item() + + f = interp1d(lc["p"], lc[vario, con]) + if p == 0 or p == 1: + return 1.0 + return 1.0 / f(p) + + +def obtainLctobinBack(p, con, vario): + + pb = np.linspace(0.0, 1.0, 11) + + if vario == 2: + i = [0.0, 1.951, 2.142, 2.247, 2.301, 2.317, 2.301, 2.246, 2.142, 1.952, 0.0] + c = [0.0, 1.188, 1.460, 1.730, 2.017, 2.284, 2.497, 2.652, 2.736, 2.689, 0.0] + d = [0.0, 2.689, 2.736, 2.652, 2.497, 2.284, 2.017, 1.730, 1.460, 1.188, 0.0] + lcBin = np.array([i, c, d]) + lcBin = lcBin / 3.0 + + if vario == 1: + i = [0.0, 3.13, 3.66, 3.94, 4.08, 4.10, 4.01, 3.84, 3.55, 3.00, 0.0] + c = [0.0, 0.85, 1.095, 1.312, 1.547, 1.762, 1.966, 2.149, 2.257, 2.186, 0.0] + d = [ + 0.0, + 2.186, + 2.2575, + 2.1495, + 1.9660, + 1.7625, + 1.5476, + 1.3128, + 1.0950, + 0.8510, + 0.0, + ] + lcBin = np.array([i, c, d]) + lcBin = lcBin / 6.0 + + f = interp1d(pb, lcBin[con - 1]) + return 1.0 / f(p) + + +# @profile +def generate_K( + Nx, + Ny, + Nz, + con, + lc, + p, + kh, + kl, + seed, + logn, + LogVariance, + vario, + datadir, + compute_lc, + uselc_bin, +): + + k = genGaussK( + Nx, Ny, Nz, con, lc, p, kh, kl, seed, logn, LogVariance, vario, uselc_bin + ) + if compute_lc == "yes": + lcG = get_lc(k, vario) + lcNst = lcG + lcBin = np.nan + if con == 2: + k = -nst(k) # normal score transform + if compute_lc == "yes": + lcNst = get_lc(k, vario) + if con == 3: + k = nst(k) + if compute_lc == "yes": + lcNst = get_lc(k, vario) + + if logn == "yes": + k = k * (LogVariance ** 0.5) + k = np.exp(k) + + else: + k = binarize(k, kh, kl, p) + if compute_lc == "yes": + lcBin = get_lc(np.where(k > kl, 1, 0), vario) + np.save(datadir + "k.npy", k) + if compute_lc == "yes": + np.savetxt( + datadir + "lc.txt", + np.array([lcG, lcNst, lcBin]), + header="lcG, lcNst, lcBin", + ) + return + + +def genGaussK( + Nx, Ny, Nz, con, lc, p, kh, kl, seed, logn, LogVariance, vario, uselc_bin +): + + typ = 0 # structure du champ: 0=normal; 1=lognormal; 2=log-10 + dx, dy, dz = 1.0, 1.0, 1.0 + var = 1 # Nbr de structure du variogramme + alpha = 1 # valeur exposant + if con == 0: + lc = 0.000001 + + if (con == 2 or con == 3) and vario == 2: + lc = lc / 0.60019978939 + if (con == 2 or con == 3) and vario == 1: + lc = lc / 0.38165155120015 + + if uselc_bin == "yes" and con != 0: + lc = lc * obtainLctobin(p, con, vario) + v1 = ( + var, + vario, + alpha, + lc, + lc, + lc, + 1, + 0, + 0, + 0, + 1, + 0, + ) # coord des vecteurs de base (1 0 0) y (0 1 0) + k = gen( + Nz, Ny, Nx, dx, dy, dz, seed, [v1], 0, 1, 0 + ) # 0, 1, 0 = mean, variance, typ #Generation of a correlated standard dsitribution N(0,1) + + return k - t0=time.time() - parser, iterables = get_config(conffile) - params = DotheLoop(job,parser, iterables ) - - binary=parser.get('Generation','binary') - uselc_bin=parser.get('Generation','lcBin') - - if binary=='yes': - logn='no' - con,lc,p,seed = params[0],params[1],params[2],params[3] - variance=0 - else: - logn='yes' - con,lc,variance,seed = params[0],params[1],params[2],params[3] - p=0 - - - Nx,Ny,Nz = int(parser.get('Generation','Nx')), int(parser.get('Generation','Ny')), int(parser.get('Generation','Nz')) - #N=int(42.666666667*lc) - #Nx,Ny,Nz = N,N,N - #print(N) - - kh,kl,vario = float(parser.get('Generation','kh')),float(parser.get('Generation','kl')), int(parser.get('Generation','variogram_type')) - compute_lc=parser.get('Generation','compute_lc') - generate_K(Nx,Ny,Nz,con,lc,p,kh,kl,seed,logn,variance,vario,datadir,compute_lc,uselc_bin) - - np.savetxt(datadir+'GenParams.txt',np.array([time.time()-t0,Nx,Ny,Nz,con,lc,p,kh,kl,seed,variance,vario]),header='Runtime,Nx,Ny,Nz,con,lc,p,kh,kl,seed,variance,vario') - - - return - - -def obtainLctobin(p,con,vario): - - - lc=np.load('./tools/generation/lc.npy',allow_pickle=True, encoding = 'latin1').item() - - f=interp1d(lc['p'],lc[vario,con]) - if p==0 or p==1: - return 1.0 - return 1.0/f(p) - -def obtainLctobinBack(p,con,vario): - - pb=np.linspace(0.0,1.0,11) - - - if vario==2: - i=[0.0, 1.951, 2.142, 2.247, 2.301, 2.317, 2.301, 2.246, 2.142, 1.952, 0.0] - c=[0.0, 1.188, 1.460, 1.730, 2.017, 2.284, 2.497, 2.652, 2.736, 2.689, 0.0] - d=[0.0,2.689, 2.736,2.652, 2.497, 2.284, 2.017, 1.730, 1.460, 1.188, 0.0] - lcBin=np.array([i,c,d]) - lcBin=lcBin/3.0 - - if vario==1: - i=[0.0,3.13, 3.66, 3.94, 4.08, 4.10, 4.01, 3.84, 3.55, 3.00,0.0] - c=[0.0,0.85, 1.095, 1.312, 1.547, 1.762, 1.966, 2.149, 2.257, 2.186,0.0] - d=[0.0,2.186, 2.2575,2.1495,1.9660,1.7625,1.5476,1.3128,1.0950,0.8510,0.0] - lcBin=np.array([i,c,d]) - lcBin=lcBin/6.0 - - f=interp1d(pb,lcBin[con-1]) - return 1.0/f(p) - - -#@profile -def generate_K(Nx,Ny,Nz,con,lc,p,kh,kl,seed,logn,LogVariance,vario,datadir,compute_lc,uselc_bin): - - - k=genGaussK(Nx,Ny,Nz,con,lc,p,kh,kl,seed,logn,LogVariance,vario,uselc_bin) - if compute_lc =='yes': - lcG=get_lc(k,vario) - lcNst=lcG - lcBin=np.nan - if con==2: - k = -nst(k) #normal score transform - if compute_lc =='yes': - lcNst=get_lc(k,vario) - if con==3: - k = nst(k) - if compute_lc =='yes': - lcNst=get_lc(k,vario) - - if logn == 'yes': - k=k*(LogVariance**0.5) - k = np.exp(k) - - else: - k = binarize(k,kh,kl,p) - if compute_lc =='yes': - lcBin=get_lc(np.where(k>kl,1,0),vario) - np.save(datadir+'k.npy',k) - if compute_lc =='yes': - np.savetxt(datadir+'lc.txt',np.array([lcG,lcNst,lcBin]),header='lcG, lcNst, lcBin') - return - -def genGaussK(Nx,Ny,Nz,con,lc,p,kh,kl,seed,logn,LogVariance,vario,uselc_bin): +def nst(kc): + kc = np.abs(kc) + kc = np.sqrt(2) * erfinv(2 * erf(kc / np.sqrt(2)) - 1) + return kc - typ=0 #structure du champ: 0=normal; 1=lognormal; 2=log-10 - dx, dy, dz = 1.0, 1.0, 1.0 - var=1 #Nbr de structure du variogramme - alpha=1 #valeur exposant - if con==0: - lc=0.000001 - if (con==2 or con==3) and vario==2: - lc=lc/0.60019978939 - if (con==2 or con==3) and vario==1: - lc=lc/0.38165155120015 +def binarize(kc, kh, kl, p): - if uselc_bin=='yes' and con!=0: - lc=lc*obtainLctobin(p,con,vario) - v1 = (var, vario, alpha, lc, lc, lc, 1, 0, 0, 0, 1, 0) # coord des vecteurs de base (1 0 0) y (0 1 0) - k=gen(Nz, Ny, Nx, dx, dy, dz, seed, [v1], 0, 1, 0) # 0, 1, 0 = mean, variance, typ #Generation of a correlated standard dsitribution N(0,1) + if kc.size < 100 ** 3: + if p > 0: + at = int((1 - p) * kc.size) # + else: + at = kc.size - 1 + t1 = np.sort(kc.reshape(-1))[at] # get permeability treshold + kc = np.where(kc < t1, kl, kh) # Binarization + t1 = 0 + else: + t1 = norm.ppf(1 - p) + kc = np.where(kc < t1, kl, kh) + return kc - return k +# CONFIG_FILE_PATH = 'config.ini' if 'CONFIG_FILE_PATH' not in os.environ else os.environ['CONFIG_FILE_PATH'] -def nst(kc): - kc=np.abs(kc) - kc=np.sqrt(2)*erfinv(2*erf(kc/np.sqrt(2))-1) - return kc - - -def binarize(kc,kh,kl,p): - - if kc.size < 100**3: - if p>0: - at=int((1-p)*kc.size) # - else: - at=kc.size-1 - t1=np.sort(kc.reshape(-1))[at] #get permeability treshold - kc=np.where(kc0: - y=(1/dim)*np.log2(nr*kc.size/(nimax*(smin**dim))) - else: - y=0 - y=int(y) - s=int((2**y) * smin) - if s 0: + y = (1 / dim) * np.log2(nr * kc.size / (nimax * (smin ** dim))) + else: + y = 0 + y = int(y) + s = int((2 ** y) * smin) + if s < smin: + s = smin + + return s diff --git a/tools/postprocessK/comp_PostKeff.py b/tools/postprocessK/comp_PostKeff.py index 7c808a8..9e2acc8 100755 --- a/tools/postprocessK/comp_PostKeff.py +++ b/tools/postprocessK/comp_PostKeff.py @@ -2,120 +2,145 @@ import numpy as np import os import time from tools.postprocessK.flow import ComputeVol, comp_Kdiss_Kaverage -#import subprocess -from tools.postprocessK.kperm.Ndar1P import PetscP -#k[x,y,z] -import json - -def comp_postKeff(parser,rundir,nr): - - - - k=np.load(rundir+'k.npy') - try: - P=np.load(rundir+'P.npy') - except: - print('no pressure file '+rundir) - return - ref=P.shape[0]//k.shape[0] - - SaveV = parser.get('K-Postprocess','SaveVfield') - if SaveV=='yes': - SaveV=True - else: - SaveV=False - - t0=time.time() - k, diss, vx,vy,vz, Px, Py, Pz = ComputeVol(k,P,SaveV) #refina k - tDissVel=time.time()-t0 - - P=0 - - - - S_min_post = int(parser.get('K-Postprocess','MinBlockSize')) - nimax =2** int(parser.get('K-Postprocess','Max_sample_size')) - compKperm =parser.get('K-Postprocess','kperm') - if compKperm=='yes': - compKperm=True - - S_min_post=S_min_post*ref - - if S_min_post==0: - sx=1 #k.shape[0] - else: - sx = get_min_nbl(k,nimax,nr,S_min_post) - - - kdiss,kave=getKpost(k, diss, vx, Px, Py, Pz,sx,rundir,ref,compKperm) - - - ttotal=time.time()-t0 - - summary = np.array([kdiss,kave,ttotal,tDissVel/ttotal]).T - np.savetxt(rundir + 'PosKeffSummary.txt',summary,header='K_diss, K_average,ttotal,tDiss/ttotal') - if SaveV: - np.save(rundir+'V.npy',np.array([vx,vy,vz])) - np.save(rundir+'D.npy',diss) - return +# import subprocess +from tools.postprocessK.kperm.Ndar1P import PetscP +# k[x,y,z] +import json -def getKpost(kf, diss, vx, Px, Py, Pz,sx,rundir,ref,compkperm): - - - - ex=int(np.log2(kf.shape[0])) - esx=int(np.log2(sx)) - - scales=2**np.arange(esx,ex) - datadir=rundir+'KpostProcess/' - try: - os.makedirs(datadir) - except: - nada=0 - - for l in scales: - nblx, nbly, nblz = kf.shape[0]//l, kf.shape[1]//l, kf.shape[2]//l - sx,sy,sz=l,l,l - if kf.shape[2]==1: - nblz=1 - sz=1 - - Kdiss,Kave=np.zeros((nblx,nbly,nblz)),np.zeros((nblx,nbly,nblz)) - - for i in range(nblx): - for j in range(nbly): - for k in range(nblz): - Kdiss[i,j,k],Kave[i,j,k]=comp_Kdiss_Kaverage(kf[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz], diss[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz], vx[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz], Px[i*sx:(i+1)*sx+1,j*sy:(j+1)*sy+1,k*sz:(k+1)*sz+1], Py[i*sx:(i+1)*sx+1,j*sy:(j+1)*sy+1,k*sz:(k+1)*sz+1], Pz[i*sx:(i+1)*sx+1,j*sy:(j+1)*sy+1,k*sz:(k+1)*sz+1]) - - - np.save(datadir+'Kd'+str(l//ref)+'.npy',Kdiss) - np.save(datadir+'Kv'+str(l//ref)+'.npy',Kave) - - Kdiss,Kave = comp_Kdiss_Kaverage(kf, diss, vx, Px, Py, Pz) - np.save(datadir+'Kd'+str(kf.shape[0]//ref)+'.npy',np.array([Kdiss])) - np.save(datadir+'Kv'+str(kf.shape[0]//ref)+'.npy',np.array([Kave])) - - - return Kdiss, Kave - - -def get_min_nbl(kc,nimax,nr,smin): - - if kc.shape[2]==1: - dim=2.0 - else: - dim=3.0 - if nr>0: - y=(1/dim)*np.log2(nr*kc.size/(nimax*(smin**dim))) - else: - y=0 - y=int(y) - s=int((2**y) * smin) - if s 0: + y = (1 / dim) * np.log2(nr * kc.size / (nimax * (smin ** dim))) + else: + y = 0 + y = int(y) + s = int((2 ** y) * smin) + if s < smin: + s = smin + + return s diff --git a/tools/postprocessK/flow.py b/tools/postprocessK/flow.py index 5c6fb71..b421f41 100755 --- a/tools/postprocessK/flow.py +++ b/tools/postprocessK/flow.py @@ -2,158 +2,199 @@ import numpy as np from scipy.sparse import diags from scipy.stats import mstats -from scipy.sparse.linalg import bicg, bicgstab, cg, dsolve #,LinearOperator, spilu, bicgstab -#from scikits.umfpack import spsolve, splu +from scipy.sparse.linalg import ( + bicg, + bicgstab, + cg, + dsolve, +) # ,LinearOperator, spilu, bicgstab + +# from scikits.umfpack import spsolve, splu import time -def getDiss(k,vx,vy,vz): - diss = (vx[1:,:,:]**2+vx[:-1,:,:]**2+vy[:,1:,:]**2+vy[:,:-1,:]**2+vz[:,:,1:]**2+vz[:,:,:-1]**2)/(2*k) - return diss + +def getDiss(k, vx, vy, vz): + diss = ( + vx[1:, :, :] ** 2 + + vx[:-1, :, :] ** 2 + + vy[:, 1:, :] ** 2 + + vy[:, :-1, :] ** 2 + + vz[:, :, 1:] ** 2 + + vz[:, :, :-1] ** 2 + ) / (2 * k) + return diss -def ComputeVol(k,P,saveV): +def ComputeVol(k, P, saveV): - k=refina(k, P.shape[0]//k.shape[0]) - Px,Py,Pz = getPfaces(k,P) - vx,vy,vz = getVfaces(k,P, Px,Py, Pz) - diss = getDiss(k,vx,vy,vz) - if saveV==False: - vy, vz= 0, 0 - else: - vy, vz= 0.5*(vy[:,1:,:]+vy[:,:-1,:]), 0.5*(vz[:,:,1:]+vz[:,:,:-1]) - vx= 0.5*(vx[1:,:,:]+vx[:-1,:,:]) + k = refina(k, P.shape[0] // k.shape[0]) + Px, Py, Pz = getPfaces(k, P) + vx, vy, vz = getVfaces(k, P, Px, Py, Pz) + diss = getDiss(k, vx, vy, vz) + if saveV == False: + vy, vz = 0, 0 + else: + vy, vz = 0.5 * (vy[:, 1:, :] + vy[:, :-1, :]), 0.5 * ( + vz[:, :, 1:] + vz[:, :, :-1] + ) + vx = 0.5 * (vx[1:, :, :] + vx[:-1, :, :]) + + return k, diss, vx, vy, vz, Px, Py, Pz - return k, diss, vx,vy,vz, Px, Py, Pz def comp_Kdiss_Kaverage(k, diss, vx, Px, Py, Pz): - mgx, mgy, mgz = np.mean(Px[-1,:,:]-Px[0,:,:])/k.shape[0],np.mean(Py[:,-1,:]-Py[:,0,:])/k.shape[1],np.mean(Pz[:,:,-1]-Pz[:,:,0])/k.shape[2] - kave=np.mean(vx)/mgx - kdiss=np.mean(diss)/(mgx**2+mgy**2+mgz**2) - return kdiss, kave + mgx, mgy, mgz = ( + np.mean(Px[-1, :, :] - Px[0, :, :]) / k.shape[0], + np.mean(Py[:, -1, :] - Py[:, 0, :]) / k.shape[1], + np.mean(Pz[:, :, -1] - Pz[:, :, 0]) / k.shape[2], + ) + kave = np.mean(vx) / mgx + kdiss = np.mean(diss) / (mgx ** 2 + mgy ** 2 + mgz ** 2) + return kdiss, kave +def getKeff(pm, k, pbc, Nz): -def getKeff(pm,k,pbc,Nz): + nx = k.shape[2] # Pasar k sin bordes de k=0 + ny = k.shape[1] - nx = k.shape[2] #Pasar k sin bordes de k=0 - ny = k.shape[1] + tz = 2 * k[1, :, :] * k[0, :, :] / (k[0, :, :] + k[1, :, :]) + q = ((pm[0, :, :] - pm[1, :, :]) * tz).sum() + area = ny * nx + l = Nz + keff = q * l / (pbc * area) + return keff, q - tz = 2*k[1,:,:]*k[0, :,:]/(k[0, :,:]+k[1,:,:]) - q=((pm[0,:,:]-pm[1,:,:])*tz).sum() - area=ny*nx - l=Nz - keff=q*l/(pbc*area) - return keff,q -def getPfaces(k,P): - nx,ny,nz=k.shape[0],k.shape[1],k.shape[2] - Px,Py,Pz= np.zeros((nx+1,ny,nz)),np.zeros((nx,ny+1,nz)),np.zeros((nx,ny,nz+1)) +def getPfaces(k, P): + nx, ny, nz = k.shape[0], k.shape[1], k.shape[2] + Px, Py, Pz = ( + np.zeros((nx + 1, ny, nz)), + np.zeros((nx, ny + 1, nz)), + np.zeros((nx, ny, nz + 1)), + ) - Px[1:-1,:,:] = (k[:-1,:,:]*P[:-1,:,:]+k[1:,:,:]*P[1:,:,:])/(k[:-1,:,:]+k[1:,:,:]) - Px[0,:,:]=nx + Px[1:-1, :, :] = (k[:-1, :, :] * P[:-1, :, :] + k[1:, :, :] * P[1:, :, :]) / ( + k[:-1, :, :] + k[1:, :, :] + ) + Px[0, :, :] = nx - Py[:,1:-1,:] = (k[:,:-1,:]*P[:,:-1,:]+k[:,1:,:]*P[:,1:,:])/(k[:,:-1,:]+k[:,1:,:]) - Py[:,0,:],Py[:,-1,:] =P[:,0,:], P[:,-1,:] + Py[:, 1:-1, :] = (k[:, :-1, :] * P[:, :-1, :] + k[:, 1:, :] * P[:, 1:, :]) / ( + k[:, :-1, :] + k[:, 1:, :] + ) + Py[:, 0, :], Py[:, -1, :] = P[:, 0, :], P[:, -1, :] - Pz[:,:,1:-1] = (k[:,:,:-1]*P[:,:,:-1]+k[:,:,1:]*P[:,:,1:])/(k[:,:,:-1]+k[:,:,1:]) - Pz[:,:,0],Pz[:,:,-1] =P[:,:,0], P[:,:,-1] + Pz[:, :, 1:-1] = (k[:, :, :-1] * P[:, :, :-1] + k[:, :, 1:] * P[:, :, 1:]) / ( + k[:, :, :-1] + k[:, :, 1:] + ) + Pz[:, :, 0], Pz[:, :, -1] = P[:, :, 0], P[:, :, -1] - return Px, Py, Pz + return Px, Py, Pz -def getVfaces(k,P, Px,Py, Pz): - nx,ny,nz=k.shape[0],k.shape[1],k.shape[2] - vx,vy,vz= np.zeros((nx+1,ny,nz)),np.zeros((nx,ny+1,nz)),np.zeros((nx,ny,nz+1)) - vx[1:,:,:] = 2*k*(Px[1:,:,:]-P) #v= k*(deltaP)/(deltaX/2) - vx[0,:,:] = 2*k[0,:,:]*(P[0,:,:]-Px[0,:,:]) +def getVfaces(k, P, Px, Py, Pz): + nx, ny, nz = k.shape[0], k.shape[1], k.shape[2] + vx, vy, vz = ( + np.zeros((nx + 1, ny, nz)), + np.zeros((nx, ny + 1, nz)), + np.zeros((nx, ny, nz + 1)), + ) + vx[1:, :, :] = 2 * k * (Px[1:, :, :] - P) # v= k*(deltaP)/(deltaX/2) + vx[0, :, :] = 2 * k[0, :, :] * (P[0, :, :] - Px[0, :, :]) - vy[:,1:,:] = 2*k*(Py[:,1:,:]-P) - vy[:,0,:] = 2*k[:,0,:]*(P[:,0,:]-Py[:,0,:]) - vz[:,:,1:] = 2*k*(Pz[:,:,1:]-P) - vz[:,:,0] = 2*k[:,:,0]*(P[:,:,0]-Pz[:,:,0]) - return vx,vy,vz + vy[:, 1:, :] = 2 * k * (Py[:, 1:, :] - P) + vy[:, 0, :] = 2 * k[:, 0, :] * (P[:, 0, :] - Py[:, 0, :]) + vz[:, :, 1:] = 2 * k * (Pz[:, :, 1:] - P) + vz[:, :, 0] = 2 * k[:, :, 0] * (P[:, :, 0] - Pz[:, :, 0]) + return vx, vy, vz def refina(k, ref): - if ref==1: - return k - - nx,ny,nz=k.shape[0],k.shape[1],k.shape[2] + if ref == 1: + return k + + nx, ny, nz = k.shape[0], k.shape[1], k.shape[2] - krx=np.zeros((ref*nx,ny,nz)) - for i in range(ref): - krx[i::ref,:,:]=k - k=0 - krxy=np.zeros((ref*nx,ny*ref,nz)) - - for i in range(ref): - krxy[:,i::ref,:]=krx - krx=0 - if nz==1: - return krxy + krx = np.zeros((ref * nx, ny, nz)) + for i in range(ref): + krx[i::ref, :, :] = k + k = 0 + krxy = np.zeros((ref * nx, ny * ref, nz)) + for i in range(ref): + krxy[:, i::ref, :] = krx + krx = 0 + if nz == 1: + return krxy - krxyz=np.zeros((ref*nx,ny*ref,nz*ref)) - for i in range(ref): - krxyz[:,:,i::ref]=krxy - krxy=0 + krxyz = np.zeros((ref * nx, ny * ref, nz * ref)) + for i in range(ref): + krxyz[:, :, i::ref] = krxy + krxy = 0 - return krxyz + return krxyz def computeT(k): - nx = k.shape[0] - ny = k.shape[1] - nz = k.shape[2] - tx = np.zeros((nx+1,ny, nz)) - ty = np.zeros((nx,ny+1, nz)) - tz = np.zeros((nx,ny, nz+1)) - tx[1:-1,:,:] = 2*k[:-1,:,:]*k[1:,:,:]/(k[:-1,:,:]+k[1:,:,:]) - ty[:,1:-1,:] = 2*k[:,:-1,:]*k[:,1:,:]/(k[:,:-1,:]+k[:,1:,:]) - tz[:,:,1:-1] = 2*k[:,:,:-1]*k[:,:,1:]/(k[:,:,:-1]+k[:,:,1:]) + nx = k.shape[0] + ny = k.shape[1] + nz = k.shape[2] + tx = np.zeros((nx + 1, ny, nz)) + ty = np.zeros((nx, ny + 1, nz)) + tz = np.zeros((nx, ny, nz + 1)) + tx[1:-1, :, :] = 2 * k[:-1, :, :] * k[1:, :, :] / (k[:-1, :, :] + k[1:, :, :]) + ty[:, 1:-1, :] = 2 * k[:, :-1, :] * k[:, 1:, :] / (k[:, :-1, :] + k[:, 1:, :]) + tz[:, :, 1:-1] = 2 * k[:, :, :-1] * k[:, :, 1:] / (k[:, :, :-1] + k[:, :, 1:]) - return tx, ty, tz + return tx, ty, tz def Rmat(k): + pbc = k.shape[0] + tx, ty, tz = computeT(k) - pbc=k.shape[0] - tx, ty , tz = computeT(k) + tx[0, :, :], tx[-1, :, :] = 2 * tx[1, :, :], 2 * tx[-2, :, :] - tx[0,:,:],tx[-1,:,:] = 2*tx[1,:,:],2*tx[-2,:,:] + rh = np.zeros((k.shape[0], k.shape[1], k.shape[2])) - rh=np.zeros((k.shape[0],k.shape[1],k.shape[2])) + rh[0, :, :] = pbc * tx[0, :, :] + rh = rh.reshape(-1) + d = ( + tz[:, :, :-1] + + tz[:, :, 1:] + + ty[:, :-1, :] + + ty[:, 1:, :] + + tx[:-1, :, :] + + tx[1:, :, :] + ).reshape(-1) + a = (-tz[:, :, :-1].reshape(-1))[1:] + # a=(tx.reshape(-1))[:-1] + b = (-ty[:, 1:, :].reshape(-1))[: -k.shape[2]] + c = -tx[1:-1, :, :].reshape(-1) - rh[0,:,:]=pbc*tx[0,:,:] - rh=rh.reshape(-1) - d=(tz[:,:,:-1]+tz[:,:,1:]+ty[:,:-1,:]+ty[:,1:,:]+tx[:-1,:,:]+tx[1:,:,:]).reshape(-1) - a=(-tz[:,:,:-1].reshape(-1))[1:] - #a=(tx.reshape(-1))[:-1] - b=(-ty[:,1:,:].reshape(-1))[:-k.shape[2]] - c=-tx[1:-1,:,:].reshape(-1) - - - return a, b, c, d, rh + return a, b, c, d, rh def PysolveP(k, solver): - a, b, c, d, rh = Rmat(k) - nx, ny, nz = k.shape[0], k.shape[1],k.shape[2] - offset = [-nz*ny,-nz, -1, 0, 1, nz, nz*ny] - km=diags(np.array([c, b, a, d, a, b, c]), offset, format='csc') - a, b, c, d = 0, 0 ,0 , 0 - lu = splu(km) - print(lu) - p = solver(km, rh) - p=p.reshape(nx, ny, nz) - keff,q = getKeff(p,k,nz,nz) - return keff -''' + a, b, c, d, rh = Rmat(k) + nx, ny, nz = k.shape[0], k.shape[1], k.shape[2] + offset = [-nz * ny, -nz, -1, 0, 1, nz, nz * ny] + km = diags(np.array([c, b, a, d, a, b, c]), offset, format="csc") + a, b, c, d = 0, 0, 0, 0 + lu = splu(km) + print(lu) + p = solver(km, rh) + p = p.reshape(nx, ny, nz) + keff, q = getKeff(p, k, nz, nz) + return keff + + +""" solvers=[bicg, bicgstab, cg, dsolve, spsolve] snames=['bicg', 'bicgstab',' cg',' dsolve',' spsolve'] @@ -168,5 +209,4 @@ for job in range(jobs): keff=PysolveP(kff, solvers[i]) print('Solver: '+snames[i]+' time: '+str(time.time()-t0)) -''' - +""" diff --git a/tools/postprocessK/kperm/Ndar1P.py b/tools/postprocessK/kperm/Ndar1P.py index 4b39d92..11dc38b 100755 --- a/tools/postprocessK/kperm/Ndar1P.py +++ b/tools/postprocessK/kperm/Ndar1P.py @@ -2,88 +2,77 @@ import numpy as np import petsc4py import math import time -#from mpi4py import MPI + +# from mpi4py import MPI from tools.postprocessK.kperm.computeFlows import * from petsc4py import PETSc -petsc4py.init('-ksp_max_it 9999999999') -from tools.postprocessK.kperm.flow import getKeff - - - - -def PetscP(datadir,ref,k,saveres): - - ref=1 - rank=0 - pn=1 - t0=time.time() - - pcomm=PETSc.COMM_SELF - if k.shape[2]==1: - refz=1 - else: - refz=ref - - nz, ny, nx=k.shape[0]*ref,k.shape[1]*ref,k.shape[2]*refz - n=nx*ny*nz - - - - K = PETSc.Mat().create(comm=pcomm) - K.setType('seqaij') - K.setSizes(((n,None),(n,None))) # Aca igual que lo que usas arriba - K.setPreallocationNNZ(nnz=(7,4)) # Idem anterior - K.setUp() - - R = PETSc.Vec().createSeq((n,None),comm=pcomm) #PETSc.COMM_WORLD - R.setUp() - - k2, Nz, nnz2=getKref(k,1,2,ref) - k, Nz, nnz=getKref(k,0,2,ref) +petsc4py.init("-ksp_max_it 9999999999") +from tools.postprocessK.kperm.flow import getKeff - pbc=float(Nz) +def PetscP(datadir, ref, k, saveres): - K,R = firstL(K,R,k,pbc) - r=(k.shape[1]-2)*(k.shape[2]-2)*nnz2 #start row - K,R =lastL(K,R,k2,r) + ref = 1 + rank = 0 + pn = 1 + t0 = time.time() - k2=0 + pcomm = PETSc.COMM_SELF + if k.shape[2] == 1: + refz = 1 + else: + refz = ref - K.assemble() - R.assemble() + nz, ny, nx = k.shape[0] * ref, k.shape[1] * ref, k.shape[2] * refz + n = nx * ny * nz + K = PETSc.Mat().create(comm=pcomm) + K.setType("seqaij") + K.setSizes(((n, None), (n, None))) # Aca igual que lo que usas arriba + K.setPreallocationNNZ(nnz=(7, 4)) # Idem anterior + K.setUp() + R = PETSc.Vec().createSeq((n, None), comm=pcomm) # PETSc.COMM_WORLD + R.setUp() - ksp = PETSc.KSP() - ksp.create(comm=pcomm) - ksp.setFromOptions() - P = R.copy() - ksp.setType(PETSc.KSP.Type.CG) - pc = PETSc.PC() + k2, Nz, nnz2 = getKref(k, 1, 2, ref) + k, Nz, nnz = getKref(k, 0, 2, ref) - pc.create(comm=pcomm) - pc.setType(PETSc.PC.Type.JACOBI) - ksp.setPC(pc) - ksp.setOperators(K) - ksp.setUp() - t1=time.time() - ksp.solve(R, P) - t2=time.time() - p=P.getArray().reshape(nz,ny,nx) + pbc = float(Nz) - if rank==0: - keff,Q=getKeff(p,k[1:-1,1:-1,1:-1],pbc,Nz) - print(keff,ref,nx,ny,nz) - return keff + K, R = firstL(K, R, k, pbc) + r = (k.shape[1] - 2) * (k.shape[2] - 2) * nnz2 # start row + K, R = lastL(K, R, k2, r) + k2 = 0 - return + K.assemble() + R.assemble() + ksp = PETSc.KSP() + ksp.create(comm=pcomm) + ksp.setFromOptions() + P = R.copy() + ksp.setType(PETSc.KSP.Type.CG) + pc = PETSc.PC() + pc.create(comm=pcomm) + pc.setType(PETSc.PC.Type.JACOBI) + ksp.setPC(pc) + ksp.setOperators(K) + ksp.setUp() + t1 = time.time() + ksp.solve(R, P) + t2 = time.time() + p = P.getArray().reshape(nz, ny, nx) -#Ver: A posteriori error estimates and adaptive solvers for porous media flows (Martin Vohralik) + if rank == 0: + keff, Q = getKeff(p, k[1:-1, 1:-1, 1:-1], pbc, Nz) + print(keff, ref, nx, ny, nz) + return keff + return +# Ver: A posteriori error estimates and adaptive solvers for porous media flows (Martin Vohralik) diff --git a/tools/postprocessK/kperm/Ndar1PBack.py b/tools/postprocessK/kperm/Ndar1PBack.py index ffaebc3..d8a43d5 100755 --- a/tools/postprocessK/kperm/Ndar1PBack.py +++ b/tools/postprocessK/kperm/Ndar1PBack.py @@ -1,135 +1,123 @@ import numpy as np -#import petsc4py + +# import petsc4py import math import time -#from mpi4py import MPI + +# from mpi4py import MPI from tools.postprocessK.kperm.computeFlows import * from petsc4py import PETSc -#petsc4py.init('-ksp_max_it 9999999999',comm=PETSc.COMM_SELF) -from tools.postprocessK.flow import getKeff - - - - - - -def PetscP(datadir,ref,k,saveres): - - #datadir='./data/'+str(job)+'/' - - - #comm=MPI.COMM_WORLD - #rank=comm.Get_rank() - ''' - size=comm.Get_size() - print(rank,size) - pcomm = MPI.COMM_WORLD.Split(color=rank, key=rank) - #print(new_comm.Get_rank()) - - #pcomm=comm.Create(newgroup) - print('entro') - print pcomm.Get_rank() - - print pcomm.Get_size() - pcomm=comm - - - rank=pcomm.rank - pn=pcomm.size - #PETSc.COMM_WORLD.PetscSubcommCreate(pcomm,PetscSubcomm *psubcomm) - print(rank,pn) - ''' - #Optpetsc = PETSc.Options() - rank=0 - pn=1 - t0=time.time() - - #comm=MPI.Comm.Create() - if k.shape[2]==1: - refz=1 - else: - refz=ref - - nz, ny, nx=k.shape[0]*ref,k.shape[1]*ref,k.shape[2]*refz - n=nx*ny*nz - - print('algo') - - K = PETSc.Mat().create(comm=PETSc.COMM_SELF) - print('algo2') - K.setType('seqaij') - print('algo3') - K.setSizes(((n,None),(n,None))) # Aca igual que lo que usas arriba - K.setPreallocationNNZ(nnz=(7,4)) # Idem anterior - - - #K = PETSc.Mat('seqaij', m=n,n=n,nz=7,comm=PETSc.COMM_WORLD) - #K = PETSc.Mat('aij', ((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) - #K = PETSc.Mat().createAIJ(((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) - #K = PETSc.Mat().createSeqAIJ(((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) - #K.setPreallocationNNZ(nnz=(7,4)) - print('ksetup') - #K.MatCreateSeqAIJ() - #K=PETSc.Mat().MatCreate(PETSc.COMM_WORLD) - - #K = PETSc.Mat().createAIJ(((n,None),(n,None)), nnz=(7,4),comm=pcomm) - - K.setUp() - print('entro2') - R = PETSc.Vec().createSeq((n,None),comm=PETSc.COMM_SELF) #PETSc.COMM_WORLD - R.setUp() - print('entro2') - k2, Nz, nnz2=getKref(k,1,2,ref) - k, Nz, nnz=getKref(k,0,2,ref) - - - pbc=float(Nz) - #print('entro3') - - K,R = firstL(K,R,k,pbc) - r=(k.shape[1]-2)*(k.shape[2]-2)*nnz2 #start row - K,R =lastL(K,R,k2,r) - - k2=0 - - K.assemble() - R.assemble() - - - print('entro3') - ksp = PETSc.KSP() - ksp.create(comm=PETSc.COMM_SELF) - ksp.setFromOptions() - print('entro4') - P = R.copy() - ksp.setType(PETSc.KSP.Type.CG) - pc = PETSc.PC() - - pc.create(comm=PETSc.COMM_SELF) - print('entro4') - pc.setType(PETSc.PC.Type.JACOBI) - ksp.setPC(pc) - ksp.setOperators(K) - ksp.setUp() - t1=time.time() - ksp.solve(R, P) - t2=time.time() - p=P.getArray().reshape(nz,ny,nx) - - - - if rank==0: - keff,Q=getKeff(p,k[1:-1,1:-1,1:-1],pbc,Nz) - - return keff - - - return - - - -#Ver: A posteriori error estimates and adaptive solvers for porous media flows (Martin Vohralik) +# petsc4py.init('-ksp_max_it 9999999999',comm=PETSc.COMM_SELF) +from tools.postprocessK.flow import getKeff +def PetscP(datadir, ref, k, saveres): + + # datadir='./data/'+str(job)+'/' + + # comm=MPI.COMM_WORLD + # rank=comm.Get_rank() + """ + size=comm.Get_size() + print(rank,size) + pcomm = MPI.COMM_WORLD.Split(color=rank, key=rank) + #print(new_comm.Get_rank()) + + #pcomm=comm.Create(newgroup) + print('entro') + print pcomm.Get_rank() + + print pcomm.Get_size() + pcomm=comm + + + rank=pcomm.rank + pn=pcomm.size + #PETSc.COMM_WORLD.PetscSubcommCreate(pcomm,PetscSubcomm *psubcomm) + print(rank,pn) + """ + # Optpetsc = PETSc.Options() + rank = 0 + pn = 1 + t0 = time.time() + + # comm=MPI.Comm.Create() + if k.shape[2] == 1: + refz = 1 + else: + refz = ref + + nz, ny, nx = k.shape[0] * ref, k.shape[1] * ref, k.shape[2] * refz + n = nx * ny * nz + + print("algo") + + K = PETSc.Mat().create(comm=PETSc.COMM_SELF) + print("algo2") + K.setType("seqaij") + print("algo3") + K.setSizes(((n, None), (n, None))) # Aca igual que lo que usas arriba + K.setPreallocationNNZ(nnz=(7, 4)) # Idem anterior + + # K = PETSc.Mat('seqaij', m=n,n=n,nz=7,comm=PETSc.COMM_WORLD) + # K = PETSc.Mat('aij', ((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) + # K = PETSc.Mat().createAIJ(((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) + # K = PETSc.Mat().createSeqAIJ(((n,None),(n,None)), nnz=(7,4),comm=PETSc.COMM_WORLD) + # K.setPreallocationNNZ(nnz=(7,4)) + print("ksetup") + # K.MatCreateSeqAIJ() + # K=PETSc.Mat().MatCreate(PETSc.COMM_WORLD) + + # K = PETSc.Mat().createAIJ(((n,None),(n,None)), nnz=(7,4),comm=pcomm) + + K.setUp() + print("entro2") + R = PETSc.Vec().createSeq((n, None), comm=PETSc.COMM_SELF) # PETSc.COMM_WORLD + R.setUp() + print("entro2") + k2, Nz, nnz2 = getKref(k, 1, 2, ref) + k, Nz, nnz = getKref(k, 0, 2, ref) + + pbc = float(Nz) + # print('entro3') + + K, R = firstL(K, R, k, pbc) + r = (k.shape[1] - 2) * (k.shape[2] - 2) * nnz2 # start row + K, R = lastL(K, R, k2, r) + + k2 = 0 + + K.assemble() + R.assemble() + + print("entro3") + ksp = PETSc.KSP() + ksp.create(comm=PETSc.COMM_SELF) + ksp.setFromOptions() + print("entro4") + P = R.copy() + ksp.setType(PETSc.KSP.Type.CG) + pc = PETSc.PC() + + pc.create(comm=PETSc.COMM_SELF) + print("entro4") + pc.setType(PETSc.PC.Type.JACOBI) + ksp.setPC(pc) + ksp.setOperators(K) + ksp.setUp() + t1 = time.time() + ksp.solve(R, P) + t2 = time.time() + p = P.getArray().reshape(nz, ny, nx) + + if rank == 0: + keff, Q = getKeff(p, k[1:-1, 1:-1, 1:-1], pbc, Nz) + + return keff + + return + + +# Ver: A posteriori error estimates and adaptive solvers for porous media flows (Martin Vohralik) diff --git a/tools/postprocessK/kperm/computeFlows.py b/tools/postprocessK/kperm/computeFlows.py index 0164842..9f4666c 100755 --- a/tools/postprocessK/kperm/computeFlows.py +++ b/tools/postprocessK/kperm/computeFlows.py @@ -2,50 +2,72 @@ import numpy as np import math -def getKref(k,rank,pn,ref): - - - Nz = k.shape[0] - nz = Nz//pn - if ref==1: - return getK(k,rank,pn) - - - if (rank>0) and (rank 0) and (rank < pn - 1): + + k = k[rank * nz - 1 : (rank + 1) * nz + 1, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[(ref - 1) : -(ref - 1), :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz, ny + 2, nx + 2)) + ki[:, 1:-1, 1:-1] = k + nnz = nz + if rank == 0: + k = k[: (rank + 1) * nz + 1, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[: -(ref - 1), :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz + 1, ny + 2, nx + 2)) + ki[1:, 1:-1, 1:-1] = k + ki[0, :, :] = ki[1, :, :] + nnz = nz + if rank == (pn - 1): + k = k[rank * nz - 1 :, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[(ref - 1) :, :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz + 1, ny + 2, nx + 2)) + ki[:-1, 1:-1, 1:-1] = k + ki[-1, :, :] = ki[-2, :, :] + nnz = (Nz // pn) * ref + return ki, Nz * ref, nnz + + +def getK(k, rank, pn): + + # k=np.load(kfile) + # nn=int(np.cbrt(k.shape[0])) + # k=k.reshape((nn,nn,nn)) + + Nz, Ny, Nx = k.shape[0], k.shape[1], k.shape[2] + nz = Nz // pn + if rank == pn - 1: + nnz = Nz - (pn - 1) * nz + ki = np.zeros((nnz + 2, Ny + 2, Nx + 2)) + else: + nnz = nz + ki = np.zeros((nz + 2, Ny + 2, Nx + 2)) + if (rank > 0) and (rank < pn - 1): + ki[:, 1:-1, 1:-1] = k[rank * nz - 1 : (rank + 1) * nz + 1, :, :] + if rank == 0: + ki[1:, 1:-1, 1:-1] = k[: (rank + 1) * nz + 1, :, :] + ki[0, :, :] = ki[1, :, :] + if rank == (pn - 1): + ki[:-1, 1:-1, 1:-1] = k[rank * nz - 1 :, :, :] + ki[-1, :, :] = ki[-2, :, :] + return ki, Nz, nz + + +""" def getK(k,rank,pn): #k=np.load(kfile) @@ -69,188 +91,307 @@ def getK(k,rank,pn): ki[:-1,1:-1,1:-1]=k[rank*nz-1:,:,:] ki[-1,:,:]=ki[-2,:,:] return ki, Nz, nz -''' -def getK(k,rank,pn): +""" - #k=np.load(kfile) - #nn=int(np.cbrt(k.shape[0])) - #k=k.reshape((nn,nn,nn)) - - Nz, Ny,Nx=k.shape[0],k.shape[1],k.shape[2] - nz=Nz//pn - if rank==pn-1: - nnz= Nz-(pn-1)*nz - ki=np.zeros((nnz+2,Ny+2,Nx+2)) - else: - nnz=nz - ki=np.zeros((nz+2,Ny+2,Nx+2)) - if (rank>0) and (rank0) and (rank0) and (rank 0) and (rank < pn - 1): + K, R = centL(K, R, k, r) + k = 0 + if rank == (pn - 1): + K, R = lastL(K, R, k, r) + k = 0 + + K.assemble() + R.assemble() + + ksp = PETSc.KSP() + ksp.create(comm=pcomm) + ksp.setTolerances(rtol=Rtol, atol=1.0e-100, max_it=999999999) + ksp.setFromOptions() + P = R.copy() + ksp.setType(PETSc.KSP.Type.CG) + pc = PETSc.PC() + pc.create(comm=pcomm) + pc.setType(PETSc.PC.Type.JACOBI) + ksp.setPC(pc) + ksp.setOperators(K) + ksp.setUp() + t1 = time.time() + ksp.solve(R, P) + t2 = time.time() + p = P.getArray().reshape(nz, ny, nx) + + if rank == 0: + keff, Q = getKeff(p, k[1:-1, 1:-1, 1:-1], pbc, Nz) + if saveres == True: + + for i in range(1, pn): + from mpi4py import MPI + + comm = MPI.COMM_WORLD + pi = comm.recv(source=i) + p = np.append(p, pi, axis=0) + + np.save(datadir + "P", p) + f = open(datadir + "RunTimes.out", "a") + f.write("ref: " + str(ref) + "\n") + f.write("Matrix creation: " + str(t1 - t0) + "\n") + f.write("Solver: " + str(t2 - t1) + "\n") + f.write("Keff: " + str(keff) + "\n") + f.write("N_cores: " + str(pn) + "\n") + f.close() + try: + res = np.loadtxt(datadir + "SolverRes.txt") + res = np.append(res, np.array([keff, ref, t2 - t0, pn])) + except: + res = np.array([keff, ref, t2 - t0, pn]) + np.savetxt( + datadir + "SolverRes.txt", res, header="Keff, ref, Runtime, N_cores" + ) + print(datadir[-3:], " keff= " + str(keff), " rtime= " + str(t2 - t0)) + return keff + + else: + if saveres == True: + from mpi4py import MPI + + comm = MPI.COMM_WORLD + comm.send(p, dest=0) + + return + + +# Ver: A posteriori error estimates and adaptive solvers for porous media flows (Martin Vohralik) + + +ddir = "./test/0/" +ref = 1 icomm = MPI.Comm.Get_parent() -print('aca') -PetscP(ddir,ref,'0',True,0.000001,1) -#icomm = MPI.Comm.Get_parent() +print("aca") +PetscP(ddir, ref, "0", True, 0.000001, 1) +# icomm = MPI.Comm.Get_parent() icomm.Disconnect() - - - - - diff --git a/tools/solver/comp_Kperm_scale.py b/tools/solver/comp_Kperm_scale.py index 89fa743..4efcd4c 100755 --- a/tools/solver/comp_Kperm_scale.py +++ b/tools/solver/comp_Kperm_scale.py @@ -3,106 +3,105 @@ import os import time from tools.solver.Ndar import PetscP -def comp_kperm_sub(parser,rundir,nr): +def comp_kperm_sub(parser, rundir, nr): - k=np.load(rundir+'k.npy') + k = np.load(rundir + "k.npy") - ref=int(parser.get('Solver',"ref")) + ref = int(parser.get("Solver", "ref")) + t0 = time.time() - t0=time.time() + S_min_post = int(parser.get("K-Postprocess", "MinBlockSize")) + nimax = 2 ** int(parser.get("K-Postprocess", "Max_sample_size")) - S_min_post = int(parser.get('K-Postprocess','MinBlockSize')) - nimax =2** int(parser.get('K-Postprocess','Max_sample_size')) - + S_min_post = S_min_post * ref - S_min_post=S_min_post*ref + if S_min_post == 0: + sx = 2 # k.shape[0] + else: + sx = get_min_nbl(k, nimax, nr, S_min_post) - if S_min_post==0: - sx=2 #k.shape[0] - else: - sx = get_min_nbl(k,nimax,nr,S_min_post) + if sx == 1: + sx = 2 - if sx==1: - sx=2 + tkperm = getKpost(k, sx, rundir, ref) - tkperm=getKpost(k, sx,rundir,ref) + ttotal = time.time() - t0 + return - ttotal=time.time()-t0 - return +def getKpost(kf, sx, rundir, ref): + ex = int(np.log2(kf.shape[0])) + esx = int(np.log2(sx)) + scales = 2 ** np.arange(esx, ex) + datadir = rundir + "KpostProcess/" + try: + os.makedirs(datadir) + except: + nada = 0 + tkperm = np.zeros((scales.shape[0])) -def getKpost(kf, sx,rundir,ref): + for il in range(scales.shape[0]): + l = scales[il] + nblx, nbly, nblz = kf.shape[0] // l, kf.shape[1] // l, kf.shape[2] // l + sx, sy, sz = l, l, l + if kf.shape[2] == 1: + nblz = 1 + sz = 1 + if l == 2: + refDeg = 2 + else: + refDeg = ref + tkperm[il] = time.time() + Kperm = np.zeros((nblx, nbly, nblz)) + try: + Kperm = np.load(datadir + "Kperm" + str(l // ref) + ".npy") + except: + for i in range(nblx): + for j in range(nbly): + for k in range(nblz): - ex=int(np.log2(kf.shape[0])) - esx=int(np.log2(sx)) + Kperm[i, j, k] = PetscP( + "", + refDeg, + kf[ + i * sx : (i + 1) * sx, + j * sy : (j + 1) * sy, + k * sz : (k + 1) * sz, + ], + False, + 1e-4, + 0, + ) - scales=2**np.arange(esx,ex) - datadir=rundir+'KpostProcess/' - try: - os.makedirs(datadir) - except: - nada=0 + tkperm[il] = time.time() - tkperm[il] + np.save(datadir + "Kperm" + str(sx) + ".npy", Kperm) - tkperm=np.zeros((scales.shape[0])) + np.savetxt(rundir + "tkperm_sub.txt", tkperm) - for il in range(scales.shape[0]): - l=scales[il] - nblx, nbly, nblz = kf.shape[0]//l, kf.shape[1]//l, kf.shape[2]//l - sx,sy,sz=l,l,l - if kf.shape[2]==1: - nblz=1 - sz=1 - if l==2: - refDeg=2 - else: - refDeg=ref + return tkperm - tkperm[il]=time.time() - Kperm = np.zeros((nblx,nbly,nblz)) - try: - Kperm=np.load(datadir+'Kperm'+str(l//ref)+'.npy') - - except: - for i in range(nblx): - for j in range(nbly): - for k in range(nblz): - - Kperm[i,j,k]=PetscP('',refDeg,kf[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],False,1e-4,0) +def get_min_nbl(kc, nimax, nr, smin): - - tkperm[il]= time.time()-tkperm[il] - np.save(datadir+'Kperm'+str(sx)+'.npy',Kperm) - - - - np.savetxt(rundir+'tkperm_sub.txt',tkperm) - - return tkperm - - -def get_min_nbl(kc,nimax,nr,smin): - - if kc.shape[2]==1: - dim=2.0 - else: - dim=3.0 - if nr>0: - y=(1/dim)*np.log2(nr*kc.size/(nimax*(smin**dim))) - else: - y=0 - y=int(y) - s=int((2**y) * smin) - if s 0: + y = (1 / dim) * np.log2(nr * kc.size / (nimax * (smin ** dim))) + else: + y = 0 + y = int(y) + s = int((2 ** y) * smin) + if s < smin: + s = smin + return s diff --git a/tools/solver/computeFlows.py b/tools/solver/computeFlows.py index 281772b..95d2d79 100755 --- a/tools/solver/computeFlows.py +++ b/tools/solver/computeFlows.py @@ -2,51 +2,73 @@ import numpy as np import math -def getKref(k,rank,pn,ref): - - - Nz = k.shape[0] - nz = Nz//pn - if ref==1: - return getK(k,rank,pn) - - - if (rank>0) and (rank 0) and (rank < pn - 1): + + k = k[rank * nz - 1 : (rank + 1) * nz + 1, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[(ref - 1) : -(ref - 1), :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz, ny + 2, nx + 2)) + ki[:, 1:-1, 1:-1] = k + # print(ki.shape) + nnz = nz - 2 + if rank == 0: + k = k[: (rank + 1) * nz + 1, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[: -(ref - 1), :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz + 1, ny + 2, nx + 2)) + ki[1:, 1:-1, 1:-1] = k + ki[0, :, :] = ki[1, :, :] + nnz = nz + if rank == (pn - 1): + k = k[rank * nz - 1 :, :, :] + k = refinaPy(k, ref) + if ref != 1: + k = k[(ref - 1) :, :, :] + nz, ny, nx = k.shape[0], k.shape[1], k.shape[2] + ki = np.zeros((nz + 1, ny + 2, nx + 2)) + ki[:-1, 1:-1, 1:-1] = k + ki[-1, :, :] = ki[-2, :, :] + nnz = (Nz // pn) * ref + return ki, Nz * ref, nnz + + +def getK(k, rank, pn): + + # k=np.load(kfile) + # nn=int(np.cbrt(k.shape[0])) + # k=k.reshape((nn,nn,nn)) + + Nz, Ny, Nx = k.shape[0], k.shape[1], k.shape[2] + nz = Nz // pn + if rank == pn - 1: + nnz = Nz - (pn - 1) * nz + ki = np.zeros((nnz + 2, Ny + 2, Nx + 2)) + else: + nnz = nz + ki = np.zeros((nz + 2, Ny + 2, Nx + 2)) + if (rank > 0) and (rank < pn - 1): + ki[:, 1:-1, 1:-1] = k[rank * nz - 1 : (rank + 1) * nz + 1, :, :] + if rank == 0: + ki[1:, 1:-1, 1:-1] = k[: (rank + 1) * nz + 1, :, :] + ki[0, :, :] = ki[1, :, :] + if rank == (pn - 1): + ki[:-1, 1:-1, 1:-1] = k[rank * nz - 1 :, :, :] + ki[-1, :, :] = ki[-2, :, :] + return ki, Nz, nz + + +""" def getK(k,rank,pn): #k=np.load(kfile) @@ -70,187 +92,306 @@ def getK(k,rank,pn): ki[:-1,1:-1,1:-1]=k[rank*nz-1:,:,:] ki[-1,:,:]=ki[-2,:,:] return ki, Nz, nz -''' -def getK(k,rank,pn): +""" - #k=np.load(kfile) - #nn=int(np.cbrt(k.shape[0])) - #k=k.reshape((nn,nn,nn)) - - Nz, Ny,Nx=k.shape[0],k.shape[1],k.shape[2] - nz=Nz//pn - if rank==pn-1: - nnz= Nz-(pn-1)*nz - ki=np.zeros((nnz+2,Ny+2,Nx+2)) - else: - nnz=nz - ki=np.zeros((nz+2,Ny+2,Nx+2)) - if (rank>0) and (rank)$' + label = r"$\log_{10}(vx/)$" - folder=j*50+i + folder = j * 50 + i - V=np.load(rdir+str(folder)+'/V.npy')[0][:,:,0] - perco=np.load(rdir+str(folder)+'/ConnectivityMetrics/1024.npy',allow_pickle=True).item()['spanning'][0,0,0] - V=np.log10(np.abs(V)) #/np.mean(np.abs(V))) - leg='p = '+str(ps[i])[:4]+'% ('+str(perco)+')' - plot_hist(V,rdir+str(folder)+'/HisTabsV',log,label,-.8,.5,leg) - plotK_imshow(V[512:1536,512:1536],rdir+str(i)+'/V',log,label,-4,1) - plt.legend(loc='upper left') - plt.savefig(rdir+str(folder)+'VelHistogramB.png') - plt.close() -''' + V = np.load(rdir + str(folder) + "/V.npy")[0][:, :, 0] + perco = np.load( + rdir + str(folder) + "/ConnectivityMetrics/1024.npy", allow_pickle=True + ).item()["spanning"][0, 0, 0] + V = np.log10(np.abs(V)) # /np.mean(np.abs(V))) + leg = "p = " + str(ps[i])[:4] + "% (" + str(perco) + ")" + plot_hist(V, rdir + str(folder) + "/HisTabsV", log, label, -0.8, 0.5, leg) + plotK_imshow(V[512:1536, 512:1536], rdir + str(i) + "/V", log, label, -4, 1) + plt.legend(loc="upper left") + plt.savefig(rdir + str(folder) + "VelHistogramB.png") + plt.close() +""" label=r'$\log_{10}(|v_x|/<|v_x|>)$' V=np.load(rdir+str(i)+'/V.npy')[0][:,:,0] @@ -80,6 +83,4 @@ plotK_imshow(V[1024:2048,512:1024],rdir+str(i)+'/Vy',log,label,0,1) -''' - - +""" diff --git a/utilities/plot_lc.py b/utilities/plot_lc.py index e407dd5..df2a937 100755 --- a/utilities/plot_lc.py +++ b/utilities/plot_lc.py @@ -2,22 +2,20 @@ import numpy as np import matplotlib.pyplot as plt +rdir = "./lc_vslcbin/" +nps = 41 +ps = np.linspace(0.1, 0.5, nps) -rdir='./lc_vslcbin/' -nps=41 -ps=np.linspace(.1,.5,nps) - -clabels=['Intermediate','high','low'] +clabels = ["Intermediate", "high", "low"] for con in range(1): - keff=np.zeros(nps) - for ip in range(nps): - folder=con*nps+ip - keff[ip]=np.loadtxt(rdir+str(folder)+'/lc.txt')[2] + keff = np.zeros(nps) + for ip in range(nps): + folder = con * nps + ip + keff[ip] = np.loadtxt(rdir + str(folder) + "/lc.txt")[2] - plt.plot(ps,keff,label=clabels[con]) + plt.plot(ps, keff, label=clabels[con]) plt.legend() plt.grid() -plt.savefig('lc2.png') - +plt.savefig("lc2.png") diff --git a/utilities/plot_val_Kpost.py b/utilities/plot_val_Kpost.py index 5d67b34..4a22580 100755 --- a/utilities/plot_val_Kpost.py +++ b/utilities/plot_val_Kpost.py @@ -2,32 +2,60 @@ import numpy as np import matplotlib.pyplot as plt +rdir = "./data/" -rdir='./data/' +clabels = [r"$K_{perm}$", r"$K_{diss}$", r"$K_{average}$"] +cases = [ + r"$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$", + r"$Lognormal \ \sigma^{2}_{\log(k)} = 7$", + r"$Binary p = 0.2; k+/k- = 10^4$", +] -clabels=[r'$K_{perm}$',r'$K_{diss}$',r'$K_{average}$'] -cases=[r'$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$',r'$Lognormal \ \sigma^{2}_{\log(k)} = 7$', r'$Binary p = 0.2; k+/k- = 10^4$'] - -scales=np.array([4,8,16,32,64,128,256,512]) -lcs=[16,16,8] +scales = np.array([4, 8, 16, 32, 64, 128, 256, 512]) +lcs = [16, 16, 8] for i in range(3): - kpost=np.zeros((len(scales),3)) - for scale in range(len(scales)): - kpost[scale,0]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/kperm/'+str(scales[scale])+'.npy')))) - kpost[scale,1]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/KpostProcess/Kd'+str(scales[scale])+'.npy')))) - kpost[scale,2]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/KpostProcess/Kv'+str(scales[scale])+'.npy')))) - plt.semilogx(scales/512.0,kpost[:,0],label=clabels[0],marker='x') - plt.semilogx(scales/512.0,kpost[:,1],label=clabels[1],marker='s') - plt.semilogx(scales/512.0,kpost[:,2],label=clabels[2],marker='^') - plt.vlines(lcs[i]/512.0,kpost[:,0].min(),kpost[:,0].max(),label=r'$lc = $'+str(lcs[i])) - plt.xlabel(r'$\lambda / L$') - plt.ylabel(r'$_G$') - plt.legend() - plt.grid() - plt.title(cases[i]) - plt.tight_layout() - plt.savefig(rdir+str(i)+'/Kpost_mean.png') - plt.close() + kpost = np.zeros((len(scales), 3)) + for scale in range(len(scales)): + kpost[scale, 0] = np.exp( + np.nanmean( + np.log(np.load(rdir + str(i) + "/kperm/" + str(scales[scale]) + ".npy")) + ) + ) + kpost[scale, 1] = np.exp( + np.nanmean( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kd" + str(scales[scale]) + ".npy" + ) + ) + ) + ) + kpost[scale, 2] = np.exp( + np.nanmean( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kv" + str(scales[scale]) + ".npy" + ) + ) + ) + ) + plt.semilogx(scales / 512.0, kpost[:, 0], label=clabels[0], marker="x") + plt.semilogx(scales / 512.0, kpost[:, 1], label=clabels[1], marker="s") + plt.semilogx(scales / 512.0, kpost[:, 2], label=clabels[2], marker="^") + plt.vlines( + lcs[i] / 512.0, + kpost[:, 0].min(), + kpost[:, 0].max(), + label=r"$lc = $" + str(lcs[i]), + ) + plt.xlabel(r"$\lambda / L$") + plt.ylabel(r"$_G$") + plt.legend() + plt.grid() + plt.title(cases[i]) + plt.tight_layout() + plt.savefig(rdir + str(i) + "/Kpost_mean.png") + plt.close() diff --git a/utilities/plot_val_Kpost_varianve.py b/utilities/plot_val_Kpost_varianve.py index 3f3c4dc..d7035e0 100755 --- a/utilities/plot_val_Kpost_varianve.py +++ b/utilities/plot_val_Kpost_varianve.py @@ -3,42 +3,90 @@ import matplotlib.pyplot as plt from Var_analytical import * -rdir='./data/' +rdir = "./data/" -clabels=[r'$K_{perm}$',r'$K_{diss}$',r'$K_{average}$',r'$K_{1/3}$','analitycal Gaussian cov'] -cases=[r'$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$',r'$Lognormal \ \sigma^{2}_{\log(k)} = 7$', r'$Binary p = 0.2; k+/k- = 10^4$'] +clabels = [ + r"$K_{perm}$", + r"$K_{diss}$", + r"$K_{average}$", + r"$K_{1/3}$", + "analitycal Gaussian cov", +] +cases = [ + r"$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$", + r"$Lognormal \ \sigma^{2}_{\log(k)} = 7$", + r"$Binary p = 0.2; k+/k- = 10^4$", +] -scales=np.array([4,8,16,32,64,128]) -variances=[0.1,7,13.572859162824695] -x=scales/512.0 -lcs=[16,16,8] -va=VarLgauss(16/2.45398,scales,3) +scales = np.array([4, 8, 16, 32, 64, 128]) +variances = [0.1, 7, 13.572859162824695] +x = scales / 512.0 +lcs = [16, 16, 8] +va = VarLgauss(16 / 2.45398, scales, 3) for i in range(3): + kpost = np.zeros((len(scales), 4)) + for scale in range(len(scales)): + kpost[scale, 0] = ( + np.nanvar( + np.log(np.load(rdir + str(i) + "/kperm/" + str(scales[scale]) + ".npy")) + ) + / variances[i] + ) + kpost[scale, 1] = ( + np.nanvar( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kd" + str(scales[scale]) + ".npy" + ) + ) + ) + / variances[i] + ) + kpost[scale, 2] = ( + np.nanvar( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kv" + str(scales[scale]) + ".npy" + ) + ) + ) + / variances[i] + ) + kpost[scale, 3] = ( + np.nanvar( + np.log( + np.load( + rdir + + str(i) + + "/KpostProcess/Kpo" + + str(scales[scale]) + + ".npy" + ) + ) + ) + / variances[i] + ) + plt.loglog(x, (x ** 3) * kpost[:, 0], label=clabels[0], marker="x") + plt.loglog(x, (x ** 3) * kpost[:, 1], label=clabels[1], marker="s") + plt.loglog(x, (x ** 3) * kpost[:, 2], label=clabels[2], marker="^") + plt.loglog(x, (x ** 3) * kpost[:, 3], label=clabels[3], marker="o") + if i == 0 or i == 1: + plt.loglog(x, (x ** 3) * va, label=clabels[4], marker="", linestyle="--") - - kpost=np.zeros((len(scales),4)) - for scale in range(len(scales)): - kpost[scale,0]=np.nanvar(np.log(np.load(rdir+str(i)+'/kperm/'+str(scales[scale])+'.npy')))/variances[i] - kpost[scale,1]=np.nanvar(np.log(np.load(rdir+str(i)+'/KpostProcess/Kd'+str(scales[scale])+'.npy')))/variances[i] - kpost[scale,2]=np.nanvar(np.log(np.load(rdir+str(i)+'/KpostProcess/Kv'+str(scales[scale])+'.npy')))/variances[i] - kpost[scale,3]=np.nanvar(np.log(np.load(rdir+str(i)+'/KpostProcess/Kpo'+str(scales[scale])+'.npy')))/variances[i] - plt.loglog(x,(x**3)*kpost[:,0],label=clabels[0],marker='x') - plt.loglog(x,(x**3)*kpost[:,1],label=clabels[1],marker='s') - plt.loglog(x,(x**3)*kpost[:,2],label=clabels[2],marker='^') - plt.loglog(x,(x**3)*kpost[:,3],label=clabels[3],marker='o') - if i==0 or i==1: - plt.loglog(x,(x**3)*va,label=clabels[4],marker='',linestyle='--') - - plt.vlines(lcs[i]/512.0,((x**3)*kpost[:,0]).min(),((x**3)*kpost[:,0]).max(),label=r'$lc = $'+str(lcs[i])) - plt.xlabel(r'$\lambda / L$') - plt.ylabel(r'$(\lambda / L)^3 \sigma^{2}_{\log(K_{eff})} / \sigma^{2}_{\log(k)}$') - plt.legend() - plt.grid() - plt.title(cases[i]) - plt.tight_layout() - plt.savefig(rdir+str(i)+'/Kpost_var.png') - plt.close() - + plt.vlines( + lcs[i] / 512.0, + ((x ** 3) * kpost[:, 0]).min(), + ((x ** 3) * kpost[:, 0]).max(), + label=r"$lc = $" + str(lcs[i]), + ) + plt.xlabel(r"$\lambda / L$") + plt.ylabel(r"$(\lambda / L)^3 \sigma^{2}_{\log(K_{eff})} / \sigma^{2}_{\log(k)}$") + plt.legend() + plt.grid() + plt.title(cases[i]) + plt.tight_layout() + plt.savefig(rdir + str(i) + "/Kpost_var.png") + plt.close() diff --git a/utilities/plot_val_Kpost_with_power (copy).py b/utilities/plot_val_Kpost_with_power (copy).py index 91f32a5..391046b 100755 --- a/utilities/plot_val_Kpost_with_power (copy).py +++ b/utilities/plot_val_Kpost_with_power (copy).py @@ -2,43 +2,64 @@ import numpy as np import matplotlib.pyplot as plt +rdir = "./data/" -rdir='./data/' +clabels = [r"$K_{perm}$", r"$K_{diss}$", r"$K_{average}$", r"$K_{1/3}$"] +names = ["Kperm", "Kdiss", "Kaverage", "Kpower"] +cases = [ + r"$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$", + r"$Lognormal \ \sigma^{2}_{\log(k)} = 7$", + r"$Binary p = 0.2; k+/k- = 10^4$", +] -clabels=[r'$K_{perm}$',r'$K_{diss}$',r'$K_{average}$',r'$K_{1/3}$'] -names=['Kperm','Kdiss','Kaverage','Kpower'] -cases=[r'$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$',r'$Lognormal \ \sigma^{2}_{\log(k)} = 7$', r'$Binary p = 0.2; k+/k- = 10^4$'] - -scales=np.array([4,8,16,32,64]) -lcs=[16,16,8] -est=3 -ranges=[(-0.5,0.5),(-5,5),(-4,4)] +scales = np.array([4, 8, 16, 32, 64]) +lcs = [16, 16, 8] +est = 3 +ranges = [(-0.5, 0.5), (-5, 5), (-4, 4)] for i in range(3): + for scale in range(len(scales)): + if est == 0: + keff = np.log( + np.load(rdir + str(i) + "/kperm/" + str(scales[scale]) + ".npy") + ) + if est == 1: + keff = np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kd" + str(scales[scale]) + ".npy" + ) + ) + + if est == 2: + keff = np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kv" + str(scales[scale]) + ".npy" + ) + ) + if est == 3: + keff = np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kpo" + str(scales[scale]) + ".npy" + ) + ) - for scale in range(len(scales)): - if est==0: - keff=np.log(np.load(rdir+str(i)+'/kperm/'+str(scales[scale])+'.npy')) - if est==1: - keff=np.log(np.load(rdir+str(i)+'/KpostProcess/Kd'+str(scales[scale])+'.npy')) - - if est==2: - keff=np.log(np.load(rdir+str(i)+'/KpostProcess/Kv'+str(scales[scale])+'.npy')) - if est==3: - keff=np.log(np.load(rdir+str(i)+'/KpostProcess/Kpo'+str(scales[scale])+'.npy')) - - - plt.hist(keff.reshape(-1),label=r'$\lambda = $'+' ' +str(scales[scale]),density=True,histtype='step',range=ranges[i]) - #plt.semilogx(scales/512.0,kpost[:,1],label=clabels[1],marker='s') - #plt.semilogx(scales/512.0,kpost[:,2],label=clabels[2],marker='^') - #plt.semilogx(scales/512.0,kpost[:,3],label=clabels[3],marker='o') - #plt.vlines(lcs[i]/512.0,kpost[:,0].min(),kpost[:,0].max(),label=r'$lc = $'+str(lcs[i])) - plt.xlabel(r'$\log(K_{eff})$') - plt.ylabel(r'$P(K_{eff})$') - plt.legend() - plt.grid() - plt.title(cases[i]+' '+str(names[est])) - plt.tight_layout() - plt.savefig(rdir+str(i)+'/Kpost_dist_scales_'+names[est]+'.png') - plt.close() + plt.hist( + keff.reshape(-1), + label=r"$\lambda = $" + " " + str(scales[scale]), + density=True, + histtype="step", + range=ranges[i], + ) + # plt.semilogx(scales/512.0,kpost[:,1],label=clabels[1],marker='s') + # plt.semilogx(scales/512.0,kpost[:,2],label=clabels[2],marker='^') + # plt.semilogx(scales/512.0,kpost[:,3],label=clabels[3],marker='o') + # plt.vlines(lcs[i]/512.0,kpost[:,0].min(),kpost[:,0].max(),label=r'$lc = $'+str(lcs[i])) + plt.xlabel(r"$\log(K_{eff})$") + plt.ylabel(r"$P(K_{eff})$") + plt.legend() + plt.grid() + plt.title(cases[i] + " " + str(names[est])) + plt.tight_layout() + plt.savefig(rdir + str(i) + "/Kpost_dist_scales_" + names[est] + ".png") + plt.close() diff --git a/utilities/plot_val_Kpost_with_power.py b/utilities/plot_val_Kpost_with_power.py index e47076b..033b690 100755 --- a/utilities/plot_val_Kpost_with_power.py +++ b/utilities/plot_val_Kpost_with_power.py @@ -2,36 +2,75 @@ import numpy as np import matplotlib.pyplot as plt +rdir = "./data/" -rdir='./data/' +clabels = [r"$K_{perm}$", r"$K_{diss}$", r"$K_{average}$", r"$K_{1/3}$"] +cases = [ + r"$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$", + r"$Lognormal \ \sigma^{2}_{\log(k)} = 7$", + r"$Binary p = 0.2; k+/k- = 10^4$", +] -clabels=[r'$K_{perm}$',r'$K_{diss}$',r'$K_{average}$',r'$K_{1/3}$'] -cases=[r'$Lognormal \ \sigma^{2}_{\log(k)} = 0.1$',r'$Lognormal \ \sigma^{2}_{\log(k)} = 7$', r'$Binary p = 0.2; k+/k- = 10^4$'] - -scales=np.array([4,8,16,32,64,128,256,512]) -lcs=[16,16,8] +scales = np.array([4, 8, 16, 32, 64, 128, 256, 512]) +lcs = [16, 16, 8] for i in range(3): - kpost=np.zeros((len(scales),4)) - - - for scale in range(len(scales)): - kpost[scale,0]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/kperm/'+str(scales[scale])+'.npy')))) - kpost[scale,1]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/KpostProcess/Kd'+str(scales[scale])+'.npy')))) - kpost[scale,2]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/KpostProcess/Kv'+str(scales[scale])+'.npy')))) - kpost[scale,3]=np.exp(np.nanmean(np.log(np.load(rdir+str(i)+'/KpostProcess/Kpo'+str(scales[scale])+'.npy')))) - plt.semilogx(scales/512.0,kpost[:,0],label=clabels[0],marker='x') - plt.semilogx(scales/512.0,kpost[:,1],label=clabels[1],marker='s') - plt.semilogx(scales/512.0,kpost[:,2],label=clabels[2],marker='^') - plt.semilogx(scales/512.0,kpost[:,3],label=clabels[3],marker='o') - plt.vlines(lcs[i]/512.0,kpost[:,0].min(),kpost[:,0].max(),label=r'$lc = $'+str(lcs[i])) - plt.xlabel(r'$\lambda / L$') - plt.ylabel(r'$_G$') - plt.legend() - plt.grid() - plt.title(cases[i]) - plt.tight_layout() - plt.savefig(rdir+str(i)+'/Kpost_mean.png') - plt.close() + kpost = np.zeros((len(scales), 4)) + + for scale in range(len(scales)): + kpost[scale, 0] = np.exp( + np.nanmean( + np.log(np.load(rdir + str(i) + "/kperm/" + str(scales[scale]) + ".npy")) + ) + ) + kpost[scale, 1] = np.exp( + np.nanmean( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kd" + str(scales[scale]) + ".npy" + ) + ) + ) + ) + kpost[scale, 2] = np.exp( + np.nanmean( + np.log( + np.load( + rdir + str(i) + "/KpostProcess/Kv" + str(scales[scale]) + ".npy" + ) + ) + ) + ) + kpost[scale, 3] = np.exp( + np.nanmean( + np.log( + np.load( + rdir + + str(i) + + "/KpostProcess/Kpo" + + str(scales[scale]) + + ".npy" + ) + ) + ) + ) + plt.semilogx(scales / 512.0, kpost[:, 0], label=clabels[0], marker="x") + plt.semilogx(scales / 512.0, kpost[:, 1], label=clabels[1], marker="s") + plt.semilogx(scales / 512.0, kpost[:, 2], label=clabels[2], marker="^") + plt.semilogx(scales / 512.0, kpost[:, 3], label=clabels[3], marker="o") + plt.vlines( + lcs[i] / 512.0, + kpost[:, 0].min(), + kpost[:, 0].max(), + label=r"$lc = $" + str(lcs[i]), + ) + plt.xlabel(r"$\lambda / L$") + plt.ylabel(r"$_G$") + plt.legend() + plt.grid() + plt.title(cases[i]) + plt.tight_layout() + plt.savefig(rdir + str(i) + "/Kpost_mean.png") + plt.close()