Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
dedbb983
Commit
dedbb983
authored
Sep 08, 2017
by
Matthieu Schaller
Browse files
Added particle type solely designated to debugging the SELF/PAIR.
parent
d9ebdcfe
Changes
12
Hide whitespace changes
Inline
Side-by-side
configure.ac
View file @
dedbb983
...
...
@@ -627,7 +627,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax default: gadget2@:>@]
[Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax
debug
default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
...
...
@@ -639,6 +639,9 @@ case "$with_hydro" in
minimal)
AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
;;
debug)
AC_DEFINE([DEBUG_INTERACTIONS_SPH], [1], [Debug SPH])
;;
hopkins)
AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
;;
...
...
examples/SedovBlast_3D/sedov.yml
View file @
dedbb983
...
...
@@ -32,4 +32,5 @@ SPH:
# Parameters related to the initial conditions
InitialConditions
:
file_name
:
./sedov.hdf5
# The file to read
h_scaling
:
3.3
src/debug.c
View file @
dedbb983
...
...
@@ -21,16 +21,15 @@
******************************************************************************/
/* Config parameters. */
#include
"../config.h"
/* This object's header. */
#include
"debug.h"
/* Some standard headers. */
#include
<float.h>
#include
<stdio.h>
#include
<unistd.h>
/* This object's header. */
#include
"debug.h"
/* Local includes. */
#include
"active.h"
#include
"cell.h"
...
...
@@ -41,7 +40,9 @@
#include
"space.h"
/* Import the right hydro definition */
#if defined(MINIMAL_SPH)
#if defined(DEBUG_INTERACTIONS_SPH)
#include
"./hydro/DebugInteractions/hydro_debug.h"
#elif defined(MINIMAL_SPH)
#include
"./hydro/Minimal/hydro_debug.h"
#elif defined(GADGET2_SPH)
#include
"./hydro/Gadget2/hydro_debug.h"
...
...
src/debug.h
View file @
dedbb983
...
...
@@ -19,6 +19,9 @@
#ifndef SWIFT_DEBUG_H
#define SWIFT_DEBUG_H
/* Config parameters. */
#include
"../config.h"
/* Includes. */
#include
"cell.h"
#include
"part.h"
...
...
src/hydro.h
View file @
dedbb983
...
...
@@ -27,7 +27,11 @@
#include
"kernel_hydro.h"
/* Import the right functions */
#if defined(MINIMAL_SPH)
#if defined(DEBUG_INTERACTIONS_SPH)
#include
"./hydro/DebugInteractions/hydro.h"
#include
"./hydro/DebugInteractions/hydro_iact.h"
#define SPH_IMPLEMENTATION "Debug SELF/PAIR"
#elif defined(MINIMAL_SPH)
#include
"./hydro/Minimal/hydro.h"
#include
"./hydro/Minimal/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH (e.g. Price 2010)"
...
...
src/hydro/DebugInteractions/hydro.h
0 → 100644
View file @
dedbb983
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_H
#define SWIFT_DEBUG_INTERACTIONS_HYDRO_H
/**
* @file DebugInteractions/hydro.h
* @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
*/
#include
"hydro_properties.h"
#include
"hydro_space.h"
#include
"part.h"
#include
<float.h>
/**
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_internal_energy
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_pressure
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_entropy
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_soundspeed
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the density of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_density
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the mass of a particle
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_mass
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the time derivative of internal energy of a particle
*
* We assume a constant density.
*
* @param p The particle of interest
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_internal_energy_dt
(
const
struct
part
*
restrict
p
)
{
return
0
.
f
;
}
/**
* @brief Returns the time derivative of internal energy of a particle
*
* We assume a constant density.
*
* @param p The particle of interest.
* @param du_dt The new time derivative of the internal energy.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_set_internal_energy_dt
(
struct
part
*
restrict
p
,
float
du_dt
)
{}
/**
* @brief Computes the hydro time-step of a given particle
*
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
*
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_compute_timestep
(
const
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
const
struct
hydro_props
*
restrict
hydro_properties
)
{
return
FLT_MAX
;
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @param p The particle to act upon.
* @param dt Physical time step of the particle during the next step.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_timestep_extra
(
struct
part
*
p
,
float
dt
)
{}
/**
* @brief Prepares a particle for the density calculation.
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the variaous density tasks
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_init_part
(
struct
part
*
restrict
p
,
const
struct
hydro_space
*
hs
)
{
for
(
int
i
=
0
;
i
<
256
;
++
i
)
p
->
ids_ngbs_density
[
i
]
=
-
1
;
p
->
num_ngb_density
=
0
;
p
->
density
.
wcount
=
0
.
f
;
p
->
density
.
wcount_dh
=
0
.
f
;
}
/**
* @brief Finishes the density calculation.
*
* Multiplies the density and number of neighbours by the appropiate constants
* and add the self-contribution term.
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_density
(
struct
part
*
restrict
p
)
{
/* Some smoothing length multiples. */
const
float
h
=
p
->
h
;
const
float
h_inv
=
1
.
0
f
/
h
;
/* 1/h */
const
float
h_inv_dim
=
pow_dimension
(
h_inv
);
/* 1/h^d */
const
float
h_inv_dim_plus_one
=
h_inv_dim
*
h_inv
;
/* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p
->
density
.
wcount
+=
kernel_root
;
p
->
density
.
wcount_dh
-=
hydro_dimension
*
kernel_root
;
/* Finish the calculation by inserting the missing h-factors */
p
->
density
.
wcount
*=
h_inv_dim
;
p
->
density
.
wcount_dh
*=
h_inv_dim_plus_one
;
}
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_part_has_no_neighbours
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{
/* Some smoothing length multiples. */
const
float
h
=
p
->
h
;
const
float
h_inv
=
1
.
0
f
/
h
;
/* 1/h */
const
float
h_inv_dim
=
pow_dimension
(
h_inv
);
/* 1/h^d */
/* Re-set problematic values */
p
->
density
.
wcount
=
kernel_root
*
kernel_norm
*
h_inv_dim
;
p
->
density
.
wcount_dh
=
0
.
f
;
}
/**
* @brief Prepare a particle for the force calculation.
*
* Computes viscosity term, conduction term and smoothing length gradient terms.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param ti_current The current time (on the timeline)
* @param timeBase The minimal time-step size
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_prepare_force
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{}
/**
* @brief Reset acceleration fields of a particle
*
* Resets all hydro acceleration and time derivative fields in preparation
* for the sums taking place in the variaous force tasks
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_reset_acceleration
(
struct
part
*
restrict
p
)
{
for
(
int
i
=
0
;
i
<
256
;
++
i
)
p
->
ids_ngbs_force
[
i
]
=
-
1
;
p
->
num_ngb_force
=
0
;
p
->
force
.
h_dt
=
0
.
f
;
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param p The particle.
* @param xp The extended data of this particle.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_reset_predicted_values
(
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
)
{}
/**
* @brief Predict additional particle fields forward in time when drifting
*
* @param p The particle
* @param xp The extended data of the particle
* @param dt The drift time-step.
* @param t0 The time at the start of the drift
* @param t1 The time at the end of the drift
* @param timeBase The minimal time-step size
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_predict_extra
(
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
float
dt
)
{}
/**
* @brief Finishes the force calculation.
*
* Multiplies the forces and accelerationsby the appropiate constants
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_force
(
struct
part
*
restrict
p
)
{}
/**
* @brief Kick the additional variables
*
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_kick_extra
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
float
dt
)
{}
/**
* @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_convert_quantities
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_first_init_part
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{
p
->
time_bin
=
0
;
xp
->
v_full
[
0
]
=
p
->
v
[
0
];
xp
->
v_full
[
1
]
=
p
->
v
[
1
];
xp
->
v_full
[
2
]
=
p
->
v
[
2
];
hydro_reset_acceleration
(
p
);
hydro_init_part
(
p
,
NULL
);
}
#endif
/* SWIFT_DEBUG_INTERACTIONS_HYDRO_H */
src/hydro/DebugInteractions/hydro_debug.h
0 → 100644
View file @
dedbb983
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
#define SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
/**
* @file DebugInteractions/hydro_debug.h
* @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_debug_particle
(
const
struct
part
*
p
,
const
struct
xpart
*
xp
)
{
printf
(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e]
\n
a=[%.3e,%.3e,%.3e],
\n
"
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, dh/dt=%.3e time_bin=%d
\n
"
,
p
->
x
[
0
],
p
->
x
[
1
],
p
->
x
[
2
],
p
->
v
[
0
],
p
->
v
[
1
],
p
->
v
[
2
],
xp
->
v_full
[
0
],
xp
->
v_full
[
1
],
xp
->
v_full
[
2
],
p
->
a_hydro
[
0
],
p
->
a_hydro
[
1
],
p
->
a_hydro
[
2
],
p
->
h
,
p
->
density
.
wcount
,
p
->
density
.
wcount_dh
,
p
->
force
.
h_dt
,
p
->
time_bin
);
}
#endif
/* SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H */
src/hydro/DebugInteractions/hydro_iact.h
0 → 100644
View file @
dedbb983
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H
#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H
/**
* @file DebugInteractions/hydro_iact.h
* @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
*/
/**
* @brief Density loop
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
wi
,
wi_dx
;
float
wj
,
wj_dx
;
/* Get r and r inverse. */
const
float
r
=
sqrtf
(
r2
);
/* Compute the kernel function for pi */
const
float
hi_inv
=
1
.
f
/
hi
;
const
float
ui
=
r
*
hi_inv
;
kernel_deval
(
ui
,
&
wi
,
&
wi_dx
);
/* Compute contribution to the number of neighbours */
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
/* Compute the kernel function for pj */
const
float
hj_inv
=
1
.
f
/
hj
;
const
float
uj
=
r
*
hj_inv
;
kernel_deval
(
uj
,
&
wj
,
&
wj_dx
);
/* Compute contribution to the number of neighbours */
pj
->
density
.
wcount
+=
wj
;
pj
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wj
+
uj
*
wj_dx
);
/* Update ngb counters */
pi
->
ids_ngbs_density
[
pi
->
num_ngb_density
]
=
pj
->
id
;
++
pi
->
num_ngb_density
;
pj
->
ids_ngbs_density
[
pj
->
num_ngb_density
]
=
pi
->
id
;
++
pj
->
num_ngb_density
;
}
/**
* @brief Density loop (non-symmetric version)
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
wi
,
wi_dx
;
/* Get r and r inverse. */
const
float
r
=
sqrtf
(
r2
);
/* Compute the kernel function */
const
float
hi_inv
=
1
.
0
f
/
hi
;
const
float
ui
=
r
*
hi_inv
;
kernel_deval
(
ui
,
&
wi
,
&
wi_dx
);
/* Compute contribution to the number of neighbours */
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
/* Update ngb counters */
pi
->
ids_ngbs_density
[
pi
->
num_ngb_density
]
=
pj
->
id
;
++
pi
->
num_ngb_density
;
}
/**
* @brief Force loop
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_force
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
/* Update ngb counters */
pi
->
ids_ngbs_force
[
pi
->
num_ngb_force
]
=
pj
->
id
;
++
pi
->
num_ngb_force
;
pj
->
ids_ngbs_force
[
pj
->
num_ngb_force
]
=
pi
->
id
;
++
pj
->
num_ngb_force
;
}
/**
* @brief Force loop (non-symmetric version)
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_force
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
/* Update ngb counters */
pi
->
ids_ngbs_force
[
pi
->
num_ngb_force
]
=
pj
->
id
;
++
pi
->
num_ngb_force
;
}
#endif
/* SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H */
src/hydro/DebugInteractions/hydro_io.h
0 → 100644
View file @
dedbb983
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H
#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H
/**
* @file DebugInteractions/hydro_io.h
* @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
*/
#include
"adiabatic_index.h"
#include
"hydro.h"
#include
"io_properties.h"
#include
"kernel_hydro.h"
/**
* @brief Specifies which particle fields to read from a dataset
*
* @param parts The particle array.
* @param list The list of i/o properties to read.
* @param num_fields The number of i/o fields to read.
*/
void
hydro_read_particles
(
struct
part
*
parts
,
struct
io_props
*
list
,
int
*
num_fields
)
{
*
num_fields
=
4
;
/* List what we want to read */
list
[
0
]
=
io_make_input_field
(
"Coordinates"
,
DOUBLE
,
3
,
COMPULSORY
,
UNIT_CONV_LENGTH
,
parts
,
x
);
list
[
1
]
=
io_make_input_field
(
"Velocities"
,
FLOAT
,
3
,
COMPULSORY
,
UNIT_CONV_SPEED
,
parts
,
v
);
list
[
2
]
=
io_make_input_field
(
"SmoothingLength"
,
FLOAT
,
1
,
COMPULSORY
,
UNIT_CONV_LENGTH
,
parts
,
h
);
list
[
3
]
=
io_make_input_field
(
"ParticleIDs"
,
ULONGLONG
,
1
,
COMPULSORY
,
UNIT_CONV_NO_UNITS
,
parts
,
id
);
}
float
convert_u
(
struct
engine
*
e
,
struct
part
*
p
)
{
return
hydro_get_internal_energy
(
p
);
}
float
convert_P
(
struct
engine
*
e
,
struct
part
*
p
)
{
return
hydro_get_pressure
(
p
);
}
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.