apps: Included impala implementation of md with short-range interactions

apps/impala/CMakeLists.txt

+cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
+find_package(AnyDSL-runtime REQUIRED)
+include(${ANYDSL_RUNTIME_CMAKE_DIR}/Runtime.cmake)
+include_directories(${ANYDSL_RUNTIME_DIR}/src)
+
+set(CLANG_FLAGS -O3 -march=native)
+set(IMPALA_FLAGS --log-level info)
+
+# examples based on the "mapping"
+anydsl_runtime_wrap(IMPALA_OBJECT_FILE
+    CLANG_FLAGS ${CLANG_FLAGS}
+    IMPALA_FLAGS ${IMPALA_FLAGS}
+    FILES intrinsics/intrinsics_cpu.impala common.impala datastructures/linked_cell/cpu.impala utilities/linked_cell/cpu.impala algorithm/linked_cell/cpu.impala initialization/test.impala potential/lennard_jones.impala integration/verlet.impala boundary/rigid_walls.impala time_integration/cpu.impala run/cpu.impala)
+
+waLBerla_add_executable ( NAME md_simulation
+    FILES md_simulation.cpp ${IMPALA_OBJECT_FILE} clib/fileIO.c clib/allocate.c
+    DEPENDS core blockforest pe domain_decomposition vtk)

apps/impala/algorithm/array/cpu.impala

+
+fn compute_force(P: ParticleSystem, force: fn(Particle, Particle, Constants) -> ()) -> ()
+{
+    let start = 0 as size_t;
+    let lut = P.lut();
+    for i in par_loop(start, P.np()) {
+        let p = P.getParticle(i);
+        p.setForces(get_null_vector());
+    }
+    /*
+    initialize(lut.data, lut.rows()*lut.cols()*sizeof[u8]() as size_t, 0i32); 
+    
+    for i in par_loop(start, P.np()){
+        for j in loop(start, P.np()) {
+            if(i != j)
+            {
+                let mut row = i;
+                let mut col = j; 
+                if (i > j) {
+                    row = j; 
+                    col = i;
+                } 
+                
+                if(row >= lut.rows()) {
+                    let tmp = col;
+                    col = row;
+                    row = tmp;
+                    col -= lut.rows();
+                    row -= lut.rows();
+
+                }
+                let flag = cmpxchg(&mut bitcast[&mut[u8]](lut.data)(row*lut.cols() + col), 0u8, 1u8);
+                if (flag(1) == true)
+                {
+                    force(P.getParticle(i), P.getParticle(j), P.constants());
+                }
+            }
+        }
+    }*/
+    
+    for i in loop(start, P.np()) {
+        for j in loop(i+(1 as size_t), P.np()) {
+            let p1 = P.getParticle(i);
+            let p2 = P.getParticle(j);
+            force(p1, p2, P.constants());
+        }
+    }
+}
+
+fn update(P: ParticleSystem, dt: real, f: fn(Particle, real) -> ()) -> () {
+    for i in loop(0 as size_t, P.np()) {
+        f(P.getParticle(i), dt);
+    }
+}
+
+fn move_particles(P: ParticleSystem) -> () {
+    for i in loop(0 as size_t, P.np()) {
+        boundary(P.getParticle(i), P.l());
+    }
+}
+
+
+fn fprint_particle_system(fname: &[u8], step: size_t, P: ParticleSystem) -> ()
+{
+    let start = 0 as size_t;
+    let fp = open_file(fname);
+    let N = P.np();
+
+    fprint_line(fp, "# vtk DataFile Version 2.0");
+    fprint_string(fp, "Step ");
+    fprint_size_t(fp, step);
+    fprint_line(fp, " data");
+    fprint_line(fp, "ASCII");
+    fprint_line(fp, "DATASET UNSTRUCTURED_GRID");
+    fprint_string(fp, "POINTS ");
+    fprint_size_t(fp, N);
+    fprint_line(fp, " double");
+    loop(start, N, |i| {
+        let p = P.getParticle(i);
+        p.getCoordinates().fprint(fp);
+        fprint_string(fp, "\n");
+    });
+    fprint_string(fp, "\n\n");
+    fprint_string(fp, "CELLS ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, " ");
+    fprint_size_t(fp, (2 as size_t)*N);
+
+    fprint_string(fp, "\n");
+    loop(start, N, |i| {
+        fprint_string(fp, "1 ");
+        fprint_size_t(fp, i);
+        fprint_string(fp, "\n");
+    });
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "CELL_TYPES ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    loop(start, N, |i| {
+        fprint_string(fp, "1");
+        fprint_string(fp, "\n");
+    });
+
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "POINT_DATA ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    @loop(start, DIM, |d| {
+        fprint_string(fp, "SCALARS velocity_dim_");
+        fprint_size_t(fp, d);
+        fprint_line(fp, " double");
+        fprint_line(fp, "LOOKUP_TABLE default");
+        $loop(start, N, |i| {
+
+            let p = P.getParticle(i);
+            let v = p.getVelocities();
+            let x = v.x();
+            fprint_double(fp, x(d));
+            fprint_string(fp, "\n");
+        });
+        fprint_string(fp, "\n");
+    });
+
+    @loop(start, DIM, |d| {
+        fprint_string(fp, "SCALARS force_dim_");
+        fprint_size_t(fp, d);
+        fprint_line(fp, " double");
+        fprint_line(fp, "LOOKUP_TABLE default");
+        $loop(start, N, |i| {
+
+            let p = P.getParticle(i);
+            let v = p.getForces();
+            let x = v.x();
+            fprint_double(fp, x(d));
+            fprint_string(fp, "\n");
+        });
+        fprint_string(fp, "\n");
+    });
+
+
+    fprint_string(fp, "SCALARS mass");
+    fprint_line(fp, " double");
+    fprint_line(fp, "LOOKUP_TABLE default");
+    loop(start, N, |i| {
+        let p = P.getParticle(i);
+        let m = p.getMass();
+        fprint_double(fp, m);
+        fprint_string(fp, "\n");
+    });
+    fprint_string(fp, "\n");
+
+
+    close_file(fp);
+}
+
+
+fn print_statistics(P: ParticleSystem, t: real) -> ()
+{
+
+    let sqr = |x : real| {x*x};
+    let mut e : real = 0.0;
+    loop(0 as size_t, P.np(), |i| {
+        let p = P.getParticle(i);
+        let v = p.getVelocities().reduce(sqr, |x,y|{x+y}, 0.0);
+        e += 0.5 * p.getMass() * v;
+    });
+    print_string("t: ");
+    print_double(t);
+    print_string("\tE: ");
+    print_double(e);
+    print_char('\n');
+}

apps/impala/algorithm/linked_cell/cpu.impala

+
+fn compute_force(P: ParticleSystem, force: fn(Particle, Particle, Constants) -> ()) -> ()
+{
+    let z = 0 as size_t;
+    let ic_start = [z,z,z];
+    let ic_end = P.nc();
+    for ic, pl1 in @iterate_over_particle_system(ic_start, ic_end, P) {
+        let p1 = get_particle_from_node(*pl1);
+        p1.setForces(get_null_vector());
+    }
+    for ic, pl1 in @iterate_over_particle_system(ic_start, ic_end, P) {
+        let p1 = get_particle_from_node(*pl1);
+        let jc_start = [if(ic(0) > z) {ic(0) - 1 as size_t} else {ic(0)}, 
+                        if(ic(1) > z) {ic(1) - 1 as size_t} else{ic(1)}, 
+                        if(ic(2) > z) {ic(2) - 1 as size_t} else {ic(2)}];
+        let jc_end =   [if(ic(0) + 1 as size_t < ic_end(0)) {ic(0) + 2 as size_t} else {ic_end(0)}, 
+                        if(ic(1) + 1 as size_t < ic_end(1)) {ic(1) + 2 as size_t} else {ic_end(1)}, 
+                        if(ic(2) + 1 as size_t < ic_end(2)) {ic(2) + 2 as size_t} else {ic_end(2)}];
+        for jc, pl2 in iterate_over_particle_system(jc_start, jc_end, P) {
+            if((pl1 as ptr_t < pl2 as ptr_t)) {
+                let p2 = get_particle_from_node(*pl2);
+                force(p1, p2, P.constants());   
+            }
+        }
+    }
+}
+
+fn update(P: ParticleSystem, dt: real, f: fn(Particle, real) -> ()) -> ()
+{
+    let z = 0 as size_t;
+    let ic_start = [z,z,z];
+    let ic_end = P.nc();
+    for ic, pl in @iterate_over_particle_system(ic_start, ic_end, P) {
+        let p = get_particle_from_node(*pl);
+        f(p, dt);
+    }
+}
+
+fn move_particles(P: ParticleSystem) -> () {
+    let z = 0 as size_t;
+    let ic_start = [z,z,z];
+    let ic_end = P.nc();
+    for ic, root in @iterate_over_particle_cells(ic_start, ic_end, P) {
+        let mut q : &mut &ParticleList = root;
+        let mut i : &ParticleList = *q;
+        
+        while(i != 0 as &ParticleList) {    
+            let p = get_particle_from_node(*i); 
+            boundary(p, P.l());
+            let kc = compute_cell_position(p, P.nc(), P.l());
+            let mut moving_node : &mut ParticleList;
+            if(ic(0) != kc(0) || ic(1) != kc(1) || ic(2) != kc(2)) {
+                moving_node = remove(q);
+                insert(P.head(index(kc, P.nc())), moving_node);
+            }
+            else {
+                q = (&(*i).next) as &mut& ParticleList;
+            }
+            i = (*q);
+        }
+    }
+
+}
+fn fprint_particle_system(fname: &[u8], step: size_t, P: ParticleSystem) -> ()
+{
+    let start = 0 as size_t;
+    let fp = open_file(fname);
+    let N = P.np();
+
+    fprint_line(fp, "# vtk DataFile Version 2.0");
+    fprint_string(fp, "Step ");
+    fprint_size_t(fp, step);
+    fprint_line(fp, " data");
+    fprint_line(fp, "ASCII");
+    fprint_line(fp, "DATASET UNSTRUCTURED_GRID");
+    fprint_string(fp, "POINTS ");
+    fprint_size_t(fp, N);
+    fprint_line(fp, " double");
+    
+    let addresses = P.addresses();
+
+    for i in loop(start, N) {
+        let p = get_particle(addresses.get(i));
+        p.getCoordinates().fprint(fp);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n\n");
+    fprint_string(fp, "CELLS ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, " ");
+    fprint_size_t(fp, (2 as size_t)*N);
+
+    fprint_string(fp, "\n");
+    for i in loop(start, N){
+        fprint_string(fp, "1 ");
+        fprint_size_t(fp, i);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "CELL_TYPES ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    for i in loop(start, N) {
+        fprint_string(fp, "1");
+        fprint_string(fp, "\n");
+    }
+
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "POINT_DATA ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    for d in loop(start, DIM) {
+        fprint_string(fp, "SCALARS velocity_dim_");
+        fprint_size_t(fp, d);
+        fprint_line(fp, " double");
+        fprint_line(fp, "LOOKUP_TABLE default");
+        for i in loop(start, N) {
+            let p = get_particle(addresses.get(i));
+            let v = p.getVelocities();
+            let x = v.x();
+            fprint_double(fp, x(d));
+            fprint_string(fp, "\n");
+        }
+        fprint_string(fp, "\n");
+    }
+    for d in loop(start, DIM) {
+        fprint_string(fp, "SCALARS force_dim_");
+        fprint_size_t(fp, d);
+        fprint_line(fp, " double");
+        fprint_line(fp, "LOOKUP_TABLE default");
+        for i in loop(start, N) {
+            let p = get_particle(addresses.get(i));
+            let v = p.getForces();
+            let x = v.x();
+            fprint_double(fp, x(d));
+            fprint_string(fp, "\n");
+        }
+        fprint_string(fp, "\n");
+    }
+
+    fprint_string(fp, "SCALARS mass");
+    fprint_line(fp, " double");
+    fprint_line(fp, "LOOKUP_TABLE default");
+    for i in loop(start, N) {
+        let p = get_particle(addresses.get(i));
+        let m = p.getMass();
+        fprint_double(fp, m);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n");
+    close_file(fp);
+}
+
+/*
+fn fprint_particle_system(fname: &[u8], step: size_t, P: ParticleSystem) -> ()
+{
+    let start = 0 as size_t;
+    let fp = open_file(fname);
+    let N = P.np();
+
+    let z = 0 as size_t;
+    let ic_start = [z,z,z];
+    let ic_end = P.nc();
+
+    fprint_line(fp, "# vtk DataFile Version 2.0");
+    fprint_string(fp, "Step ");
+    fprint_size_t(fp, step);
+    fprint_line(fp, " data");
+    fprint_line(fp, "ASCII");
+    fprint_line(fp, "DATASET UNSTRUCTURED_GRID");
+    fprint_string(fp, "POINTS ");
+    fprint_size_t(fp, N);
+    fprint_line(fp, " double");
+
+    for ic, pl in iterate_over_particle_system(ic_start, ic_end, P) {
+        let p = get_particle_from_cell(*pl);
+        p.getCoordinates().fprint(fp);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n\n");
+    fprint_string(fp, "CELLS ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, " ");
+    fprint_size_t(fp, (2 as size_t)*N);
+
+    fprint_string(fp, "\n");
+    for i in loop(start, N){
+        fprint_string(fp, "1 ");
+        fprint_size_t(fp, i);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "CELL_TYPES ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    for i in loop(start, N) {
+        fprint_string(fp, "1");
+        fprint_string(fp, "\n");
+    }
+
+    fprint_string(fp, "\n\n");
+
+
+    fprint_string(fp, "POINT_DATA ");
+    fprint_size_t(fp, N);
+    fprint_string(fp, "\n");
+    for d in loop(start, DIM) {
+        fprint_string(fp, "SCALARS velocity_dim_");
+        fprint_size_t(fp, d);
+        fprint_line(fp, " double");
+        fprint_line(fp, "LOOKUP_TABLE default");
+        for ic, pl in iterate_over_particle_system(ic_start, ic_end, P) {
+            let p = get_particle_from_cell(*pl);
+            let v = p.getVelocities();
+            let x = v.x();
+            fprint_double(fp, x(d));
+            fprint_string(fp, "\n");
+        }
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "SCALARS mass");
+    fprint_line(fp, " double");
+    fprint_line(fp, "LOOKUP_TABLE default");
+    for ic, pl in iterate_over_particle_system(ic_start, ic_end, P) {
+        let p = get_particle_from_cell(*pl);
+        let m = p.getMass();
+        fprint_double(fp, m);
+        fprint_string(fp, "\n");
+    }
+    fprint_string(fp, "\n");
+
+
+    close_file(fp);
+    
+}
+*/
+
+fn print_statistics(P: ParticleSystem, t: real) -> ()
+{
+
+    let sqr = |x : real| {x*x};
+    let e : real = 0.0;
+    /*
+    loop(0 as size_t, P.np(), |i| {
+        let v = P.velocities().get(i).reduce(sqr, |x,y|{x+y}, 0.0);
+        e += 0.5 * P.masses().get(i) * v;
+    });*/
+    print_string("t: ");
+    print_double(t);
+    print_string("\tE: ");
+    print_double(e);
+    print_char('\n');
+}

apps/impala/boundary/rigid_walls.impala

+fn boundary(p: Particle, l: Vector) -> () {
+    let mut x = p.getCoordinates().x();
+    let mut change = false;
+    for d in @loop(0 as size_t, DIM) {
+        if(x(d) < 0.0) {
+            x(d) = 0.0+EPS;
+            change = true;
+        }
+        else if(x(d) > l.x()(d)) {
+            x(d) = l.x()(d)-EPS;
+            change = true;
+        }
+    }
+    if(change == true) {
+        p.setCoordinates(get_vector(x));
+    }
+}   

apps/impala/clean_build

+#!/bin/sh
+rm -rf Makefile CMakeFiles cmake_install.cmake CMakeCache.txt CTestTestfile.cmake *.o *.ll md

apps/impala/clib/allocate.c

+#include <stdlib.h>
+#include <string.h>
+
+static size_t ALIGNMENT = 64ul;
+unsigned char * allocate(size_t size) {
+    return (unsigned char *)aligned_alloc(ALIGNMENT, size);
+}
+
+
+unsigned char * allocate_and_initialize(size_t size, int val) {
+    void *ptr = aligned_alloc(ALIGNMENT, size);
+    memset(ptr, val, size);
+    return (unsigned char *)ptr;
+}
+
+unsigned char * reallocate(unsigned char *ptr, size_t size) {
+    return (unsigned char *)realloc((void *)ptr, size);
+}
+
+void initialize(unsigned char * ptr, size_t size, int val) {
+    memset((void *)ptr, val, size);
+}
+
+void deallocate(unsigned char *ptr) {
+    free(ptr);
+}

apps/impala/clib/allocate.h

+unsigned char * allocate(size_t size);
+unsigned char * allocate_and_initialize(size_t size, int val);
+unsigned char * reallocate(unsigned char *ptr, size_t size);
+void initialize(unsigned char *ptr, size_t size, int val);
+void deallocate(unsigned char *ptr);

apps/impala/clib/fileIO.c

+#include <stdio.h>
+#include <stdint.h>
+
+
+
+void size_t_to_string(size_t i, unsigned char *buf, size_t bufsize) {
+    snprintf(buf, bufsize, "%lu", i);
+}
+
+void generate_filename(size_t i, unsigned char *buf, size_t bufsize) {
+    snprintf(buf, bufsize, "../impala_vtk/particles%lu.vtk", i); 
+}
+
+uintptr_t open_file(unsigned char const *fname) {
+    FILE *fp = fopen(fname, "w");
+    if(fp == NULL) {
+        fprintf(stderr, "open_file: Could not open file %s", fname);
+        perror("");
+        return 0;
+    }
+    else {
+        return (uintptr_t)fp;
+    }
+}
+void close_file(uintptr_t fp) {
+    fclose((FILE *)fp);
+}
+
+void fprint_string(uintptr_t fp, unsigned char const *str) {
+    fputs(str, (FILE *)fp);
+}
+
+void fprint_line(uintptr_t fp, unsigned char const *str) {
+    fprintf((FILE *)fp, "%s\n", str);
+}
+void fprint_double(uintptr_t fp, double d) {
+    fprintf((FILE *)fp, "%f", d);
+}
+
+void fprint_float(uintptr_t fp, float f) {
+    fprintf((FILE *)fp, "%f", f);
+}
+
+void fprint_int(uintptr_t fp, int i) {
+    fprintf((FILE *)fp, "%d", i);
+}
+
+
+void fprint_size_t(uintptr_t fp, size_t i) {
+    fprintf((FILE *)fp, "%lu", i);
+}
+
+void fprint_char(uintptr_t fp, unsigned char c) {
+    fputc((int)c, (FILE *)fp);
+}

apps/impala/clib/fileIO.h

+#include <stdio.h>
+#include <stdint.h>
+
+
+
+void size_t_to_string(size_t i, unsigned char *buf, size_t bufsize); 
+
+void generate_filename(size_t i, unsigned char *buf, size_t bufsize);
+
+uintptr_t open_file(unsigned char const *fname);
+
+void close_file(uintptr_t fp);
+
+void fprint_string(uintptr_t fp, unsigned char const *str);
+
+void fprint_line(uintptr_t fp, unsigned char const *str);
+void fprint_double(uintptr_t fp, double d);
+
+void fprint_float(uintptr_t fp, float f);
+
+void fprint_int(uintptr_t fp, int i);
+
+void fprint_size_t(uintptr_t fp, size_t i);
+
+void fprint_char(uintptr_t fp, unsigned char c);
+

apps/impala/common.impala

+extern "C"
+{
+    fn size_t_to_string(size_t, &mut[u8], size_t) -> ();
+    fn generate_filename(size_t , &mut[u8], size_t) -> ();
+    fn open_file(&[u8]) -> u64;
+    fn close_file(u64) -> ();
+    fn fprint_string(u64, &[u8]) -> ();
+    fn fprint_line(u64, &[u8]) -> ();
+    fn fprint_double(u64, f64) -> ();
+    fn fprint_float(u64, f32) -> ();
+    fn fprint_int(u64, i32) -> ();
+    fn fprint_size_t(u64, size_t) -> ();
+    fn fprint_char(u64, u8) -> ();
+    fn allocate(size_t) -> &u8;
+    fn allocate_and_initialize(size_t, i32) -> &u8;
+    fn initialize(&u8, size_t, i32) -> ();
+    fn deallocate(&u8) -> ();
+}
+
+type size_t = u64;
+type ptr_t = u64;
+type real = f64;
+static INF = 1e308f64;
+static EPS = 1e-12;
+static DIM = 3 as size_t;
+static SIZEOFPARTICLE = ((1 as size_t) + (4 as size_t) * DIM) * (sizeof[real]() as size_t); 
+
+
+struct Vector
+{
+    x: fn() -> [real * 3],