diff --git a/CMakeLists.txt b/CMakeLists.txt index e2bdfa7c..592db256 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,7 @@ OPTION(USE_MPI "use mpi, if found" ON) # use -DUSE_MPI=off to force disable use OPTION(BUILD_SHARED_LIBS "Builds shared libraries if ON, static libraries if OFF" ON) OPTION(SET_RPATH "set the rpath of the executable to the found libraries" ON) OPTION(WITH_NSOLVER "compile support for the nsolver module" ON) +OPTION(WITH_ILC "compile programs for inclined layer convection" OFF) OPTION(WITH_PYTHON "build a python wrapper, using boost-python" OFF) OPTION(WITH_HDF5CXX "enable legacy file format using hdf5 cxx" OFF) OPTION(WITH_GTEST "enable gtest unit testing" ON) @@ -190,6 +191,15 @@ add_subdirectory(programs) add_subdirectory(tools) add_subdirectory(examples) +if (WITH_ILC) + if (NOT WITH_NSOLVER) + message(FATAL_ERROR "Compiling the ILC module requires -DWITH_NSOLVER=ON") + endif() + add_subdirectory(modules/ilc) + add_subdirectory(modules/ilc/programs) + add_subdirectory(modules/ilc/examples) +endif () + # Check if we want to build the python wrapper and have boost-python # If Python and boost.python are there, build python interface if (WITH_PYTHON) @@ -204,6 +214,9 @@ endif () enable_testing() add_subdirectory(tests) +if (WITH_ILC) + add_subdirectory(modules/ilc/tests) +endif () if (WITH_PYTHON) add_subdirectory(python-wrapper/tests) endif() @@ -292,6 +305,12 @@ else () message(" nsolver: disabled") endif () +if (WITH_ILC) + message(" Inclined layer convection: enabled") +else () + message(" Inclined layer convection: disabled") +endif () + if (WITH_PYTHON) message(" Python wrapper: enabled") else () diff --git a/modules/ilc/CMakeLists.txt b/modules/ilc/CMakeLists.txt new file mode 100644 index 00000000..7c38a993 --- /dev/null +++ b/modules/ilc/CMakeLists.txt @@ -0,0 +1,32 @@ +# +# This file is a part of channelflow version 2.0 https://channelflow.ch. +# License is GNU GPL version 2 or later: ./LICENCE +# +# Include and configure the Inclined Layer Convection module of channelflow. +# +set( + ilc_SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/ilcflags.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/obe.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ilc.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ilcdsi.cpp +) + +set( + ilc_HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/ilcflags.h + ${CMAKE_CURRENT_SOURCE_DIR}/obe.h + ${CMAKE_CURRENT_SOURCE_DIR}/ilc.h + ${CMAKE_CURRENT_SOURCE_DIR}/ilcdsi.h +) + +# Define the target with appropriate dependencies +install_channelflow_library(ilc) +target_link_fftw(ilc PRIVATE) +target_link_libraries(ilc PUBLIC chflow) +if (WITH_NSOLVER) + target_link_libraries(ilc PUBLIC nsolver) +endif () + +# Install header files +install(FILES ${ilc_HEADERS} DESTINATION include/modules/ilc) diff --git a/modules/ilc/examples/CMakeLists.txt b/modules/ilc/examples/CMakeLists.txt new file mode 100644 index 00000000..6a5aedb9 --- /dev/null +++ b/modules/ilc/examples/CMakeLists.txt @@ -0,0 +1,6 @@ +set(ilc_EXAMPLES benchmark_ilc) + +foreach (program ${ilc_EXAMPLES}) + install_channelflow_application(${program} examples) + target_link_libraries(${program}_app PUBLIC ilc) +endforeach (program) diff --git a/modules/ilc/examples/benchmark_ilc.cpp b/modules/ilc/examples/benchmark_ilc.cpp new file mode 100644 index 00000000..a156458e --- /dev/null +++ b/modules/ilc/examples/benchmark_ilc.cpp @@ -0,0 +1,166 @@ +/** + * The classic benchmark program, like for pure shear flow + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include "channelflow/flowfield.h" +#include "modules/ilc/ilc.h" + +#ifdef HAVE_MPI +#include +#endif +#include "channelflow/cfmpi.h" + +using namespace std; +using namespace chflow; + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + int taskid = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &taskid); +#endif + if (taskid == 0) + cout << "Starting channelflow benchmark" << endl; + + string purpose( + "This program loads the field uinit from the harddisk and integrates\nit for 10 time units. The resulting " + "field is compared to the field\nufinal. Wall-clock time elapsed for each timeunit as well as an\naverage " + "are calculated."); + ArgList args(argc, argv, purpose); + int nproc0 = args.getint("-np0", "-nproc0", 0, "number of processes for transpose/parallel ffts"); + int nproc1 = args.getint("-np1", "-nproc1", 0, "number of processes for slice fft"); + bool fftwmeasure = args.getflag("-fftwmeasure", "--fftwmeasure", "use fftw_measure instead of fftw_patient"); + bool fftwwisdom = args.getflag("-fftwwisdom", "--fftwwisdom", "try loading fftw wisdom"); + bool saveresult = args.getflag("-s", "--saveresult", "save resulting field to uresult.nc"); + string dir = args.getstr("-d", "--directory", ".", "directory where fields are found and stored"); + args.check(); + + if (taskid == 0) + cout << "Creating CfMPI object..." << flush; + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + if (taskid == 0) + cout << "done" << endl; + + if (taskid == 0) + cout << "Loading FlowFields..." << endl; + FlowField u(dir + "/uinit", cfmpi); + FlowField temp(dir + "/tinit", cfmpi); + if (taskid == 0) + cout << "done" << endl; + + if (u.taskid() == 0) + cout << "Calculating l2norm..." << flush; + Real ul2n = L2Norm(u); + Real tl2n = L2Norm(temp); + if (taskid == 0) + cout << "done" << endl; + + if (u.taskid() == 0) { + cout << "================================================================\n"; + cout << purpose << endl << endl; + cout << "Distribution of processes is " << u.nproc0() << "x" << u.nproc1() << endl; + } + if (fftwwisdom && u.taskid() == 0) { + cout << "Loading fftw wisdom" << endl; + fftw_loadwisdom(); + } + if (u.taskid() == 0) + cout << "Optimizing FFTW..." << flush; + { + FlowField utmp(u); + if (fftwmeasure) + utmp.optimizeFFTW(FFTW_MEASURE); + else + utmp.optimizeFFTW(FFTW_PATIENT); + fftw_savewisdom(); + } + if (u.taskid() == 0) + cout << "done" << endl; + ul2n = L2Norm(u); + tl2n = L2Norm(temp); + if (u.taskid() == 0) + cout << " L2Norm(u) == " << ul2n << endl; + if (u.taskid() == 0) + cout << " L2Norm(T) == " << tl2n << endl; + + // Define integration parameters + const int n = 40; // take n steps between printouts + const Real dt = 1.0 / n; // integration timestep + + // Define DNS parameters + ILCFlags flags; + + // Run at default flags. If you change them, recreate the test files. + flags.dt = dt; + flags.verbosity = Silent; + + if (u.taskid() == 0) + cout << "Building FlowField q..." << flush; + vector fields = { + u, temp, FlowField(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b(), u.cfmpi(), Spectral, Spectral)}; + if (u.taskid() == 0) + cout << "done" << endl; + if (u.taskid() == 0) + cout << "Building dns..." << flush; + ILC ilc(fields, flags); + if (u.taskid() == 0) + cout << "done" << endl; + + Real avtime = 0; + int i = 0; + Real T = 10; + for (Real t = 0; t <= T; t += 1) { + timeval start, end; + gettimeofday(&start, 0); + Real cfl = ilc.CFL(fields[0]); + if (fields[0].taskid() == 0) + cout << " t == " << t << endl; + if (fields[0].taskid() == 0) + cout << " CFL == " << cfl << endl; + Real ul2n = L2Norm(fields[0]); + Real tl2n = L2Norm(fields[1]); + if (fields[0].taskid() == 0) + cout << " L2Norm(u) == " << ul2n << endl; + if (fields[0].taskid() == 0) + cout << " L2Norm(T) == " << tl2n << endl; + + // Take n steps of length dt + ilc.advance(fields, n); + + gettimeofday(&end, 0); + Real sec = (Real)(end.tv_sec - start.tv_sec); + Real ms = (((Real)end.tv_usec) - ((Real)start.tv_usec)); + Real timeused = sec + ms / 1000000.; + if (fields[0].taskid() == 0) + cout << "duration for this timeunit: " << timeused << endl; + if (t != 0) { + avtime += timeused; + i++; + } + if (fields[0].taskid() == 0) + cout << endl; + } + + if (fields[0].taskid() == 0) { + cout << "Average time/timeunit: " << avtime / i << "s" << endl; + ofstream fout("benchmark_results", ios::app); + fout << "np0 x np1 == " << cfmpi->nproc0() << " x " << cfmpi->nproc1() << endl; + fout << "fftw_flag == " << (fftwmeasure ? "fftw_measure" : "fftw_patient") << endl; + fout << "Average time/timeunit: " << avtime / i << "s" << endl << endl; + fout.close(); + } + // fftw_mpi_gather_wisdom(MPI_COMM_WORLD); + fftw_savewisdom(); + if (saveresult) + fields[0].save(dir + "/uresult.nc"); + } + cfMPI_Finalize(); + return 0; +} diff --git a/modules/ilc/ilc.cpp b/modules/ilc/ilc.cpp new file mode 100644 index 00000000..4fd16e27 --- /dev/null +++ b/modules/ilc/ilc.cpp @@ -0,0 +1,297 @@ +/** + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include "modules/ilc/ilc.h" + +namespace chflow { + +int field2vector_size(const FlowField& u, const FlowField& temp) { + int Kx = u.kxmaxDealiased(); + int Kz = u.kzmaxDealiased(); + int Ny = u.Ny(); + // FIXME: Determine array size + // The original formula was + // int N = 4* ( Kx+2*Kx*Kz+Kz ) * ( Ny-3 ) +2* ( Ny-2 ); + // but I've not been able to twist my head enough to adapt it do distributed FlowFields. + // Since it doesn't take that much time, we now perform the loops twice, once empty to + // determine cfarray sizes and once with the actual data copying. // Tobias + int N = 0; + if (u.taskid() == u.task_coeff(0, 0)) + N += 2 * (Ny - 2); + for (int kx = 1; kx <= Kx; ++kx) + if (u.taskid() == u.task_coeff(u.mx(kx), 0)) + N += 2 * (Ny - 2) + 2 * (Ny - 4); + for (int kz = 1; kz <= Kz; ++kz) + if (u.taskid() == u.task_coeff(0, u.mz(kz))) + N += 2 * (Ny - 2) + 2 * (Ny - 4); + for (int kx = -Kx; kx <= Kx; ++kx) { + if (kx == 0) + continue; + int mx = u.mx(kx); + for (int kz = 1; kz <= Kz; ++kz) { + int mz = u.mz(kz); + if (u.taskid() == u.task_coeff(mx, mz)) { + N += 2 * (Ny - 2) + 2 * (Ny - 4); + } + } + } + + if (temp.taskid() == temp.task_coeff(0, 0)) + N += Ny; + for (int kx = 1; kx <= Kx; ++kx) + if (temp.taskid() == temp.task_coeff(temp.mx(kx), 0)) + N += 2 * Ny; + for (int kz = 1; kz <= Kz; ++kz) + if (temp.taskid() == temp.task_coeff(0, temp.mz(kz))) + N += 2 * Ny; + for (int kx = -Kx; kx <= Kx; kx++) { + if (kx == 0) + continue; + for (int kz = 1; kz <= Kz; kz++) { + if (temp.taskid() == temp.task_coeff(temp.mx(kx), temp.mz(kz))) { + N += 2 * Ny; + } + } + } + + return N; +} + +/** \brief Turn the two flowfields for velocity and temperature into one Eigen vector + * \param[in] u velocity field + * \param[in] temp temperature field + * \param[in] a vector for the linear algebra + * + * The vectorization of u is analog to the field2vector, temparture is piped entirely + * into the vector (a single independent dimensions) + */ +void field2vector(const FlowField& u, const FlowField& temp, Eigen::VectorXd& a) { + Eigen::VectorXd b; + assert(temp.xzstate() == Spectral && temp.ystate() == Spectral); + field2vector(u, b); + int Kx = temp.kxmaxDealiased(); + int Kz = temp.kzmaxDealiased(); + int Ny = temp.Ny(); + + int n = field2vector_size(u, temp); // b.size() +6*Kx*Kz*Ny; + + if (a.size() < n) + a.resize(n, true); + setToZero(a); + int pos = b.size(); + a.topRows(pos) = b; + + // (0,:,0) + for (int ny = 0; ny < Ny; ++ny) + if (temp.taskid() == temp.task_coeff(0, 0)) + a(pos++) = Re(temp.cmplx(0, ny, 0, 0)); + + for (int kx = 1; kx <= Kx; ++kx) { + int mx = temp.mx(kx); + if (temp.taskid() == temp.task_coeff(mx, 0)) { + for (int ny = 0; ny < Ny; ++ny) { + a(pos++) = Re(temp.cmplx(mx, ny, 0, 0)); + a(pos++) = Im(temp.cmplx(mx, ny, 0, 0)); + } + } + } + for (int kz = 1; kz <= Kz; ++kz) { + int mz = temp.mz(kz); + if (temp.taskid() == temp.task_coeff(0, mz)) { + for (int ny = 0; ny < Ny; ++ny) { + a(pos++) = Re(temp.cmplx(0, ny, mz, 0)); + a(pos++) = Im(temp.cmplx(0, ny, mz, 0)); + } + } + } + for (int kx = -Kx; kx <= Kx; kx++) { + if (kx == 0) + continue; + int mx = temp.mx(kx); + for (int kz = 1; kz <= Kz; kz++) { + int mz = temp.mz(kz); + if (u.taskid() == u.task_coeff(mx, mz)) { + for (int ny = 0; ny < Ny; ny++) { + a(pos++) = Re(temp.cmplx(mx, ny, mz, 0)); + a(pos++) = Im(temp.cmplx(mx, ny, mz, 0)); + } + } + } + } +} + +/** \brief Turn one Eigen vector into the two flowfields for velocity and temperature + * \param[in] a vector for the linear algebra + * \param[in] u velocity field + * \param[in] temp temperature field + * + * The vectorization of u is analog to the field2vector, temperature is piped entirely + * into the vector (a single independent dimension) + */ +void vector2field(const Eigen::VectorXd& a, FlowField& u, FlowField& temp) { + assert(temp.xzstate() == Spectral && temp.ystate() == Spectral); + temp.setToZero(); + Eigen::VectorXd b; + int N = field2vector_size(u); + int Kx = temp.kxmaxDealiased(); + int Kz = temp.kzmaxDealiased(); + int Ny = temp.Ny(); + b = a.topRows(N); + vector2field(b, u); + + int pos = N; + double reval, imval; + Complex val; + + if (temp.taskid() == temp.task_coeff(0, 0)) + for (int ny = 0; ny < Ny; ++ny) + temp.cmplx(0, ny, 0, 0) = Complex(a(pos++), 0); + + for (int kx = 1; kx <= Kx; ++kx) { + int mx = temp.mx(kx); + if (temp.taskid() == temp.task_coeff(mx, 0)) { + for (int ny = 0; ny < Ny; ++ny) { + reval = a(pos++); + imval = a(pos++); + temp.cmplx(mx, ny, 0, 0) = Complex(reval, imval); + } + } + + // ------------------------------------------------------ + // Now copy conjugates of u(kx,ny,0,i) to u(-kx,ny,0,i). These are + // redundant modes stored only for the convenience of FFTW. + int mxm = temp.mx(-kx); + int send_id = temp.task_coeff(mx, 0); + int rec_id = temp.task_coeff(mxm, 0); + for (int ny = 0; ny < Ny; ++ny) { + if (temp.taskid() == send_id && send_id == rec_id) { // all is on the same process -> just copy + temp.cmplx(mxm, ny, 0, 0) = conj(temp.cmplx(mx, ny, 0, 0)); + } +#ifdef HAVE_MPI // send_id != rec_id requires multiple processes + else { // Transfer the conjugates via MPI + if (temp.taskid() == send_id) { + Complex tmp0 = conj(temp.cmplx(mx, ny, 0, 0)); + MPI_Send(&tmp0, 1, MPI_DOUBLE_COMPLEX, rec_id, 0, MPI_COMM_WORLD); + } + if (u.taskid() == rec_id) { + Complex tmp0; + MPI_Status status; + MPI_Recv(&tmp0, 1, MPI_DOUBLE_COMPLEX, send_id, 0, MPI_COMM_WORLD, &status); + temp.cmplx(mxm, ny, 0, 0) = tmp0; + } + } +#endif + } + } + for (int kz = 1; kz <= Kz; ++kz) { + int mz = temp.mz(kz); + if (temp.taskid() == temp.task_coeff(0, mz)) { + for (int ny = 0; ny < Ny; ++ny) { + reval = a(pos++); + imval = a(pos++); + temp.cmplx(0, ny, mz, 0) = Complex(reval, imval); + } + } + } + for (int kx = -Kx; kx <= Kx; kx++) { + if (kx == 0) + continue; + int mx = temp.mx(kx); + for (int kz = 1; kz <= Kz; kz++) { + int mz = temp.mz(kz); + if (u.taskid() == u.task_coeff(mx, mz)) { + for (int ny = 0; ny < Ny; ny++) { + reval = a(pos++); + imval = a(pos++); + val = Complex(reval, imval); + temp.cmplx(mx, ny, mz, 0) = val; + } + } + } + } + temp.setPadded(true); +} + +// ILC::ILC() +// : +// DNS(){ +// +// } + +ILC::ILC(const std::vector& fields, const ILCFlags& flags) + : // base class constructor with no arguments is called automatically (see DNS::DNS()) + main_obe_(0), + init_obe_(0) { + main_obe_ = newOBE(fields, flags); + // creates DNSAlgo with ptr of "nse"-daughter type "obe" + main_algorithm_ = newAlgorithm(fields, main_obe_, flags); + if (!main_algorithm_->full() && flags.initstepping != flags.timestepping) { + ILCFlags initflags = flags; + initflags.timestepping = flags.initstepping; + initflags.dt = flags.dt; + init_obe_ = newOBE(fields, flags); + + // creates DNSAlgo with ptr of "nse"-daughter type "obe" + init_algorithm_ = newAlgorithm(fields, init_obe_, initflags); + // Safety check + if (init_algorithm_->Ninitsteps() != 0) + std::cerr << "DNS::DNS(fields, flags) :\n" + << flags.initstepping << " can't initialize " << flags.timestepping + << " since it needs initialization itself.\n"; + } +} + +ILC::~ILC() {} + +const ChebyCoeff& ILC::Tbase() const { + if (main_obe_) + return main_obe_->Tbase(); + else if (init_obe_) + return init_obe_->Tbase(); + else { + std::cerr << "Error in ILC::Tbase():Tbase is currently undefined" << std::endl; + exit(1); + return init_obe_->Tbase(); // to make compiler happy + } +} + +const ChebyCoeff& ILC::Ubase() const { + if (main_obe_) + return main_obe_->Ubase(); + else if (init_obe_) + return init_obe_->Ubase(); + else { + std::cerr << "Error in ILC::Ubase(): Ubase is currently undefined" << std::endl; + exit(1); + return init_obe_->Ubase(); // to make compiler happy + } +} + +const ChebyCoeff& ILC::Wbase() const { + if (main_obe_) + return main_obe_->Wbase(); + else if (init_obe_) + return init_obe_->Wbase(); + else { + std::cerr << "Error in ILC::Wbase(): Wbase is currently undefined" << std::endl; + exit(1); + return init_obe_->Wbase(); // to make compiler happy + } +} + +std::shared_ptr ILC::newOBE(const std::vector& fields, const ILCFlags& flags) { + std::shared_ptr obe(new OBE(fields, flags)); + return obe; +} + +std::shared_ptr ILC::newOBE(const std::vector& fields, const std::vector& base, + const ILCFlags& flags) { + std::shared_ptr obe(new OBE(fields, base, flags)); + return obe; +} + +} // namespace chflow \ No newline at end of file diff --git a/modules/ilc/ilc.h b/modules/ilc/ilc.h new file mode 100644 index 00000000..6efd0f91 --- /dev/null +++ b/modules/ilc/ilc.h @@ -0,0 +1,82 @@ +/** + * Main interface to handle DNS of ILC + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#ifndef ILC_H +#define ILC_H + +#include "channelflow/dns.h" +#include "channelflow/dnsalgo.h" +#include "modules/ilc/ilcflags.h" +#include "modules/ilc/obe.h" + +namespace chflow { + +int field2vector_size(const FlowField& u, const FlowField& temp); + +/** \brief Turn the two flowfields for velocity and temperature into one Eigen vector + * \param[in] u velocity field + * \param[in] temp temperature field + * \param[in] x vector for the linear algebra + * + * The vectorization of u is analog to the field2vector, temparture is piped entirely + * into the vector (a single independent dimensions) + */ +void field2vector(const FlowField& u, const FlowField& temp, Eigen::VectorXd& x); + +/** \brief Turn one Eigen vector into the two flowfields for velocity and temperature + * \param[in] x vector for the linear algebra + * \param[in] u velocity field + * \param[in] temp temperature field + * + * The vectorization of u is analog to the field2vector, temperature is piped entirely + * into the vector (a single independent dimension) + */ +void vector2field(const Eigen::VectorXd& x, FlowField& u, FlowField& temp); + +/** \brief extension of the DNSFlags class + * + * VEDNSFlags class, holds all additional parameters for viscoelastic fluids + */ + +/** \brief wrapper class of DNSAlgorithm and ILC + * + * + */ +class ILC : public DNS { + public: + // ILC (); + // ILC (const ILC & ilc); + ILC(const std::vector& fields, const ILCFlags& flags); + // ILC (const vector & fields, const vector & base, + // const ILCFlags & flags); + + virtual ~ILC(); + + ILC& operator=(const ILC& ilc); + + // virtual void advance (vector & fields, int nSteps = 1); + // + // virtual void reset_dt (Real dt); + // virtual void printStack () const; + + const ChebyCoeff& Tbase() const; + const ChebyCoeff& Ubase() const; + const ChebyCoeff& Wbase() const; + + protected: + std::shared_ptr main_obe_; + std::shared_ptr init_obe_; + + std::shared_ptr newOBE(const std::vector& fields, const ILCFlags& flags); + std::shared_ptr newOBE(const std::vector& fields, const std::vector& base, + const ILCFlags& flags); +}; + +} // namespace chflow +#endif // ILC_H \ No newline at end of file diff --git a/modules/ilc/ilcdsi.cpp b/modules/ilc/ilcdsi.cpp new file mode 100644 index 00000000..56e8093e --- /dev/null +++ b/modules/ilc/ilcdsi.cpp @@ -0,0 +1,1266 @@ +/** + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include "modules/ilc/ilcdsi.h" +#ifdef HAVE_MPI +#include +#endif + +#include "cfbasics/mathdefs.h" +#include "channelflow/diffops.h" +// #include "viscoelastic/veutils.h" +#include "modules/ilc/ilc.h" + +using namespace std; + +namespace chflow { + +/*utility functions*/ + +std::vector ilcstats(const FlowField& u, const FlowField& temp, const ILCFlags flags) { + double l2n = L2Norm(u); + if (std::isnan(l2n)) { + cferror("L2Norm(u) is nan"); + } + + FlowField u_tot = totalVelocity(u, flags); + FlowField temp_tot = totalTemperature(temp, flags); + + std::vector stats; + stats.push_back(L2Norm(u)); + stats.push_back(heatcontent(temp_tot, flags)); + stats.push_back(L2Norm(u_tot)); + stats.push_back(L2Norm(temp_tot)); + stats.push_back(L2Norm3d(u)); + stats.push_back(Ecf(u)); + stats.push_back(getUbulk(u)); + stats.push_back(getWbulk(u)); + stats.push_back(wallshear(u_tot)); + stats.push_back(buoyPowerInput(u_tot, temp_tot, flags)); + stats.push_back(dissipation(u_tot, flags)); + stats.push_back(heatinflux(temp_tot, flags)); + stats.push_back(L2Norm(temp)); + stats.push_back(Nusselt_plane(u_tot, temp_tot, flags)); + return stats; +} + +string ilcfieldstatsheader(const ILCFlags flags) { + stringstream header; + header << setw(14) << "L2(u')" << setw(10) << "(y=" << flags.ystats << ")" // change position with L2(T') + << setw(14) << "L2(u)" << setw(14) << "L2(T)" << setw(14) << "e3d" << setw(14) << "ecf" << setw(14) + << "ubulk" << setw(14) << "wbulk" << setw(14) << "wallshear" << setw(14) << "buoyPowIn" << setw(14) + << "totalDiss" << setw(14) << "heatinflux" << setw(14) << "L2(T')" << setw(10) << "Nu(y=" << flags.ystats + << ")"; + return header.str(); +} + +string ilcfieldstatsheader_t(const string tname, const ILCFlags flags) { + stringstream header; + header << setw(6) << "#(" << tname << ")" << ilcfieldstatsheader(flags); + return header.str(); +} + +string ilcfieldstats(const FlowField& u, const FlowField& temp, const ILCFlags flags) { + std::vector stats = ilcstats(u, temp, flags); + // Return string + stringstream s; + for (uint i = 0; i < stats.size(); i++) { + s << setw(14) << stats[i]; + } + return s.str(); +} + +string ilcfieldstats_t(const FlowField& u, const FlowField& temp, const Real t, const ILCFlags flags) { + std::vector stats = ilcstats(u, temp, flags); + // Return string + stringstream s; + s << setw(8) << t; + for (uint i = 0; i < stats.size(); i++) { + s << setw(14) << stats[i]; + } + return s.str(); +} + +FlowField totalVelocity(const FlowField& velo, const ILCFlags flags) { + // copy + FlowField u(velo); + FlowField tmp(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b(), u.cfmpi()); + + // get base flow + ILC ilc({u, tmp, tmp}, flags); + ChebyCoeff Ubase = ilc.Ubase(); + ChebyCoeff Wbase = ilc.Wbase(); + + // add base flow (should be identical to code in function temperatureNL in OBE + for (int ny = 0; ny < u.Ny(); ++ny) { + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, ny, 0, 0) += Complex(Ubase(ny), 0.0); + u.cmplx(0, ny, 0, 2) += Complex(Wbase(ny), 0.0); + } + } + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, 0, 0, 1) -= Complex(flags.Vsuck, 0.); + } + + return u; +} + +FlowField totalTemperature(const FlowField& temp, const ILCFlags flags) { + // copy + FlowField T(temp); + FlowField tmp(T.Nx(), T.Ny(), T.Nz(), 3, T.Lx(), T.Lz(), T.a(), T.b(), T.cfmpi()); + + // get base flow + ILC ilc({tmp, T, T}, flags); + ChebyCoeff Tbase = ilc.Tbase(); + + // add base flow (should be identical to code in function temperatureNL in OBE + for (int ny = 0; ny < T.Ny(); ++ny) { + if (T.taskid() == T.task_coeff(0, 0)) + T.cmplx(0, ny, 0, 0) += Complex(Tbase(ny), 0.0); + } + + return T; +} + +Real buoyPowerInput(const FlowField& utot, const FlowField& ttot, const ILCFlags flags, bool relative) { + // calculate the bouyancy force along the velocity field to get + // the power input to the kinetic energy equation + + // get parameters + Real nu = flags.nu; + Real sing = sin(flags.gammax); + Real cosg = cos(flags.gammax); + Real grav = flags.grav; + Real laminarInput = grav * sing * sing / (720 * nu); // normalized by Volume + + // prepare loop over field + FlowField u(utot); + FlowField T(ttot); + FlowField xinput(T.Nx(), T.Ny(), T.Nz(), T.Nd(), T.Lx(), T.Lz(), T.a(), T.b(), T.cfmpi(), Physical, Physical); + FlowField yinput(T.Nx(), T.Ny(), T.Nz(), T.Nd(), T.Lx(), T.Lz(), T.a(), T.b(), T.cfmpi(), Physical, Physical); + lint Nz = u.Nz(); + lint nxlocmin = u.nxlocmin(); + lint nxlocmax = u.nxlocmin() + u.Nxloc(); + lint nylocmin = u.nylocmin(); + lint nylocmax = u.nylocmax(); + + // sum up buoyancy term + u.makePhysical(); + T.makePhysical(); +#ifdef HAVE_MPI + for (lint nx = nxlocmin; nx < nxlocmax; ++nx) + for (lint nz = 0; nz < Nz; ++nz) + for (lint ny = nylocmin; ny < nylocmax; ++ny) { + xinput(nx, ny, nz, 0) = u(nx, ny, nz, 0) * T(nx, ny, nz, 0); + yinput(nx, ny, nz, 0) = u(nx, ny, nz, 1) * T(nx, ny, nz, 0); + } +#else + for (lint ny = nylocmin; ny < nylocmax; ++ny) + for (lint nx = nxlocmin; nx < nxlocmax; ++nx) + for (lint nz = 0; nz < Nz; ++nz) { + xinput(nx, ny, nz, 0) = u(nx, ny, nz, 0) * T(nx, ny, nz, 0); + yinput(nx, ny, nz, 0) = u(nx, ny, nz, 1) * T(nx, ny, nz, 0); + } +#endif + xinput.makeSpectral(); + yinput.makeSpectral(); + + // calculate the input mean with cheby profile (code is taken from OBE::initConstraint) + ChebyCoeff xprof(T.My(), T.a(), T.b(), Spectral); + ChebyCoeff yprof(T.My(), T.a(), T.b(), Spectral); + Real xtmp = 0; + Real ytmp = 0; + for (int ny = 0; ny < T.My(); ++ny) { + if (utot.taskid() == 0) { + xtmp = sing * Re(xinput.cmplx(0, ny, 0, 0)); + ytmp = cosg * Re(yinput.cmplx(0, ny, 0, 0)); + } +#ifdef HAVE_MPI + MPI_Bcast(&xtmp, 1, MPI_DOUBLE, utot.task_coeff(0, 0), *utot.comm_world()); + MPI_Bcast(&ytmp, 1, MPI_DOUBLE, utot.task_coeff(0, 0), *utot.comm_world()); +#endif + xprof[ny] = xtmp; + yprof[ny] = ytmp; + } + + Real buoyancyInput = xprof.mean() + yprof.mean(); + // difference between full input and laminar input + if (relative && abs(laminarInput) > 1e-12) { + buoyancyInput *= 1.0 / laminarInput; + buoyancyInput -= 1; + } + + return buoyancyInput; +} + +Real dissipation(const FlowField& utot, const ILCFlags flags, bool normalize, bool relative) { + Real diss = flags.nu * dissipation(utot, normalize); + + // analytic laminar dissipation (for standard base flow) + Real sing = sin(flags.gammax); + Real laminarDiss = flags.grav * sing * sing / (720 * flags.nu); // normalized by Volume + // difference between full input and laminar input + if (relative && abs(laminarDiss) > 1e-12) { + diss *= 1.0 / laminarDiss; + diss -= 1; + } + + return diss; +} + +Real heatinflux(const FlowField& temp, const ILCFlags flags, bool normalize, bool relative) { + // with reference to wallshear, but only lower wall + assert(temp.ystate() == Spectral); + Real I = 0; + if (temp.taskid() == temp.task_coeff(0, 0)) { + ChebyCoeff tprof = Re(temp.profile(0, 0, 0)); + ChebyCoeff dTdy = diff(tprof); + I = flags.kappa * abs(dTdy.eval_a()); + } +#ifdef HAVE_MPI + MPI_Bcast(&I, 1, MPI_DOUBLE, temp.task_coeff(0, 0), temp.cfmpi()->comm_world); +#endif + if (!normalize) + I *= 2 * temp.Lx() * temp.Lz(); + if (relative) { + Real kappa = flags.kappa; + Real deltaT = flags.tlowerwall - flags.tupperwall; + Real H = temp.b() - temp.a(); + I *= 1.0 / kappa / deltaT * H; + I -= 1.0; + } + return I; +} + +Real heatcontent(const FlowField& ttot, const ILCFlags flags) { + assert(ttot.ystate() == Spectral); + int N = 100; + Real dy = (flags.ystats - ttot.a()) / (N - 1); + Real avt = 0; // average temperature + if (ttot.taskid() == ttot.task_coeff(0, 0)) { + ChebyCoeff tprof = Re(ttot.profile(0, 0, 0)); + for (int i = 0; i < N; ++i) { + Real y = ttot.a() + i * dy; + avt += tprof.eval(y); + } + avt *= 1.0 / N; + } +#ifdef HAVE_MPI + MPI_Bcast(&avt, 1, MPI_DOUBLE, ttot.task_coeff(0, 0), ttot.cfmpi()->comm_world); +#endif + return avt; +} + +Real Nusselt_plane(const FlowField& utot, const FlowField& ttot, const ILCFlags flags, bool relative) { + assert(utot.ystate() == Spectral && ttot.ystate() == Spectral); + + // calculate product for advective heat transport + FlowField u(utot); + FlowField T(ttot); + FlowField vt(T.Nx(), T.Ny(), T.Nz(), T.Nd(), T.Lx(), T.Lz(), T.a(), T.b(), T.cfmpi(), Physical, Physical); + lint Nz = u.Nz(); + lint nxlocmin = u.nxlocmin(); + lint nxlocmax = u.nxlocmin() + u.Nxloc(); + lint nylocmin = u.nylocmin(); + lint nylocmax = u.nylocmax(); + + // loop to form product + u.makePhysical(); + T.makePhysical(); +#ifdef HAVE_MPI + for (lint nx = nxlocmin; nx < nxlocmax; ++nx) + for (lint nz = 0; nz < Nz; ++nz) + for (lint ny = nylocmin; ny < nylocmax; ++ny) { + vt(nx, ny, nz, 0) = u(nx, ny, nz, 1) * T(nx, ny, nz, 0); + } +#else + for (lint ny = nylocmin; ny < nylocmax; ++ny) + for (lint nx = nxlocmin; nx < nxlocmax; ++nx) + for (lint nz = 0; nz < Nz; ++nz) { + vt(nx, ny, nz, 0) = u(nx, ny, nz, 1) * T(nx, ny, nz, 0); + } +#endif + vt.makeSpectral(); + + // calculate Nusselt number + Real Nu = 0; + Real y = flags.ystats; + Real kappa = flags.kappa; + if (ttot.taskid() == ttot.task_coeff(0, 0)) { + ChebyCoeff tprof = Re(ttot.profile(0, 0, 0)); + ChebyCoeff dTdy = diff(tprof); + ChebyCoeff vtprof = Re(vt.profile(0, 0, 0)); + Nu = vtprof.eval(y) - + kappa * dTdy.eval(y); // Formula 10 in Chilla&Schumacher 2012 (together with normalization below) + } + +#ifdef HAVE_MPI + MPI_Bcast(&Nu, 1, MPI_DOUBLE, ttot.task_coeff(0, 0), ttot.cfmpi()->comm_world); +#endif + + if (relative) { + Real deltaT = flags.tlowerwall - flags.tupperwall; + Real H = utot.b() - utot.a(); + Nu *= 1.0 / kappa / deltaT * H; + Nu -= 1.0; + } + return Nu; +} + +/* Begin of ilcDSI class*/ + +ilcDSI::ilcDSI() {} + +ilcDSI::ilcDSI(ILCFlags& ilcflags, FieldSymmetry sigma, PoincareCondition* h, TimeStep dt, bool Tsearch, bool xrelative, + bool zrelative, bool Tnormalize, Real Unormalize, const FlowField& u, const FlowField& temp, ostream* os) + : cfDSI(ilcflags, sigma, h, dt, Tsearch, xrelative, zrelative, Tnormalize, Unormalize, u, os), + ilcflags_(ilcflags) {} + +Eigen::VectorXd ilcDSI::eval(const Eigen::VectorXd& x) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + + Real T; + extractVectorILC(x, u, temp, sigma_, T); + + FlowField Gu(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField Gtemp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + G(u, temp, T, h_, sigma_, Gu, Gtemp, ilcflags_, dt_, Tnormalize_, Unormalize_, fcount_, CFL_, *os_); + Eigen::VectorXd Gx(Eigen::VectorXd::Zero(x.rows())); + // Galpha *= 1./vednsflags_.b_para; + field2vector(Gu, Gtemp, Gx); // This does not change the size of Gx and automatically leaves the last entries zero + + return Gx; +} + +Eigen::VectorXd ilcDSI::eval(const Eigen::VectorXd& x0, const Eigen::VectorXd& x1, bool symopt) { + FlowField u0(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField u1(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp0(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp1(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + Real T0, T1; + FieldSymmetry sigma0, sigma1; + extractVectorILC(x0, u0, temp0, sigma0, T0); + extractVectorILC(x1, u1, temp1, sigma1, T1); + + FlowField Gu(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField Gtemp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + + f(u0, temp0, T0, h_, Gu, Gtemp, ilcflags_, dt_, fcount_, CFL_, *os_); + if (symopt) { + Gu *= sigma0; + if (sigma0.sy() == -1) { + // wall-normal mirroring in velocity requires sign change in temperature + FieldSymmetry inv(-1); + sigma0 *= inv; + } + Gtemp *= sigma0; + } + Gu -= u1; + Gtemp -= temp1; + + // normalize + if (Tnormalize_) { + Gu *= 1.0 / T0; + Gtemp *= 1.0 / T0; + } + if (Unormalize_ != 0.0) { + Real funorm = L2Norm3d(Gu); + Gu *= 1. / sqrt(abs(funorm * (Unormalize_ - funorm))); + // u should stay off zero, so normalize with u for now - temp should also stay away from zero + Gtemp *= 1. / sqrt(abs(funorm * (Unormalize_ - funorm))); + } + + Eigen::VectorXd Gx(Eigen::VectorXd::Zero(x0.rows())); + field2vector(Gu, Gtemp, Gx); // This does not change the size of Gx and automatically leaves the last entries zero + + return Gx; +} + +void ilcDSI::save(const Eigen::VectorXd& x, const string filebase, const string outdir, const bool fieldsonly) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + + u.save(outdir + "u" + filebase); + temp.save(outdir + "t" + filebase); + + if (!fieldsonly) { + string fs = ilcfieldstats(u, temp, ilcflags_); + if (u.taskid() == 0) { + if (xrelative_ || zrelative_ || !sigma.isIdentity()) + sigma.save(outdir + "sigma" + filebase); + if (Tsearch_) + chflow::save(T, outdir + "T" + filebase); + // sigma.save (outdir+"sigmaconverge.asc", ios::app); + ofstream fout((outdir + "fieldconverge.asc").c_str(), ios::app); + long pos = fout.tellp(); + if (pos == 0) + fout << ilcfieldstatsheader() << endl; + fout << fs << endl; + fout.close(); + ilcflags_.save(outdir); + } + } +} + +string ilcDSI::stats(const Eigen::VectorXd& x) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + return ilcfieldstats_t(u, temp, mu_, ilcflags_); +} + +pair ilcDSI::stats_minmax(const Eigen::VectorXd& x) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField Gu(u); + FlowField Gtemp(temp); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + + std::vector stats = ilcstats(u, temp, ilcflags_); + std::vector minstats(stats); + std::vector maxstats(stats); + + // quick hack to avoid new interface or creating simple f() for ILC + TimeStep dt = TimeStep(ilcflags_.dt, 0, 1, 1, 0, 0, false); + int fcount = 0; + PoincareCondition* h = 0; + Real CFL = 0.0; + std::ostream muted_os(0); + Real timep = T / 100.0; + + *os_ << "Using flag -orbOut: Calculate minmax-statistics of periodic orbit." << endl; + for (int t = 0; t < 100; t++) { + f(u, temp, timep, h, Gu, Gtemp, ilcflags_, dt, fcount, CFL, muted_os); + stats = ilcstats(Gu, Gtemp, ilcflags_); + for (uint i = 0; i < stats.size(); i++) { + minstats[i] = (minstats[i] < stats[i]) ? minstats[i] : stats[i]; + maxstats[i] = (maxstats[i] > stats[i]) ? maxstats[i] : stats[i]; + } + u = Gu; + temp = Gtemp; + } + // Return string + stringstream smin; + stringstream smax; + smin << setw(8) << mu_; + smax << setw(8) << mu_; + for (uint i = 0; i < stats.size(); i++) { + smin << setw(14) << minstats[i]; + smax << setw(14) << maxstats[i]; + } + + pair minmax; + minmax = make_pair(smin.str(), smax.str()); + return minmax; +} + +string ilcDSI::statsHeader() { return ilcfieldstatsheader_t(ilc_cPar2s(ilc_cPar_), ilcflags_); } + +/// after finding new solution fix phases +void ilcDSI::phaseShift(Eigen::VectorXd& x) { + if (xphasehack_ || zphasehack_) { + FlowField unew(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField tnew(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, unew, tnew, sigma, T); + // vector2field (x,unew); + const int phasehackcoord = 0; // Those values were fixed in continuesoln anyway + const parity phasehackparity = Odd; + const Real phasehackguess = 0.0; + + if (zphasehack_) { + FieldSymmetry tau = zfixphasehack(unew, phasehackguess, phasehackcoord, phasehackparity); + cout << "fixing z phase of potential solution with phase shift tau == " << tau << endl; + unew *= tau; + tnew *= tau; + } + if (xphasehack_) { + FieldSymmetry tau = xfixphasehack(unew, phasehackguess, phasehackcoord, phasehackparity); + cout << "fixing x phase of potential solution with phase shift tau == " << tau << endl; + unew *= tau; + tnew *= tau; + } + if (uUbasehack_) { + cout << "fixing u+Ubase decomposition so that = 0 at walls (i.e. Ubase balances mean pressure " + "gradient))" + << endl; + Real ubulk = Re(unew.profile(0, 0, 0)).mean(); + if (abs(ubulk) < 1e-15) + ubulk = 0.0; + + ChebyCoeff Ubase = laminarProfile(ilcflags_.nu, ilcflags_.constraint, ilcflags_.dPdx, + ilcflags_.Ubulk - ubulk, ilcflags_.Vsuck, unew.a(), unew.b(), + ilcflags_.ulowerwall, ilcflags_.uupperwall, unew.Ny()); + + fixuUbasehack(unew, Ubase); + } + makeVectorILC(unew, tnew, sigma, T, x); + } +} + +void ilcDSI::phaseShift(Eigen::MatrixXd& y) { + if (xphasehack_ || zphasehack_) { + FlowField unew(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField tnew(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + Eigen::VectorXd yvec; + FieldSymmetry sigma; + Real T; + + const int phasehackcoord = 0; // Those values were fixed in continuesoln anyway + const parity phasehackparity = Odd; + const Real phasehackguess = 0.0; + + FieldSymmetry taux(0.0, 0.0); + FieldSymmetry tauz(0.0, 0.0); + + extractVectorILC(y.col(0), unew, tnew, sigma, T); + + if (xphasehack_) { + taux = xfixphasehack(unew, phasehackguess, phasehackcoord, phasehackparity); + cout << "fixing x phase of potential solution with phase shift tau == " << taux << endl; + } + if (zphasehack_) { + tauz = zfixphasehack(unew, phasehackguess, phasehackcoord, phasehackparity); + cout << "fixing z phase of potential solution with phase shift tau == " << tauz << endl; + } + + for (int i = 0; i < y.cols(); i++) { + extractVectorILC(y.col(i), unew, tnew, sigma, T); + unew *= taux; + tnew *= taux; + unew *= tauz; + tnew *= tauz; + makeVectorILC(unew, tnew, sigma, T, yvec); + y.col(i) = yvec; + } + } +} + +Real ilcDSI::extractT(const Eigen::VectorXd& x) { // inefficient hack + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + return T; +} + +Real ilcDSI::extractXshift(const Eigen::VectorXd& x) { // inefficient hack + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + return sigma.ax(); +} + +Real ilcDSI::extractZshift(const Eigen::VectorXd& x) { // inefficient hack + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + extractVectorILC(x, u, temp, sigma, T); + return sigma.az(); +} + +void ilcDSI::makeVectorILC(const FlowField& u, const FlowField& temp, const FieldSymmetry& sigma, const Real T, + Eigen::VectorXd& x) { + if (u.Nd() != 3) + cferror("ilcDSI::makeVector(): u.Nd() = " + i2s(u.Nd()) + " != 3"); + if (temp.Nd() != 1) + cferror("ilcDSI::makeVector(): temp.Nd() = " + i2s(temp.Nd()) + " != 1"); + int taskid = u.taskid(); + + int uunk = field2vector_size(u, temp); // # of variables for u and alpha unknonwn + const int Tunk = (Tsearch_ && taskid == 0) ? uunk : -1; // index for T unknown + const int xunk = (xrelative_ && taskid == 0) ? uunk + Tsearch_ : -1; + const int zunk = (zrelative_ && taskid == 0) ? uunk + Tsearch_ + xrelative_ : -1; + int Nunk = (taskid == 0) ? uunk + Tsearch_ + xrelative_ + zrelative_ : uunk; + if (x.rows() < Nunk) + x.resize(Nunk); + field2vector(u, temp, x); + if (taskid == 0) { + if (Tsearch_) + x(Tunk) = T; + if (xrelative_) + x(xunk) = sigma.ax(); + if (zrelative_) + x(zunk) = sigma.az(); + } +} + +void ilcDSI::extractVectorILC(const Eigen::VectorXd& x, FlowField& u, FlowField& temp, FieldSymmetry& sigma, Real& T) { + int uunk = field2vector_size(u, temp); // number of components in x that corresond to u and alpha + vector2field(x, u, temp); + const int Tunk = uunk + Tsearch_ - 1; + const int xunk = uunk + Tsearch_ + xrelative_ - 1; + const int zunk = uunk + Tsearch_ + xrelative_ + zrelative_ - 1; + Real ax = 0; + Real az = 0; + if (u.taskid() == 0) { + T = Tsearch_ ? x(Tunk) : Tinit_; + ax = xrelative_ ? x(xunk) : axinit_; + az = zrelative_ ? x(zunk) : azinit_; + } +#ifdef HAVE_MPI + MPI_Bcast(&T, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&ax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&az, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); +#endif + sigma = FieldSymmetry(sigma_.sx(), sigma_.sy(), sigma_.sz(), ax, az, sigma_.s()); +} + +Eigen::VectorXd ilcDSI::xdiff(const Eigen::VectorXd& a) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + vector2field(a, u, temp); + Eigen::VectorXd dadx(a.size()); + dadx.setZero(); + u = chflow::xdiff(u); + temp = chflow::xdiff(temp); + field2vector(u, temp, dadx); + dadx *= 1. / L2Norm(dadx); + return dadx; +} + +Eigen::VectorXd ilcDSI::zdiff(const Eigen::VectorXd& a) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + vector2field(a, u, temp); + Eigen::VectorXd dadz(a.size()); + dadz.setZero(); + u = chflow::zdiff(u); + temp = chflow::zdiff(temp); + field2vector(u, temp, dadz); + dadz *= 1. / L2Norm(dadz); + return dadz; +} + +Eigen::VectorXd ilcDSI::tdiff(const Eigen::VectorXd& a, Real epsDt) { + FlowField u(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField temp(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FieldSymmetry sigma; + Real T; + // quick hack to avoid new interface or creating simple f() for ILC + TimeStep dt = TimeStep(epsDt, 0, 1, 1, 0, 0, false); + int fcount = 0; + PoincareCondition* h = 0; + Real CFL = 0.0; + FlowField edudtf(u); + FlowField edtempdtf(temp); + std::ostream muted_os(0); + // vector2field (a, u, temp); + extractVectorILC(a, u, temp, sigma, T); + // use existing f() instead of simple + f(u, temp, epsDt, h, edudtf, edtempdtf, ilcflags_, dt, fcount, CFL, muted_os); + // f (temp, 1,epsDt, edtempdtf, dnsflags_, *os_); + edudtf -= u; + edtempdtf -= temp; + Eigen::VectorXd dadt(a.size()); + field2vector(edudtf, edtempdtf, dadt); + dadt *= 1. / L2Norm(dadt); + return dadt; +} + +void ilcDSI::updateMu(Real mu) { + DSI::updateMu(mu); + if (ilc_cPar_ == ilc_continuationParameter::none) { + cfDSI::updateMu(mu); + } else if (ilc_cPar_ == ilc_continuationParameter::Ra) { + Real Prandtl = ilcflags_.nu / ilcflags_.kappa; + ilcflags_.nu = sqrt(Prandtl / mu); + ilcflags_.kappa = 1.0 / sqrt(Prandtl * mu); + } else if (ilc_cPar_ == ilc_continuationParameter::Pr) { + Real Rayleigh = 1 / (ilcflags_.nu * ilcflags_.kappa); + ilcflags_.nu = sqrt(mu / Rayleigh); + ilcflags_.kappa = 1.0 / sqrt(mu * Rayleigh); + } else if (ilc_cPar_ == ilc_continuationParameter::gx) { + ilcflags_.gammax = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::gz) { + ilcflags_.gammaz = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::gxEps) { + Real gx = ilcflags_.gammax; + Real Ra_old = 1 / (ilcflags_.nu * ilcflags_.kappa); + Real Prandtl = ilcflags_.nu / ilcflags_.kappa; + Real Rc = 1707.76; + Real Rc2 = 8053.1; + Real gc2 = 77.7567; + // get old threshold + Real Rac_old = 0; + if (gx / pi * 180.0 < gc2) + Rac_old = Rc / cos(gx); + else + Rac_old = 1.0 / 41.0 * pow(gx / pi * 180.0 - gc2, 3) + 5.0 / 14.0 * pow(gx / pi * 180.0 - gc2, 2) + + 29 * (gx / pi * 180.0 - gc2) + Rc2; + // get new threshold + Real Rac_new = 0; + if (mu / pi * 180.0 < gc2) + Rac_new = Rc / cos(mu); + else + Rac_new = 1.0 / 41.0 * pow(mu / pi * 180.0 - gc2, 3) + 5.0 / 14.0 * pow(mu / pi * 180.0 - gc2, 2) + + 29 * (mu / pi * 180.0 - gc2) + Rc2; + Real Ra_new = Ra_old / Rac_old * Rac_new; + cout << "Normalized Rayleigh number is epsilon = " << Ra_new / Rac_new - 1 << endl; + ilcflags_.nu = sqrt(Prandtl / Ra_new); + ilcflags_.kappa = 1.0 / sqrt(Prandtl * Ra_new); + ilcflags_.gammax = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Grav) { + ilcflags_.grav = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Tref) { + ilcflags_.t_ref = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::P) { + ilcflags_.dPdx = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Uw) { + ilcflags_.Uwall = mu; + ilcflags_.ulowerwall = -mu * cos(ilcflags_.theta); + ilcflags_.uupperwall = mu * cos(ilcflags_.theta); + ilcflags_.wlowerwall = -mu * sin(ilcflags_.theta); + ilcflags_.wupperwall = mu * sin(ilcflags_.theta); + ; + } else if (ilc_cPar_ == ilc_continuationParameter::UwGrav) { + ilcflags_.Uwall = mu; + ilcflags_.ulowerwall = -mu * cos(ilcflags_.theta); + ilcflags_.uupperwall = mu * cos(ilcflags_.theta); + ilcflags_.wlowerwall = -mu * sin(ilcflags_.theta); + ilcflags_.wupperwall = mu * sin(ilcflags_.theta); + ; + ilcflags_.grav = 1 - mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Rot) { + ilcflags_.rotation = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Theta) { + ilcflags_.theta = mu; + ilcflags_.ulowerwall = -ilcflags_.Uwall * cos(mu); + ilcflags_.uupperwall = ilcflags_.Uwall * cos(mu); + ilcflags_.wlowerwall = -ilcflags_.Uwall * sin(mu); + ilcflags_.wupperwall = ilcflags_.Uwall * sin(mu); + ; + } else if (ilc_cPar_ == ilc_continuationParameter::ThArc) { + Real xleg = Lz_ / tan(ilcflags_.theta); // hypothetical Lx + Lz_ *= sin(mu) / sin(ilcflags_.theta); // rescale Lz for new angle at const diagonal (of xleg x Lz_) + Lx_ *= Lz_ / tan(mu) / xleg; + ilcflags_.theta = mu; + ilcflags_.ulowerwall = -ilcflags_.Uwall * cos(mu); + ilcflags_.uupperwall = ilcflags_.Uwall * cos(mu); + ilcflags_.wlowerwall = -ilcflags_.Uwall * sin(mu); + ilcflags_.wupperwall = ilcflags_.Uwall * sin(mu); + ; + } else if (ilc_cPar_ == ilc_continuationParameter::ThLx) { + Lx_ *= tan(ilcflags_.theta) / tan(mu); + ilcflags_.theta = mu; + ilcflags_.ulowerwall = -ilcflags_.Uwall * cos(mu); + ilcflags_.uupperwall = ilcflags_.Uwall * cos(mu); + ilcflags_.wlowerwall = -ilcflags_.Uwall * sin(mu); + ilcflags_.wupperwall = ilcflags_.Uwall * sin(mu); + ; + } else if (ilc_cPar_ == ilc_continuationParameter::ThLz) { + Lz_ *= tan(mu) / tan(ilcflags_.theta); + ilcflags_.theta = mu; + ilcflags_.ulowerwall = -ilcflags_.Uwall * cos(mu); + ilcflags_.uupperwall = ilcflags_.Uwall * cos(mu); + ilcflags_.wlowerwall = -ilcflags_.Uwall * sin(mu); + ilcflags_.wupperwall = ilcflags_.Uwall * sin(mu); + ; + } else if (ilc_cPar_ == ilc_continuationParameter::Lx) { + Lx_ = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Lz) { + Lz_ = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Aspect) { + Real aspect = Lx_ / Lz_; + Real update = mu / aspect; + + // do only half the adjustment in x, the other in z (i.e. equivalent to Lx_new = 0.5Lx +0.5Lx * update) + Lx_ *= 1 + (update - 1) / 2.0; + Lz_ = Lx_ / mu; + } else if (ilc_cPar_ == ilc_continuationParameter::Diag) { + Real aspect = Lx_ / Lz_; + Real theta = atan(aspect); + Real diagonal = sqrt(Lx_ * Lx_ + Lz_ * Lz_); + Real update = mu - diagonal; + // Lx_ += sqrt ( (update*update * aspect*aspect) / (1 + aspect*aspect)); + Lx_ += update * sin(theta); + Lz_ = Lx_ / aspect; + } else if (ilc_cPar_ == ilc_continuationParameter::Vs) { + ilcflags_.Vsuck = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::VsNu) { + ilcflags_.nu = mu; + ilcflags_.Vsuck = mu; + } else if (ilc_cPar_ == ilc_continuationParameter::VsH) { + Real nu = ilcflags_.nu; + Real Vs = nu * (1 - exp(-mu)); + Real H = nu * mu / Vs; + ya_ = 0; + yb_ = H; + ilcflags_.Vsuck = Vs; + } else if (ilc_cPar_ == ilc_continuationParameter::H) { + ya_ = 0; + yb_ = mu; + } else { + throw invalid_argument("ilcDSI::updateMu(): continuation parameter is unknown"); + } +} + +void ilcDSI::chooseMuILC(string muName) { + ilc_continuationParameter ilc_cPar = s2ilc_cPar(muName); + + if (ilc_cPar == ilc_continuationParameter::none) + ilcDSI::chooseMu(muName); + else + chooseMuILC(ilc_cPar); +} + +void ilcDSI::chooseMuILC(ilc_continuationParameter mu) { + ilc_cPar_ = mu; + Real Rayleigh; + Real Prandtl; + switch (mu) { + case ilc_continuationParameter::Ra: + Rayleigh = 1 / (ilcflags_.nu * ilcflags_.kappa); + updateMu(Rayleigh); + break; + case ilc_continuationParameter::Pr: + Prandtl = ilcflags_.nu / ilcflags_.kappa; + updateMu(Prandtl); + break; + case ilc_continuationParameter::gx: + updateMu(ilcflags_.gammax); + break; + case ilc_continuationParameter::gz: + updateMu(ilcflags_.gammaz); + break; + case ilc_continuationParameter::gxEps: + updateMu(ilcflags_.gammax); + break; + case ilc_continuationParameter::Grav: + updateMu(ilcflags_.grav); + break; + case ilc_continuationParameter::Tref: + updateMu(ilcflags_.t_ref); + break; + case ilc_continuationParameter::P: + updateMu(ilcflags_.dPdx); + break; + case ilc_continuationParameter::Uw: + updateMu(ilcflags_.uupperwall / cos(ilcflags_.theta)); + break; + case ilc_continuationParameter::UwGrav: + updateMu(ilcflags_.uupperwall / cos(ilcflags_.theta)); + break; + case ilc_continuationParameter::Rot: + updateMu(ilcflags_.rotation); + break; + case ilc_continuationParameter::Theta: + updateMu(ilcflags_.theta); + break; + case ilc_continuationParameter::ThArc: + updateMu(ilcflags_.theta); + break; + case ilc_continuationParameter::ThLx: + updateMu(ilcflags_.theta); + break; + case ilc_continuationParameter::ThLz: + updateMu(ilcflags_.theta); + break; + case ilc_continuationParameter::Lx: + updateMu(Lx_); + break; + case ilc_continuationParameter::Lz: + updateMu(Lz_); + break; + case ilc_continuationParameter::Aspect: + updateMu(Lx_ / Lz_); + break; + case ilc_continuationParameter::Diag: + updateMu(sqrt(Lx_ * Lx_ + Lz_ * Lz_)); + break; + case ilc_continuationParameter::Vs: + updateMu(ilcflags_.Vsuck); + break; + case ilc_continuationParameter::VsNu: + updateMu(ilcflags_.Vsuck); + break; + case ilc_continuationParameter::VsH: + updateMu(-log(1 - ilcflags_.Vsuck / ilcflags_.nu)); + break; + case ilc_continuationParameter::H: + updateMu(yb_ - ya_); + break; + case ilc_continuationParameter::none: + throw invalid_argument( + "ilcDSI::chooseMu(): continuation parameter is none, we should not reach this point"); + default: + throw invalid_argument("ilcDSI::chooseMu(): continuation parameter is unknown"); + } +} + +ilc_continuationParameter ilcDSI::s2ilc_cPar(string muname) { + std::transform(muname.begin(), muname.end(), muname.begin(), ::tolower); // why is the string made lower case? + if (muname == "ra") + return ilc_continuationParameter::Ra; + else if (muname == "pr") + return ilc_continuationParameter::Pr; + else if (muname == "gx") + return ilc_continuationParameter::gx; + else if (muname == "gz") + return ilc_continuationParameter::gz; + else if (muname == "gxeps") + return ilc_continuationParameter::gxEps; + else if (muname == "grav") + return ilc_continuationParameter::Grav; + else if (muname == "tref") + return ilc_continuationParameter::Tref; + else if (muname == "p") + return ilc_continuationParameter::P; + else if (muname == "uw") + return ilc_continuationParameter::Uw; + else if (muname == "uwgrav") + return ilc_continuationParameter::UwGrav; + else if (muname == "rot") + return ilc_continuationParameter::Rot; + else if (muname == "theta") + return ilc_continuationParameter::Theta; + else if (muname == "tharc") + return ilc_continuationParameter::ThArc; + else if (muname == "thlx") + return ilc_continuationParameter::ThLx; + else if (muname == "thlz") + return ilc_continuationParameter::ThLz; + else if (muname == "lx") + return ilc_continuationParameter::Lx; + else if (muname == "lz") + return ilc_continuationParameter::Lz; + else if (muname == "aspect") + return ilc_continuationParameter::Aspect; + else if (muname == "diag") + return ilc_continuationParameter::Diag; + else if (muname == "vs") + return ilc_continuationParameter::Vs; + else if (muname == "vsnu") + return ilc_continuationParameter::VsNu; + else if (muname == "vsh") + return ilc_continuationParameter::VsH; + else if (muname == "h") + return ilc_continuationParameter::H; + else + // cout << "ilcDSI::s2ilc_cPar(): ilc_continuation parameter '"+muname+"' is unknown, defaults to 'none'" << + // endl; + return ilc_continuationParameter::none; +} + +string ilcDSI::printMu() { return ilc_cPar2s(ilc_cPar_); } + +string ilcDSI::ilc_cPar2s(ilc_continuationParameter ilc_cPar) { + if (ilc_cPar == ilc_continuationParameter::none) + return cfDSI::cPar2s(cPar_); + else if (ilc_cPar == ilc_continuationParameter::Ra) + return "Ra"; + else if (ilc_cPar == ilc_continuationParameter::Pr) + return "Pr"; + else if (ilc_cPar == ilc_continuationParameter::gx) + return "gx"; + else if (ilc_cPar == ilc_continuationParameter::gz) + return "gz"; + else if (ilc_cPar == ilc_continuationParameter::gxEps) + return "gxeps"; + else if (ilc_cPar == ilc_continuationParameter::Grav) + return "Grav"; + else if (ilc_cPar == ilc_continuationParameter::Tref) + return "Tref"; + else if (ilc_cPar == ilc_continuationParameter::P) + return "P"; + else if (ilc_cPar == ilc_continuationParameter::Uw) + return "Uw"; + else if (ilc_cPar == ilc_continuationParameter::UwGrav) + return "Uw(Grav)"; + else if (ilc_cPar == ilc_continuationParameter::Rot) + return "Rot"; + else if (ilc_cPar == ilc_continuationParameter::Theta) + return "Theta"; + else if (ilc_cPar == ilc_continuationParameter::ThArc) + return "Theta(DiagArc)"; + else if (ilc_cPar == ilc_continuationParameter::ThLx) + return "Theta(Lx)"; + else if (ilc_cPar == ilc_continuationParameter::ThLz) + return "Theta(Lz)"; + else if (ilc_cPar == ilc_continuationParameter::Lx) + return "Lx"; + else if (ilc_cPar == ilc_continuationParameter::Lz) + return "Lz"; + else if (ilc_cPar == ilc_continuationParameter::Aspect) + return "Aspect"; + else if (ilc_cPar == ilc_continuationParameter::Diag) + return "Diag"; + else if (ilc_cPar == ilc_continuationParameter::Vs) + return "Vs"; + else if (ilc_cPar == ilc_continuationParameter::VsNu) + return "VsNu"; + else if (ilc_cPar == ilc_continuationParameter::VsH) + return "VsH"; + else if (ilc_cPar == ilc_continuationParameter::H) + return "H"; + else + throw invalid_argument("ilcDSI::ilc_cPar2s(): continuation parameter is not convertible to string"); +} + +void ilcDSI::saveParameters(string searchdir) { + // cfDSI::saveParameters (searchdir); + ilcflags_.save(searchdir); +} + +void ilcDSI::saveEigenvec(const Eigen::VectorXd& ev, const string label, const string outdir) { + FlowField efu(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField eft(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + vector2field(ev, efu, eft); + efu *= 1.0 / L2Norm(efu); + eft *= 1.0 / L2Norm(eft); + efu.save(outdir + "efu" + label); + eft.save(outdir + "eft" + label); +} + +void ilcDSI::saveEigenvec(const Eigen::VectorXd& evA, const Eigen::VectorXd& evB, const string label1, + const string label2, const string outdir) { + FlowField efAu(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField efBu(Nx_, Ny_, Nz_, Nd_, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField efAt(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + FlowField efBt(Nx_, Ny_, Nz_, 1, Lx_, Lz_, ya_, yb_, cfmpi_); + vector2field(evA, efAu, efAt); + vector2field(evB, efBu, efBt); + Real cu = 1.0 / sqrt(L2Norm2(efAu) + L2Norm2(efBu)); + Real ct = 1.0 / sqrt(L2Norm2(efAt) + L2Norm2(efBt)); + efAu *= cu; + efBu *= cu; + efAt *= ct; + efBt *= ct; + efAu.save(outdir + "efu" + label1); + efBu.save(outdir + "efu" + label2); + efAt.save(outdir + "eft" + label1); + efBt.save(outdir + "eft" + label2); +} + +/* OUTSIDE CLASS */ + +// G(x) = G(u,sigma) = (sigma f^T(u) - u) for orbits +void G(const FlowField& u, const FlowField& temp, Real& T, PoincareCondition* h, const FieldSymmetry& sigma, + FlowField& Gu, FlowField& Gtemp, const ILCFlags& ilcflags, const TimeStep& dt, bool Tnormalize, Real Unormalize, + int& fcount, Real& CFL, ostream& os) { + f(u, temp, T, h, Gu, Gtemp, ilcflags, dt, fcount, CFL, os); + Real funorm = L2Norm3d(Gu); + Gu *= sigma; + Gu -= u; + if (sigma.sy() == -1) { + // wall-normal mirroring in velocity requires sign change in temperature + FieldSymmetry tsigma(sigma); + FieldSymmetry inv(-1); + tsigma *= inv; + Gtemp *= tsigma; + } else { + Gtemp *= sigma; + } + Gtemp -= temp; + if (Tnormalize) { + Gu *= 1.0 / T; + Gtemp *= 1.0 / T; + } + if (Unormalize != 0.0) { + Gu *= 1. / sqrt(abs(funorm * (Unormalize - funorm))); + // u should stay off zero, so normalize with u for now - temp should also stay away from zero + Gtemp *= 1. / sqrt(abs(funorm * (Unormalize - funorm))); + } +} + +void f(const FlowField& u, const FlowField& temp, Real& T, PoincareCondition* h, FlowField& f_u, FlowField& f_temp, + const ILCFlags& ilcflags_, const TimeStep& dt_, int& fcount, Real& CFL, ostream& os) { + if (!isfinite(L2Norm(u))) { + os << "error in f: u is not finite. exiting." << endl; + exit(1); + } + ILCFlags flags(ilcflags_); + flags.logstream = &os; + TimeStep dt(dt_); + vector fields = {u, temp, FlowField(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b(), u.cfmpi())}; + + // f_u = u; + // f_temp = temp; + // No Poincare section, just integration to time T + if (h == 0) { + if (T < 0) { + os << "f: negative integration time T == " << T << endl + << "returning f(u,T) == (1+abs(T))*u" << endl + << "returning f(temp,T) == (1+abs(T))*temp" << endl; + fields[0] *= 1 + abs(T); + fields[1] *= 1 + abs(T); + return; + } + // Special case #1: no time steps + if (T == 0) { + os << "f: T==0, no integration, returning u and temp" << endl; + return; + } + dt.adjust_for_T(T, false); + flags.dt = dt; + // Adjust dt for CFL if necessary + ILC ilc(fields, flags); + ilc.advance(fields, 1); + if (dt.variable()) { + dt.adjust(ilc.CFL(fields[0]), false); + ilc.reset_dt(dt); + } + // t == current time in integration + // T == total integration time + // dT == CFL-check interval + // dt == DNS time-step + // N == T/dT, n == dT/dt; + // T == N dT, dT == n dt + // t == s dT (s is loop index) + + os << "f^T: " << flush; + for (int s = 1; s <= dt.N(); ++s) { + Real t = s * dt.dT(); + CFL = ilc.CFL(fields[0]); + if (s % 10 == 0) + os << iround(t) << flush; + else if (s % 2 == 0) { + if (CFL < dt.CFLmin()) + os << '<' << flush; + else if (CFL > dt.CFLmax()) + os << '>' << flush; + else + os << '.' << flush; + } + ilc.advance(fields, dt.n()); + if (dt.variable() && dt.adjust(CFL, false)) + ilc.reset_dt(dt); + } + + } + // Poincare section computation: return Poincare crossing nearest to t=T, with Tmin < t < Tmax. + else { + cout << "Poincare sectioning not yet implemented (markd as experimental)." << endl; + exit(1); + /* // Adjust dt for CFL if necessary + DNSPoincare dns (f_u, h, flags); + if (dt.variable()) { + dns.advance (f_u, p, 1); + dt.adjust (dns.CFL()); + dns.reset_dt (dt,u); + f_u = u; + } + // Collect all Poincare crossings between Tfudgemin and Tfudgemax + // If we don't find one in that range, go as far as Tlastchance + Real dTfudge = 1.05; + Real Tfudgemin = lesser (0.90*T, T - dTfudge*dt.dT()); + Real Tfudgemax = Greater (1.02*T, T + dTfudge*dt.dT()); + Real Tlastchance = 10*T; + + vector ucross; + vector tcross; + vector scross; + int s=0; + int crosssign = 0; // look for crossings in either direction + + os << "f^t: " << flush; + + for (Real t=0; t<=Tlastchance; t += dt.dT(), ++s) { + + CFL = dns.CFL(); + + if (s % 10 == 0) os << iround (t); + else if (s % 2 == 0) { + if (CFL > dt.CFLmax()) os << '>'; + else if (CFL < dt.CFLmin()) os << '<'; + else os << '.'; + os << flush; + } + + // Collect any Poincare crossings + bool crossed = dns.advanceToSection (f_u, p, dt.n(), crosssign, Tfudgemin); + if (crossed && t >= Tfudgemin) { + ucross.push_back (dns.ucrossing()); + tcross.push_back (dns.tcrossing()); + scross.push_back (dns.scrossing()); + } + + // If we've found at least one crossing within the fudge range, stop. + // Otherwise continue trying until Tlastchance + if (ucross.size() > 0 && t >= Tfudgemax) + break; + } + + if (ucross.size() <1) { + + os << "\nError in f(u, T, f_u, flags, dt, fcount, CFL, os) :\n"; + os << "the integration did not reach the Poincare section.\n"; + os << "Returning laminar solution and a b.s. value for the crossing time.\n"; + os << "I hope you can do something useful with them." << endl; + f_u.setToZero(); + T = dns.time(); + ++fcount; + return; + } + os << " " << flush; + + // Now select the crossing that is closest to the estimated crossing time + FlowField ubest = ucross[0]; + Real Tbest = tcross[0]; + int sbest = scross[0]; + int nbest = 0; + + for (uint n=1; n +#include +#include +#include +#include +#include "cfbasics/cfvector.h" +#include "channelflow/cfdsi.h" +#include "channelflow/cfmpi.h" +#include "channelflow/chebyshev.h" +#include "channelflow/flowfield.h" +#include "channelflow/periodicfunc.h" +#include "channelflow/symmetry.h" +#include "channelflow/tausolver.h" +#include "channelflow/utilfuncs.h" +#include "modules/ilc/ilc.h" +#include "nsolver/nsolver.h" + +using namespace std; + +namespace chflow { + +enum class ilc_continuationParameter { + none, + Ra, + Pr, + gx, + gz, + gxEps, + Grav, + Tref, + P, + Uw, + UwGrav, + Rot, + Theta, + ThArc, + ThLx, + ThLz, + Lx, + Lz, + Aspect, + Diag, + Vs, + VsNu, + VsH, + H +}; + +// Real GMRESHookstep_vector (FlowField& u, FlowField& alpha, Real& T, FieldSymmetry& sigma, +// PoincareCondition* h, +// const nsolver::hookstepSearchFlags& searchflags, +// DNSFlags& dnsflags, VEDNSFlags& vednsflags, TimeStep& dt, Real& CFL, Real Unormalize); + +std::vector ilcstats(const FlowField& u, const FlowField& temp, const ILCFlags flags = ILCFlags()); + +// header for fieldstats +string ilcfieldstatsheader(const ILCFlags flags = ILCFlags()); + +// header for fieldstats with parameter t +string ilcfieldstatsheader_t(std::string muname, const ILCFlags flags = ILCFlags()); +string ilcfieldstats(const FlowField& u, const FlowField& temp, const ILCFlags flags = ILCFlags()); +string ilcfieldstats_t(const FlowField& u, const FlowField& temp, Real t, const ILCFlags flags = ILCFlags()); +FlowField totalVelocity(const FlowField& velo, const ILCFlags flags); +FlowField totalTemperature(const FlowField& temp, const ILCFlags flags); +Real buoyPowerInput(const FlowField& u, const FlowField& temp, const ILCFlags flags, bool relative = true); +Real dissipation(const FlowField& u, const ILCFlags flags, bool normalize = true, bool relative = true); +Real heatinflux(const FlowField& temp, const ILCFlags flags, bool normalize = true, bool relative = true); +Real heatcontent(const FlowField& ttot, const ILCFlags flags); +Real Nusselt_plane(const FlowField& u, const FlowField& temp, const ILCFlags flags, bool relative = true); + +class ilcDSI : public cfDSI { + public: + /** \brief default constructor */ + ilcDSI(); + virtual ~ilcDSI() {} + + /** \brief Initialize ilcDSI */ + ilcDSI(ILCFlags& ilcflags, FieldSymmetry sigma, PoincareCondition* h, TimeStep dt, bool Tsearch, bool xrelative, + bool zrelative, bool Tnormalize, Real Unormalize, const FlowField& u, const FlowField& temp, + std::ostream* os = &std::cout); + + Eigen::VectorXd eval(const Eigen::VectorXd& x) override; + Eigen::VectorXd eval(const Eigen::VectorXd& x0, const Eigen::VectorXd& x1, + bool symopt) override; // needed for multishooting + void save(const Eigen::VectorXd& x, const string filebase, const string outdir = "./", + const bool fieldsonly = false) override; + + string stats(const Eigen::VectorXd& x) override; + pair stats_minmax(const Eigen::VectorXd& x) override; + string statsHeader() override; + void phaseShift(Eigen::VectorXd& x) override; + void phaseShift(Eigen::MatrixXd& y) override; + Real extractT(const Eigen::VectorXd& x) override; + Real extractXshift(const Eigen::VectorXd& x) override; + Real extractZshift(const Eigen::VectorXd& x) override; + + void makeVectorILC(const FlowField& u, const FlowField& temp, const FieldSymmetry& sigma, const Real T, + Eigen::VectorXd& x); + void extractVectorILC(const Eigen::VectorXd& x, FlowField& u, FlowField& temp, FieldSymmetry& sigma, Real& T); + + /// \name Compute derivatives of the two FlowFields contained in this vector + Eigen::VectorXd xdiff(const Eigen::VectorXd& a) override; + Eigen::VectorXd zdiff(const Eigen::VectorXd& a) override; + Eigen::VectorXd tdiff(const Eigen::VectorXd& a, Real epsDt) override; + /// \name Handle continuation parameter + void updateMu(Real mu) override; + void chooseMuILC(std::string muName); + void chooseMuILC(ilc_continuationParameter mu); + string printMu() override; // document + void saveParameters(string searchdir) override; + ilc_continuationParameter s2ilc_cPar(std::string muName); + string ilc_cPar2s(ilc_continuationParameter cPar); + + // Save real eigenvectors + void saveEigenvec(const Eigen::VectorXd& x, const string label, const string outdir) override; + // Save complex conjugate eigenvectors pair + void saveEigenvec(const Eigen::VectorXd& x1, const Eigen::VectorXd& x2, const string label1, const string label2, + const string outdir) override; + + protected: + ILCFlags ilcflags_; + ilc_continuationParameter ilc_cPar_ = ilc_continuationParameter::none; +}; + +// G(x) = G(u,sigma) = (sigma f^T(u) - u) for orbits +void G(const FlowField& u, const FlowField& temp, Real& T, PoincareCondition* h, const FieldSymmetry& sigma, + FlowField& Gu, FlowField& Gtemp, const ILCFlags& ilcflags, const TimeStep& dt, bool Tnormalize, Real Unormalize, + int& fcount, Real& CFL, ostream& os); +void f(const FlowField& u, const FlowField& temp, Real& T, PoincareCondition* h, FlowField& f_u, FlowField& f_temp, + const ILCFlags& ilcflags_, const TimeStep& dt_, int& fcount, Real& CFL, ostream& os); + +} // namespace chflow + +#endif diff --git a/modules/ilc/ilcflags.cpp b/modules/ilc/ilcflags.cpp new file mode 100644 index 00000000..db8356d7 --- /dev/null +++ b/modules/ilc/ilcflags.cpp @@ -0,0 +1,140 @@ +/** + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include "modules/ilc/ilcflags.h" + +namespace chflow { + +ILCFlags::ILCFlags(Real nu_, Real kappa_, Real alpha_, Real grav_, Real rho_ref_, Real t_ref_, Real gammax_, + Real gammaz_, Real ulowerwall_, Real uupperwall_, Real wlowerwall_, Real wupperwall_, + Real tlowerwall_, Real tupperwall_, Real ystats_) + : kappa(kappa_), + alpha(alpha_), + grav(grav_), + rho_ref(rho_ref_), + t_ref(t_ref_), + gammax(gammax_), + gammaz(gammaz_), + tlowerwall(tlowerwall_), + tupperwall(tupperwall_), + ystats(ystats_) { + nu = nu_; + ulowerwall = ulowerwall_; + uupperwall = uupperwall_; + wlowerwall = wlowerwall_; + wupperwall = wupperwall_; + + // freefall = sqrt(grav*alpha*(tlowerwall-tupperwall)); + // Treference = 0.5*(tlowerwall+tupperwall); +} + +ILCFlags::ILCFlags(ArgList& args, const bool laurette) { + // ILC system parameters + args.section("System parameters"); + const Real Rayleigh_ = args.getreal("-Ra", "--Rayleigh", 4000, "pseudo-Rayleigh number == 1/(nu kp)"); + const Real Prandtl_ = args.getreal("-Pr", "--Prandtl", 1, "Prandtl number == nu/kp"); + const Real gammax_ = args.getreal("-gx", "--gammax", 45, "inclination angle of x-dim in deg"); + const Real gammaz_ = args.getreal("-gz", "--gammaz", 0, "inclination angle of z-dim in deg"); + const Real nuarg_ = + args.getreal("-nu", "--nu", 1, "kinematic viscosity (takes precedence, with kp, over Ra and Pr, if nonunity)"); + const Real kparg_ = args.getreal("-kp", "--kappa", 1, + "thermal conductivity (takes precedence, with nu, over Ra and Pr, if nonunity)"); + + // define Channelflow boundary conditions from arglist + args2BC(args); + + // add ILC boundary conditions + const Real deltaT_ = + args.getreal("-delT", "--DeltaT", 1, "temperature difference between the walls: DeltaT = T_lower - T_upper"); + const Real Tref_ = args.getreal("-Tref", "--Treference", 0, + "reference temperature of the surrounding fluid reservoir for thermal expansion"); + const Real grav_ = args.getreal( + "-grav", "--gravity", 1.0, + "gravity coupling of the velocity fluctuations to the temperature fluctuations (no change of base or Ra)"); + const Real ilcUwall_ = + args.getreal("-ilcUw", "--ilcUwall", 0.0, "in ILC, magnitude of imposed wall velocity, +/-Uwall at y = +/-h"); + + // define Channelflow numerics from arglist + args2numerics(args, laurette); + + // also needed for ILC but better show at the end + const std::string tsymmstr = args.getstr("-tsymms", "--tempsymmetries", "", + "constrain temp(t) to invariant " + "symmetric subspace, argument is the filename for a file " + "listing the generators of the isotropy group"); + const Real ystats_ = args.getreal("-ys", "--ystats", 0, "y-coordinate of height dependent statistics, e.g. Nu(y)"); + + // set flags + nu = (nuarg_ != 1) ? nuarg_ : sqrt(Prandtl_ / Rayleigh_); + kappa = (kparg_ != 1) ? kparg_ : 1.0 / sqrt(Prandtl_ * Rayleigh_); + gammax = gammax_ / 180.0 * pi; + gammaz = gammaz_ / 180.0 * pi; + t_ref = Tref_; + tlowerwall = 0.5 * deltaT_; + tupperwall = -0.5 * deltaT_; + grav = grav_; + ystats = ystats_; + + // overwrite the wall velocity with ilcUwall flag + Uwall = ilcUwall_; + ulowerwall = -ilcUwall_ * cos(theta); + uupperwall = ilcUwall_ * cos(theta); + wlowerwall = -ilcUwall_ * sin(theta); + wupperwall = ilcUwall_ * sin(theta); + + if (tsymmstr.length() > 0) { + SymmetryList tsymms(tsymmstr); + tempsymmetries = tsymms; + } + + save(); +} + +void ILCFlags::save(const std::string& savedir) const { + DNSFlags::save(savedir); + if (mpirank() == 0) { + std::string filename = appendSuffix(savedir, "ilcflags.txt"); + std::ofstream os(filename.c_str()); + if (!os.good()) + cferror("ILCFlags::save(savedir) : can't open file " + filename); + os.precision(16); + os.setf(std::ios::left); + os << std::setw(REAL_IOWIDTH) << kappa << " %kappa\n" + << std::setw(REAL_IOWIDTH) << gammax << " %gammax\n" + << std::setw(REAL_IOWIDTH) << gammaz << " %gammaz\n" + << std::setw(REAL_IOWIDTH) << alpha << " %alpha\n" + << std::setw(REAL_IOWIDTH) << grav << " %grav\n" + << std::setw(REAL_IOWIDTH) << rho_ref << " %rho_ref\n" + << std::setw(REAL_IOWIDTH) << t_ref << " %t_ref\n" + << std::setw(REAL_IOWIDTH) << tupperwall << " %tupperwall\n" + << std::setw(REAL_IOWIDTH) << tlowerwall << " %tlowerwall\n" + << std::setw(REAL_IOWIDTH) << ystats << " %ystats\n"; + os.unsetf(std::ios::left); + } +} + +void ILCFlags::load(int taskid, const std::string indir) { + DNSFlags::load(taskid, indir); + std::ifstream is; + if (taskid == 0) { + is.open(indir + "ilcflags.txt"); + if (!is.good()) + cferror(" ILCFlags::load(taskid, flags, dt, indir): can't open file " + indir + "ilcflags.txt"); + } + kappa = getRealfromLine(taskid, is); + gammax = getRealfromLine(taskid, is); + gammaz = getRealfromLine(taskid, is); + alpha = getRealfromLine(taskid, is); + grav = getRealfromLine(taskid, is); + rho_ref = getRealfromLine(taskid, is); + t_ref = getRealfromLine(taskid, is); + tupperwall = getRealfromLine(taskid, is); + tlowerwall = getRealfromLine(taskid, is); + ystats = getRealfromLine(taskid, is); +} + +} // namespace chflow diff --git a/modules/ilc/ilcflags.h b/modules/ilc/ilcflags.h new file mode 100644 index 00000000..ed95d8c6 --- /dev/null +++ b/modules/ilc/ilcflags.h @@ -0,0 +1,64 @@ +/** + * Control parameters for time-integration within the ILC module + * ILCFlags specifies all relevant parameters for integrating the Oberbeck-Boussinesq equations in doubly periodic + * channel domains. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#ifndef ILCFLAGS_H +#define ILCFLAGS_H + +#include "channelflow/dnsflags.h" +#include "channelflow/utilfuncs.h" + +namespace chflow { + +/** \brief extension of the DNSFlags class for ILC + * + * ILCFlags class, holds all additional parameters for convective shear flows + */ +class ILCFlags : public DNSFlags { + // Is derived from DNSFlags to keep saving and loading consistent + + public: + ILCFlags(Real nu = 0.025, // values are just below primary thermal instability + Real kappa = 0.025, Real alpha = 1.0, Real grav = 1.0, Real rho_ref = 1.0, Real t_ref = 0.0, + Real gammax = 0.0, Real gammaz = 0.0, Real ulowerwall = 0.0, Real uupperwall = 0.0, Real wlowerwall = 0.0, + Real wupperwall = 0.0, Real tlowerwall = 0.5, Real tupperwall = -0.5, Real ystats = 0); + // ILCFlags (const ILCFlags& ilcflags); + // ILCFlags (const DNSFlags& flags); + // ILCFlags (const std::string& filebase); + ILCFlags(ArgList& args, const bool laurette = false); + + /** \brief The infamous virtual destructor */ + virtual ~ILCFlags() = default; + + Real kappa; + Real alpha; + Real grav; + Real rho_ref; + Real t_ref; + Real gammax; + Real gammaz; + + Real tlowerwall; + Real tupperwall; + + Real ystats; + + // Real freefall; + // Real Treference; + + cfarray tempsymmetries; // restrict temp(t) to these symmetries + + // override DNSFlags::save/load. Both methods include a call to the parent method + virtual void save(const std::string& savedir = "") const override; + virtual void load(int taskid, const std::string indir) override; +}; + +} // namespace chflow +#endif // ILCFLAGS_H \ No newline at end of file diff --git a/modules/ilc/obe.cpp b/modules/ilc/obe.cpp new file mode 100644 index 00000000..50c1ecb5 --- /dev/null +++ b/modules/ilc/obe.cpp @@ -0,0 +1,810 @@ +/** + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include "modules/ilc/obe.h" + +namespace chflow { + +void momentumNL(const FlowField& u, const FlowField& T, ChebyCoeff Ubase, ChebyCoeff Wbase, FlowField& f, + FlowField& tmp, ILCFlags flags) { + // compute the nonlinear term of NSE in the usual Channelflow style + navierstokesNL(u, Ubase, Wbase, f, tmp, flags); + + // substract the linear temperature coupling term + // f -= sin(gamma)*T e_x + cos(gamma)*T e_y + Real gammax = flags.gammax; + Real grav = flags.grav; +#ifdef HAVE_MPI + for (int mz = f.mzlocmin(); mz < f.mzlocmin() + f.Mzloc(); mz++) + for (int mx = f.mxlocmin(); mx < f.mxlocmin() + f.Mxloc(); mx++) + for (int ny = 0; ny < f.Ny(); ny++) { + f.cmplx(mx, ny, mz, 0) -= grav * sin(gammax) * T.cmplx(mx, ny, mz, 0); + f.cmplx(mx, ny, mz, 1) -= grav * cos(gammax) * T.cmplx(mx, ny, mz, 0); + } +#else + for (int ny = 0; ny < f.Ny(); ny++) + for (int mx = f.mxlocmin(); mx < f.mxlocmin() + f.Mxloc(); mx++) + for (int mz = f.mzlocmin(); mz < f.mzlocmin() + f.Mzloc(); mz++) { + f.cmplx(mx, ny, mz, 0) -= grav * sin(gammax) * T.cmplx(mx, ny, mz, 0); + f.cmplx(mx, ny, mz, 1) -= grav * cos(gammax) * T.cmplx(mx, ny, mz, 0); + } +#endif + + // dealiasing modes + if (flags.dealias_xz()) + f.zeroPaddedModes(); +} + +void temperatureNL(const FlowField& u_, const FlowField& T_, ChebyCoeff Ubase, ChebyCoeff Wbase, ChebyCoeff Tbase, + FlowField& f, FlowField& tmp, ILCFlags flags) { + FlowField& u = const_cast(u_); + FlowField& T = const_cast(T_); + + // f += Base; + for (int ny = 0; ny < u.Ny(); ++ny) { + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, ny, 0, 0) += Complex(Ubase(ny), 0.0); + u.cmplx(0, ny, 0, 2) += Complex(Wbase(ny), 0.0); + } + if (u.taskid() == u.task_coeff(0, 0)) + T.cmplx(0, ny, 0, 0) += Complex(Tbase(ny), 0.0); + } + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, 0, 0, 1) -= Complex(flags.Vsuck, 0.); + } + + // compute the nonlinearity (temperature advection (u*grad)T ) analogous to the convectiveNL and store in f + dotgradScalar(u, T, f, tmp); + + // f -= Base; + for (int ny = 0; ny < u.Ny(); ++ny) { + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, ny, 0, 0) -= Complex(Ubase(ny), 0.0); + u.cmplx(0, ny, 0, 2) -= Complex(Wbase(ny), 0.0); + } + if (u.taskid() == u.task_coeff(0, 0)) + T.cmplx(0, ny, 0, 0) -= Complex(Tbase(ny), 0.0); + } + if (u.taskid() == u.task_coeff(0, 0)) { + u.cmplx(0, 0, 0, 1) += Complex(flags.Vsuck, 0.); + } + + // dealiasing modes + if (flags.dealias_xz()) + f.zeroPaddedModes(); +} + +OBE::OBE(const std::vector& fields, const ILCFlags& flags) + : NSE(fields, flags), + heatsolver_(0), // heatsolvers are allocated when reset_lambda is called for the first time + flags_(flags), + gsingx_(0), + gcosgx_(0), + Tref_(flags.t_ref), + Tbase_(), + Tbaseyy_(), + Pbasey_(), + Pbasex_(0), + Cu_(), + Cw_(), + Ct_(), + nonzCu_(false), + nonzCw_(false), + nonzCt_(false), + Tk_(Nyd_, a_, b_, Spectral), + Rtk_(Nyd_, a_, b_, Spectral), + baseflow_(false), + constraint_(false) { + gsingx_ = flags.grav * sin(flags_.gammax); + gcosgx_ = flags.grav * cos(flags_.gammax); + + // set member variables for base flow + createILCBaseFlow(); + Pbasex_ = hydrostaticPressureGradientX(flags_); + Pbasey_ = hydrostaticPressureGradientY(Tbase_, flags_); + + // set member variables for contraint + initILCConstraint(fields[0]); + + // define the constant terms of the OBE + createConstants(); +} + +OBE::OBE(const std::vector& fields, const std::vector& base, const ILCFlags& flags) + : NSE(fields, base, flags), + heatsolver_(0), // heatsolvers are allocated when reset_lambda is called for the first time + flags_(flags), + gsingx_(0), + gcosgx_(0), + Tref_(flags.t_ref), + Tbase_(base[2]), + Tbaseyy_(), + Pbasey_(), + Pbasex_(0), + Cu_(), + Cw_(), + Ct_(), + nonzCu_(false), + nonzCw_(false), + nonzCt_(false), + Tk_(Nyd_, a_, b_, Spectral), + Rtk_(Nyd_, a_, b_, Spectral), + baseflow_(false), + constraint_(false) { + gsingx_ = flags.grav * sin(flags_.gammax); + gcosgx_ = flags.grav * cos(flags_.gammax); + + baseflow_ = true; // base flow is passed to constructor + Pbasex_ = hydrostaticPressureGradientX(flags_); + Pbasey_ = hydrostaticPressureGradientY(Tbase_, flags_); + + // set member variables for contraint + initILCConstraint(fields[0]); + + // define the constant terms of the OBE + createConstants(); +} + +OBE::~OBE() { + if (tausolver_) { + for (uint j = 0; j < lambda_t_.size(); ++j) { + for (int mx = 0; mx < Mxloc_; ++mx) { + delete[] tausolver_[j][mx]; // undo new #3 + tausolver_[j][mx] = 0; + } + delete[] tausolver_[j]; // undo new #2 + tausolver_[j] = 0; + } + delete[] tausolver_; // undo new #1 + tausolver_ = 0; + } + + if (heatsolver_) { + for (uint j = 0; j < lambda_t_.size(); ++j) { + for (int mx = 0; mx < Mxloc_; ++mx) { + delete[] heatsolver_[j][mx]; // undo new #3 + heatsolver_[j][mx] = 0; + } + delete[] heatsolver_[j]; // undo new #2 + heatsolver_[j] = 0; + } + delete[] heatsolver_; // undo new #1 + heatsolver_ = 0; + } +} + +void OBE::nonlinear(const std::vector& infields, std::vector& outfields) { + // The first entry in vector must be velocity FlowField, the second a temperature FlowField. + // Pressure as third entry in in/outfields is not touched. + momentumNL(infields[0], infields[1], Ubase_, Wbase_, outfields[0], tmp_, flags_); + temperatureNL(infields[0], infields[1], Ubase_, Wbase_, Tbase_, outfields[1], tmp_, flags_); +} + +void OBE::linear(const std::vector& infields, std::vector& outfields) { + // Method takes input fields {u,press} and computes the linear terms for velocity output field {u} + + // Make sure user does not expect a pressure output. Outfields should be created outside OBE with OBE::createRHS() + assert(infields.size() == (outfields.size() + 1)); + const int kxmax = infields[0].kxmax(); + const int kzmax = infields[0].kzmax(); + const Real nu = flags_.nu; + const Real kappa = flags_.kappa; + + // Loop over Fourier modes. 2nd derivative and summation of linear term + // is most sufficient on ComplexChebyCoeff. Therefore, the old loop structure is kept. + for (lint mx = mxlocmin_; mx < mxlocmin_ + Mxloc_; ++mx) { + const int kx = infields[0].kx(mx); + for (lint mz = mzlocmin_; mz < mzlocmin_ + Mzloc_; ++mz) { + const int kz = infields[0].kz(mz); + + // Skip last and aliased modes + if ((kx == kxmax || kz == kzmax) || (flags_.dealias_xz() && isAliasedMode(kx, kz))) + break; + + // Compute linear temperature terms + //================================ + // Goal is to compute + // L = nu uj" - kappa2 nu uj - grad qj + C + + // Extract relevant Fourier modes of uj and qj + for (int ny = 0; ny < Nyd_; ++ny) { + uk_.set(ny, nu * infields[0].cmplx(mx, ny, mz, 0)); + vk_.set(ny, nu * infields[0].cmplx(mx, ny, mz, 1)); + wk_.set(ny, nu * infields[0].cmplx(mx, ny, mz, 2)); + Pk_.set(ny, infields[2].cmplx(mx, ny, mz, 0)); // pressure is 3. entry in fields + } + + // (1) Put nu uj" into in R. (Pyk_ is used as tmp workspace) + diff2(uk_, Ruk_, Pyk_); + diff2(vk_, Rvk_, Pyk_); + diff2(wk_, Rwk_, Pyk_); + + // (2) Put qn' into Pyk (compute y-comp of pressure gradient). + diff(Pk_, Pyk_); + + // (3) Summation of all derivative terms and assignment to output velocity field. + const Real k2 = 4 * pi * pi * (square(kx / Lx_) + square(kz / Lz_)); + const Complex Dx = infields[0].Dx(mx); + const Complex Dz = infields[0].Dz(mz); + for (int ny = 0; ny < Nyd_; ++ny) { + outfields[0].cmplx(mx, ny, mz, 0) = Ruk_[ny] - k2 * uk_[ny] - Dx * Pk_[ny]; + outfields[0].cmplx(mx, ny, mz, 1) = Rvk_[ny] - k2 * vk_[ny] - Pyk_[ny]; + outfields[0].cmplx(mx, ny, mz, 2) = Rwk_[ny] - k2 * wk_[ny] - Dz * Pk_[ny]; + } + + // (4) Add const. terms + if (kx == 0 && kz == 0) { + // L includes const dissipation term of Ubase and Wbase: nu Uyy, nu Wyy + // and bouyancy term of Temperature profile: sin(gamma)(Tbase_-Tref) + // absolute bouyancy term cancels with Pbasex_ (both constant) + if (nonzCu_ || nonzCw_) { + for (int ny = 0; ny < My_; ++ny) { + outfields[0].cmplx(mx, ny, mz, 0) += Cu_[ny]; + outfields[0].cmplx(mx, ny, mz, 2) += Cw_[ny]; + } + } + // Add base pressure gradient depending on the constraint + if (flags_.constraint == PressureGradient) { + // dPdx is supplied as dPdxRef + outfields[0].cmplx(mx, 0, mz, 0) -= Complex(dPdxRef_, 0); + outfields[0].cmplx(mx, 0, mz, 2) -= Complex(dPdzRef_, 0); + } else { // const bulk velocity + // actual dPdx is unknown but defined by constraint of bulk velocity + // Determine actual dPdx from Ubase + u. + Real Ly = b_ - a_; + diff(uk_, Ruk_); + diff(wk_, Rwk_); + Real dPdxAct = Re(Ruk_.eval_b() - Ruk_.eval_a()) / Ly; + Real dPdzAct = Re(Rwk_.eval_b() - Rwk_.eval_a()) / Ly; + ChebyCoeff Ubasey = diff(Ubase_); + ChebyCoeff Wbasey = diff(Wbase_); + if (Ubase_.length() != 0) + dPdxAct += nu * (Ubasey.eval_b() - Ubasey.eval_a()) / Ly; + if (Wbase_.length() != 0) + dPdzAct += nu * (Wbasey.eval_b() - Wbasey.eval_a()) / Ly; + // add press. gradient to linear term + outfields[0].cmplx(mx, 0, mz, 0) -= Complex(dPdxAct, 0); + outfields[0].cmplx(mx, 0, mz, 2) -= Complex(dPdzAct, 0); + } + } // End of const. terms + + // Compute linear heat equation terms + //================================ + // Goal is to compute + // L = kappa T" - k2 kappa T - Ubase dT/dx + C + + // Extract relevant Fourier modes of T: use Pk and Rxk + for (int ny = 0; ny < Nyd_; ++ny) + Tk_.set(ny, kappa * infields[1].cmplx(mx, ny, mz, 0)); + + // (1) Put kappa T" into in R. (Pyk_ is used as tmp workspace) + diff2(Tk_, Rtk_, Pyk_); + + // (2+3) Summation of both diffusive terms and linear advection of temperature. + // k2 and Dx were defined above for the current Fourier mode + for (int ny = 0; ny < Nyd_; ++ny) + // from nonlinear: - Dx*Tk_[ny]*Complex(Ubase_[ny],0)/kappa; + outfields[1].cmplx(mx, ny, mz, 0) = Rtk_[ny] - k2 * Tk_[ny]; + + // (4) Add const. terms + if (kx == 0 && kz == 0) { + // L includes const diffusion term of Tbase: kappa Tbase_yy + if (nonzCt_) + for (int ny = 0; ny < My_; ++ny) + outfields[1].cmplx(mx, ny, mz, 0) += Ct_[ny]; + } + } + } // End of loop over Fourier modes +} + +void OBE::solve(std::vector& outfields, const std::vector& rhs, const int s) { + // Method takes a right hand side {u} and solves for output fields {u,press} + + // Make sure user provides correct RHS which can be created outside NSE with NSE::createRHS() + assert(outfields.size() == (rhs.size() + 1)); + const int kxmax = outfields[0].kxmax(); + const int kzmax = outfields[0].kzmax(); + + // Update each Fourier mode with solution of the implicit problem + for (lint ix = 0; ix < Mxloc_; ++ix) { + const lint mx = ix + mxlocmin_; + const int kx = outfields[0].kx(mx); + + for (lint iz = 0; iz < Mzloc_; ++iz) { + const lint mz = iz + mzlocmin_; + const int kz = outfields[0].kz(mz); + + // Skip last and aliased modes + if ((kx == kxmax || kz == kzmax) || (flags_.dealias_xz() && isAliasedMode(kx, kz))) + break; + + // Construct ComplexChebyCoeff + for (int ny = 0; ny < Nyd_; ++ny) { + Ruk_.set(ny, rhs[0].cmplx(mx, ny, mz, 0)); + Rvk_.set(ny, rhs[0].cmplx(mx, ny, mz, 1)); + Rwk_.set(ny, rhs[0].cmplx(mx, ny, mz, 2)); + // negative RHS because HelmholtzSolver solves the negative problem + Rtk_.set(ny, -rhs[1].cmplx(mx, ny, mz, 0)); + } + + // Solve the tau equations for momentum + //============================= + if (kx != 0 || kz != 0) + tausolver_[s][ix][iz].solve(uk_, vk_, wk_, Pk_, Ruk_, Rvk_, Rwk_); + // solve(ix,iz,uk_,vk_,wk_,Pk_, Ruk_,Rvk_,Rwk_); + else { // kx,kz == 0,0 + // LHS includes also the constant terms C which can be added to RHS + if (nonzCu_ || nonzCw_) { + for (int ny = 0; ny < My_; ++ny) { + Ruk_.re[ny] += Cu_.re[ny]; + Rwk_.re[ny] += Cw_.re[ny]; + } + } + + if (flags_.constraint == PressureGradient) { + // pressure is supplied, put on RHS of tau eqn + Ruk_.re[0] -= dPdxRef_; + Rwk_.re[0] -= dPdzRef_; + tausolver_[s][ix][iz].solve(uk_, vk_, wk_, Pk_, Ruk_, Rvk_, Rwk_); + // solve(ix,iz,uk_, vk_, wk_, Pk_, Ruk_,Rvk_,Rwk_); + // Bulk vel is free variable determined from soln of tau eqn //TODO: write method that computes + // UbulkAct everytime it is needed + + } else { // const bulk velocity + // bulk velocity is supplied, use alternative tau solver + + // Use tausolver with additional variable and constraint: + // free variable: dPdxAct at next time-step, + // constraint: UbulkBase + mean(u) = UbulkRef. + tausolver_[s][ix][iz].solve(uk_, vk_, wk_, Pk_, dPdxAct_, dPdzAct_, Ruk_, Rvk_, Rwk_, + UbulkRef_ - UbulkBase_, WbulkRef_ - WbulkBase_); + // solve(ix,iz,uk_, vk_, wk_, Pk_, dPdxAct_, dPdzAct_, + // Ruk_, Rvk_, Rwk_, + // UbulkRef_ - UbulkBase_, + // WbulkRef_ - WbulkBase_); + + // test if UbulkRef == UbulkAct = UbulkBase_ + uk_.re.mean() + assert((UbulkRef_ - UbulkBase_ - uk_.re.mean()) < 1e-15); + // test if WbulkRef == WbulkAct = WbulkBase_ + wk_.re.mean() + assert((WbulkRef_ - WbulkBase_ - wk_.re.mean()) < 1e-15); + } + } + // Load solutions into u and p. + // Because of FFTW complex symmetries + // The 0,0 mode must be real. + // For Nx even, the kxmax,0 mode must be real + // For Nz even, the 0,kzmax mode must be real + // For Nx,Nz even, the kxmax,kzmax mode must be real + if ((kx == 0 && kz == 0) || (outfields[0].Nx() % 2 == 0 && kx == kxmax && kz == 0) || + (outfields[0].Nz() % 2 == 0 && kz == kzmax && kx == 0) || + (outfields[0].Nx() % 2 == 0 && outfields[0].Nz() % 2 == 0 && kx == kxmax && kz == kzmax)) { + for (int ny = 0; ny < Nyd_; ++ny) { + outfields[0].cmplx(mx, ny, mz, 0) = Complex(Re(uk_[ny]), 0.0); + outfields[0].cmplx(mx, ny, mz, 1) = Complex(Re(vk_[ny]), 0.0); + outfields[0].cmplx(mx, ny, mz, 2) = Complex(Re(wk_[ny]), 0.0); + outfields[2].cmplx(mx, ny, mz, 0) = Complex(Re(Pk_[ny]), 0.0); + } + } + // The normal case, for general kx,kz + else + for (int ny = 0; ny < Nyd_; ++ny) { + outfields[0].cmplx(mx, ny, mz, 0) = uk_[ny]; + outfields[0].cmplx(mx, ny, mz, 1) = vk_[ny]; + outfields[0].cmplx(mx, ny, mz, 2) = wk_[ny]; + outfields[2].cmplx(mx, ny, mz, 0) = Pk_[ny]; + } + + // Solve the helmholtz problem for the heat equation + //============================= + if (kx != 0 || kz != 0) { + heatsolver_[s][ix][iz].solve(Tk_.re, Rtk_.re, 0.0, 0.0); // BC are considered through the base profile + heatsolver_[s][ix][iz].solve(Tk_.im, Rtk_.im, 0.0, 0.0); + } else { // kx,kz == 0,0 + // LHS includes also the constant term C=kappa Tbase_yy, which can be added to RHS + if (nonzCt_) { + for (int ny = 0; ny < My_; ++ny) + Rtk_.re[ny] -= Ct_.re[ny]; + } + + heatsolver_[s][ix][iz].solve(Tk_.re, Rtk_.re, 0.0, 0.0); + heatsolver_[s][ix][iz].solve(Tk_.im, Rtk_.im, 0.0, 0.0); + } + + // Load solution into T. + // Because of FFTW complex symmetries + // The 0,0 mode must be real. + // For Nx even, the kxmax,0 mode must be real + // For Nz even, the 0,kzmax mode must be real + // For Nx,Nz even, the kxmax,kzmax mode must be real + if ((kx == 0 && kz == 0) || (outfields[1].Nx() % 2 == 0 && kx == kxmax && kz == 0) || + (outfields[1].Nz() % 2 == 0 && kz == kzmax && kx == 0) || + (outfields[1].Nx() % 2 == 0 && outfields[1].Nz() % 2 == 0 && kx == kxmax && kz == kzmax)) { + for (int ny = 0; ny < Nyd_; ++ny) + outfields[1].cmplx(mx, ny, mz, 0) = Complex(Re(Tk_[ny]), 0.0); + + } + // The normal case, for general kx,kz + else + for (int ny = 0; ny < Nyd_; ++ny) + outfields[1].cmplx(mx, ny, mz, 0) = Tk_[ny]; + + } // End of loop over Fourier modes + } +} + +void OBE::reset_lambda(std::vector lambda_t) { + lambda_t_ = lambda_t; + + if ((tausolver_ == 0) && (heatsolver_ == 0)) { // TauSolver need to be constructed + // Allocate memory for [Nsubsteps x Mx_ x Mz_] Tausolver cfarray + tausolver_ = new TauSolver**[lambda_t.size()]; // new #1 + heatsolver_ = new HelmholtzSolver**[lambda_t.size()]; + for (uint j = 0; j < lambda_t.size(); ++j) { + tausolver_[j] = new TauSolver*[Mxloc_]; // new #2 + heatsolver_[j] = new HelmholtzSolver*[Mxloc_]; + for (int mx = 0; mx < Mxloc_; ++mx) { + tausolver_[j][mx] = new TauSolver[Mzloc_]; // new #3 + heatsolver_[j][mx] = new HelmholtzSolver[Mzloc_]; + } + } + } else if ((tausolver_ == 0) && (heatsolver_ != 0)) { + tausolver_ = new TauSolver**[lambda_t.size()]; // new #1 + for (uint j = 0; j < lambda_t.size(); ++j) { + tausolver_[j] = new TauSolver*[Mxloc_]; // new #2 + for (int mx = 0; mx < Mxloc_; ++mx) + tausolver_[j][mx] = new TauSolver[Mzloc_]; // new #3 + } + } else if ((tausolver_ != 0) && (heatsolver_ == 0)) { + heatsolver_ = new HelmholtzSolver**[lambda_t.size()]; // new #1 + for (uint j = 0; j < lambda_t.size(); ++j) { + heatsolver_[j] = new HelmholtzSolver*[Mxloc_]; // new #2 + for (int mx = 0; mx < Mxloc_; ++mx) + heatsolver_[j][mx] = new HelmholtzSolver[Mzloc_]; // new #3 + } + } + + // Configure tausolvers + // FlowField u=fields_[0]; + const Real nu = flags_.nu; + const Real kappa = flags_.kappa; + const Real c = 4.0 * square(pi); + // const int kxmax = u.kxmax(); + // const int kzmax = u.kzmax(); + for (uint j = 0; j < lambda_t.size(); ++j) { + for (int mx = 0; mx < Mxloc_; ++mx) { + int kx = kxloc_[mx]; + for (int mz = 0; mz < Mzloc_; ++mz) { + int kz = kzloc_[mz]; + Real lambda_tau = lambda_t[j] + nu * c * (square(kx / Lx_) + square(kz / Lz_)); + Real lambda_heat = lambda_t[j] + kappa * c * (square(kx / Lx_) + square(kz / Lz_)); + + if ((kx != kxmax_ || kz != kzmax_) && (!flags_.dealias_xz() || !isAliasedMode(kx, kz))) { + tausolver_[j][mx][mz] = + TauSolver(kx, kz, Lx_, Lz_, a_, b_, lambda_tau, nu, Nyd_, flags_.taucorrection); + heatsolver_[j][mx][mz] = HelmholtzSolver(Nyd_, a_, b_, lambda_heat, kappa); + } + } + } + } +} + +const ChebyCoeff& OBE::Tbase() const { return Tbase_; } +const ChebyCoeff& OBE::Ubase() const { return Ubase_; } +const ChebyCoeff& OBE::Wbase() const { return Wbase_; } + +std::vector OBE::createRHS(const std::vector& fields) const { + return {fields[0], fields[1]}; // return only velocity and temperature fields +} + +std::vector> OBE::createSymmVec() const { + // velocity symmetry + cfarray usym = SymmetryList(1); // unity operation + if (flags_.symmetries.length() > 0) + usym = flags_.symmetries; + // temperature symmetry + cfarray tsym = SymmetryList(1); // unity operation + if (flags_.tempsymmetries.length() > 0) + tsym = flags_.tempsymmetries; + // pressure symmetry + cfarray psym = SymmetryList(1); // unity operation + return {usym, tsym, psym}; +} + +void OBE::createILCBaseFlow() { + Real ulowerwall = flags_.ulowerwall; + Real uupperwall = flags_.uupperwall; + Real wlowerwall = flags_.wlowerwall; + Real wupperwall = flags_.wupperwall; + + switch (flags_.baseflow) { + case ZeroBase: + Ubase_ = ChebyCoeff(My_, a_, b_, Spectral); + Wbase_ = ChebyCoeff(My_, a_, b_, Spectral); + Tbase_ = ChebyCoeff(My_, a_, b_, Spectral); + break; + case LinearBase: + std::cerr << "error in OBE::createBaseFlow :\n"; + std::cerr << "LinearBase is not defined in ILC.\n"; + break; + case ParabolicBase: + std::cerr << "error in OBE::createBaseFlow :\n"; + std::cerr << "ParabolicBase is not defined in ILC.\n"; + break; + case SuctionBase: + std::cerr << "error in OBE::createBaseFlow :\n"; + std::cerr << "ParabolicBase is not defined in ILC.\n"; + break; + case LaminarBase: + Ubase_ = laminarVelocityProfile(flags_.gammax, flags_.dPdx, flags_.Ubulk, ulowerwall, uupperwall, a_, b_, + My_, flags_); + + // choosing angle zero removes buoyancy effect + Wbase_ = + laminarVelocityProfile(0.0, flags_.dPdz, flags_.Wbulk, wlowerwall, wupperwall, a_, b_, My_, flags_); + + Tbase_ = linearTemperatureProfile(a_, b_, My_, flags_); + + break; + case ArbitraryBase: + std::cerr << "error in NSE::createBaseFlow :\n"; + std::cerr << "flags.baseflow is ArbitraryBase.\n"; + std::cerr << "Please provide {Ubase, Wbase} when constructing DNS.\n"; + cferror(""); + default: + std::cerr << "error in NSE::createBaseFlow :\n"; + std::cerr << "flags.baseflow should be ZeroBase, LinearBase, ParabolicBase, LaminarBase, SuctionBase.\n"; + std::cerr << "Other cases require use of the DNS::DNS(fields, base, flags) constructor.\n"; + cferror(""); + } + + baseflow_ = true; +} + +void OBE::createConstants() { + if ((!baseflow_) || (!constraint_)) + std::cerr << "OBE::createConstants: Not all base flow members have been initialized." << std::endl; + + ComplexChebyCoeff c(My_, a_, b_, Spectral); + Real nu = flags_.nu; + Real kappa = flags_.kappa; + Real eta = 1.0 / (flags_.alpha * (flags_.tlowerwall - flags_.tupperwall)); + Real t_ref = flags_.t_ref; + + // constant u-term: + Real hydrostatx = gsingx_ * (Tbase_[0] - eta - t_ref) - Pbasex_; + if (abs(hydrostatx) > 1e-15) + *flags_.logstream << "ILC with hydrostatic dPdx." << std::endl; + if (Ubaseyy_.length() > 0) { + c[0] = Complex(nu * Ubaseyy_[0] + hydrostatx, 0); + if ((abs(c[0]) > 1e-15) && !nonzCu_) + nonzCu_ = true; + for (int ny = 1; ny < My_; ++ny) { + c[ny] = Complex(nu * Ubaseyy_[ny] + gsingx_ * Tbase_[ny], 0); + if ((abs(c[ny]) > 1e-15) && !nonzCu_) + nonzCu_ = true; + } + } + if (nonzCu_) { + Cu_ = c; + *flags_.logstream << "ILC with nonzero U-const." << std::endl; + } + + // check that constant v-term is zero: + Real hydrostaty = gcosgx_ * (Tbase_[0] - eta - t_ref) - Pbasey_[0]; + if (abs(hydrostaty) > 1e-15) + std::cerr << "Wall-normal hydrostatic pressure is unballanced" << std::endl; + for (int ny = 1; ny < My_; ++ny) + if (abs(gcosgx_ * Tbase_[ny] - Pbasey_[ny]) > 1e-15) + std::cerr << "Wall-normal hydrostatic pressure is unballanced" << std::endl; + + // constant w-term: + if (Wbaseyy_.length() > 0) { + for (int ny = 0; ny < My_; ++ny) { + c[ny] = Complex(nu * Wbaseyy_[ny], 0); + if ((abs(c[ny]) > 1e-15) && !nonzCw_) + nonzCw_ = true; + } + } + if (nonzCw_) { + Cw_ = c; + *flags_.logstream << "ILC with nonzero W-const." << std::endl; + } + + // constant t-term: + if (Tbaseyy_.length() > 0) { + for (int ny = 1; ny < My_; ++ny) { + c[ny] = Complex(kappa * Tbaseyy_[ny], 0); + if ((abs(c[ny]) > 1e-15) && !nonzCt_) + nonzCt_ = true; + } + } + if (nonzCt_) { + Ct_ = c; + *flags_.logstream << "ILC with nonzero T-const." << std::endl; + } +} + +void OBE::initILCConstraint(const FlowField& u) { + if (!baseflow_) + std::cerr << "OBE::initConstraint: Base flow has not been created." << std::endl; + + // Calculate Ubaseyy_ and related quantities + UbulkBase_ = Ubase_.mean(); + ChebyCoeff Ubasey = diff(Ubase_); + Ubaseyy_ = diff(Ubasey); + WbulkBase_ = Wbase_.mean(); + ChebyCoeff Wbasey = diff(Wbase_); + Wbaseyy_ = diff(Wbasey); + + ChebyCoeff Tbasey = diff(Tbase_); + Tbaseyy_ = diff(Tbasey); + + // Determine actual Ubulk and dPdx from initial data Ubase + u. + Real Ly = b_ - a_; + ChebyCoeff u00(My_, a_, b_, Spectral); + Real reucmplx = 0; + for (int ny = 0; ny < My_; ++ny) { + if (u.taskid() == 0) + reucmplx = Re(u.cmplx(0, ny, 0, 0)); +#ifdef HAVE_MPI + MPI_Bcast(&reucmplx, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); +#endif + u00[ny] = reucmplx; + } + ChebyCoeff du00dy = diff(u00); + UbulkAct_ = UbulkBase_ + u00.mean(); + dPdxAct_ = flags_.nu * (du00dy.eval_b() - du00dy.eval_a()) / Ly; + ChebyCoeff w00(My_, a_, b_, Spectral); + for (int ny = 0; ny < My_; ++ny) { + if (u.taskid() == 0) + reucmplx = Re(u.cmplx(0, ny, 0, 2)); +#ifdef HAVE_MPI + MPI_Bcast(&reucmplx, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); +#endif + w00[ny] = reucmplx; + } + ChebyCoeff dw00dy = diff(w00); + WbulkAct_ = WbulkBase_ + w00.mean(); + dPdzAct_ = flags_.nu * (dw00dy.eval_b() - dw00dy.eval_a()) / Ly; + +#ifdef HAVE_MPI + MPI_Bcast(&dPdxAct_, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); + MPI_Bcast(&UbulkAct_, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); + MPI_Bcast(&dPdzAct_, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); + MPI_Bcast(&WbulkAct_, 1, MPI_DOUBLE, u.task_coeff(0, 0), *u.comm_world()); +#endif + + if (Ubase_.length() != 0) + dPdxAct_ += flags_.nu * (Ubasey.eval_b() - Ubasey.eval_a()) / Ly; + if (Wbase_.length() != 0) + dPdzAct_ += flags_.nu * (Wbasey.eval_b() - Wbasey.eval_a()) / Ly; + + if (flags_.constraint == BulkVelocity) { + UbulkRef_ = flags_.Ubulk; + WbulkRef_ = flags_.Wbulk; + } else { + dPdxAct_ = flags_.dPdx; + dPdxRef_ = flags_.dPdx; + dPdzAct_ = flags_.dPdz; + dPdzRef_ = flags_.dPdz; + } + + constraint_ = true; +} + +ChebyCoeff laminarVelocityProfile(Real gamma, Real dPdx, Real Ubulk, Real Ua, Real Ub, Real a, Real b, int Ny, + ILCFlags flags) { + MeanConstraint constraint = flags.constraint; + Real nu = flags.nu; + Real Vsuck = flags.Vsuck; + Real dT = flags.tlowerwall - flags.tupperwall; + Real mT = 0.5 * (flags.tlowerwall + flags.tupperwall); + Real Tref = flags.t_ref; + + ChebyCoeff u(Ny, a, b, Spectral); + Real grav = flags.grav; + Real q = grav * sin(gamma) / nu / 16.0; + + if (constraint == BulkVelocity) { + cferror("Using ILC with constraint BulkVelocity is not implemented yet"); + } else { + if (Vsuck < 1e-14) { + if (dPdx < 1e-14) { + u[0] = q * (mT - Tref) + 0.5 * (Ub + Ua); // See documentation about base solution + u[1] = -q / 12.0 * dT + 0.5 * (Ub - Ua); + u[2] = -q * (mT - Tref); + u[3] = q / 12.0 * dT; + } else { + cferror("Using ILC with nonzero dPdx is not implemented yet"); + } + } else { + cferror("Using ILC with SuctionVelocity is not implemented yet"); + } + } + return u; +} + +ChebyCoeff linearTemperatureProfile(Real a, Real b, int Ny, ILCFlags flags) { + MeanConstraint constraint = flags.constraint; + + ChebyCoeff T(Ny, a, b, Spectral); + Real Vsuck = flags.Vsuck; + Real dPdx = flags.dPdx; + Real Ta = flags.tlowerwall; + Real Tb = flags.tupperwall; + + if (constraint == BulkVelocity) { + cferror("Using ILC with constraint BulkVelocity is not implemented yet"); + } else { + if (Vsuck < 1e-14) { + if (dPdx < 1e-14) { + T[0] = 0.5 * (Tb + Ta); // the boundary conditions are given in units of Delta_T=Ta-Tb + T[1] = 0.5 * (Tb - Ta); + } else { + cferror("Using ILC with nonzero dPdx is not implemented yet"); + } + } else { + cferror("Using ILC with SuctionVelocity is not implemented yet"); + } + } + return T; +} + +Real hydrostaticPressureGradientX(ILCFlags flags) { + MeanConstraint constraint = flags.constraint; + + Real Px = 0.0; + Real Vsuck = flags.Vsuck; + Real dPdx = flags.dPdx; + Real sgamma = sin(flags.gammax); + Real eta = 1.0 / (flags.alpha * (flags.tlowerwall - flags.tupperwall)); + Real grav = flags.grav; + + if (constraint == BulkVelocity) { + cferror("Using ILC with constraint BulkVelocity is not implemented yet"); + } else { + if (Vsuck < 1e-14) { + if (dPdx < 1e-14) { + Px = -grav * eta * sgamma; + } else { + cferror("Using ILC with nonzero dPdx is not implemented yet"); + } + } else { + cferror("Using ILC with SuctionVelocity is not implemented yet"); + } + } + return Px; +} + +ChebyCoeff hydrostaticPressureGradientY(ChebyCoeff Tbase, ILCFlags flags) { + MeanConstraint constraint = flags.constraint; + + ChebyCoeff Py(Tbase); + Real Vsuck = flags.Vsuck; + Real dPdx = flags.dPdx; + Real cgamma = cos(flags.gammax); + Real eta = 1.0 / (flags.alpha * (flags.tlowerwall - flags.tupperwall)); + Real grav = flags.grav; + + if (constraint == BulkVelocity) { + cferror("Using ILC with constraint BulkVelocity is not implemented yet"); + } else { + if (Vsuck < 1e-14) { + if (dPdx < 1e-14) { + Py[0] = grav * cgamma * (Tbase[0] - eta - flags.t_ref); // TODO: implement general loop + Py[1] = grav * cgamma * Tbase[1]; + } else { + cferror("Using ILC with nonzero dPdx is not implemented yet"); + } + } else { + cferror("Using ILC with SuctionVelocity is not implemented yet"); + } + } + return Py; +} + +} // namespace chflow \ No newline at end of file diff --git a/modules/ilc/obe.h b/modules/ilc/obe.h new file mode 100644 index 00000000..6a141f77 --- /dev/null +++ b/modules/ilc/obe.h @@ -0,0 +1,105 @@ +/** + * System class of Oberbeck-Boussinesq equations for standard channel flows + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#ifndef OBE_H +#define OBE_H + +#include "channelflow/diffops.h" +#include "channelflow/flowfield.h" +#include "channelflow/nse.h" +#include "channelflow/tausolver.h" +#include "modules/ilc/ilcflags.h" + +namespace chflow { + +// nonlinear term of NSE plus the linear coupling term to the temperature equation +void momentumNL(const FlowField& u, const FlowField& T, ChebyCoeff Ubase, ChebyCoeff Wbase, FlowField& f, + FlowField& tmp, ILCFlags flags); + +// nonlinear term of heat equation plus the linear coupling term to the momentum equation +void temperatureNL(const FlowField& u, const FlowField& T, ChebyCoeff Ubase, ChebyCoeff Wbase, ChebyCoeff Tbase, + FlowField& f, FlowField& tmp, ILCFlags flags); + +class OBE : public NSE { + public: + // OBE (); + // OBE (const NSE& nse); + OBE(const std::vector& fields, const ILCFlags& flags); + OBE(const std::vector& fields, const std::vector& base, const ILCFlags& flags); + virtual ~OBE(); + + void nonlinear(const std::vector& infields, std::vector& outfields) override; + void linear(const std::vector& infields, std::vector& outfields) override; + + // calls a tausolver for each Fourier mode + void solve(std::vector& outfields, const std::vector& infields, const int i = 0) override; + + // redefines the tausolver objects with new time-stepping constant (allocates memory for tausolver at first use) + void reset_lambda(const std::vector lambda_t) override; + + // vector of RHS is smaller than of fields because of missing pressure equation + std::vector createRHS(const std::vector& fields) const override; + + // returns vector of symmetries confining the vector of fields to a subspace + std::vector> createSymmVec() const override; + + const ChebyCoeff& Tbase() const; + const ChebyCoeff& Ubase() const override; + const ChebyCoeff& Wbase() const override; + + protected: + HelmholtzSolver*** heatsolver_; // 3d cfarray of tausolvers, indexed by [i][mx][mz] for substep, Fourier Mode x,z + + ILCFlags flags_; // User-defined integration parameters + Real gsingx_; + Real gcosgx_; + Real Tref_; + + // additional base solution profiles + ChebyCoeff Tbase_; // temperature base profile (physical) + ChebyCoeff Tbaseyy_; // 2. deriv. of temperature base profile + ChebyCoeff Pbasey_; // wall normal pressure gradient (y-dependent) + Real Pbasex_; // along wall pressure gradient (const) + + // constant terms + ComplexChebyCoeff Cu_; // constant + ComplexChebyCoeff Cw_; + ComplexChebyCoeff Ct_; + bool nonzCu_; + bool nonzCw_; + bool nonzCt_; + + ComplexChebyCoeff Tk_; + ComplexChebyCoeff Rtk_; + + private: + void createILCBaseFlow(); + void initILCConstraint(const FlowField& u); // method called only at construction + void createConstants(); + + bool baseflow_; + bool constraint_; +}; + +// Construct laminar flow profile for given flow parameters. +// [a,b] == y position of [lower, upper] walls +// [ua,ub] == in-plane speed of [lower, upper] walls +// constraint == is mean pressure gradient fixed, or mean (bulk) velocity? +// dPdx, Ubulk == value of fixed pressure gradient or fixed Ubulk velocity +// Vsuck == suction velocity at walls (asymptotic suction boundary layer) +ChebyCoeff laminarVelocityProfile(Real gamma, Real dPdx, Real Ubulk, Real ua, Real ub, Real a, Real b, int Ny, + ILCFlags flags); + +ChebyCoeff linearTemperatureProfile(Real a, Real b, int Ny, ILCFlags flags); + +Real hydrostaticPressureGradientX(ILCFlags flags); +ChebyCoeff hydrostaticPressureGradientY(ChebyCoeff Tbase, ILCFlags flags); + +} // namespace chflow +#endif \ No newline at end of file diff --git a/modules/ilc/programs/CMakeLists.txt b/modules/ilc/programs/CMakeLists.txt new file mode 100644 index 00000000..9496bb12 --- /dev/null +++ b/modules/ilc/programs/CMakeLists.txt @@ -0,0 +1,15 @@ +set( + ilc_APPS + ilc_simulateflow + ilc_findsoln + ilc_continuesoln + ilc_findeigenvals + ilc_edgetracking + ilc_addbaseflow +) + +foreach (program ${ilc_APPS}) + install_channelflow_application(${program} bin) + target_link_libraries(${program}_app PUBLIC ilc) +endforeach (program) + diff --git a/modules/ilc/programs/ilc_addbaseflow.cpp b/modules/ilc/programs/ilc_addbaseflow.cpp new file mode 100644 index 00000000..70fb278b --- /dev/null +++ b/modules/ilc/programs/ilc_addbaseflow.cpp @@ -0,0 +1,65 @@ +/** + * This program adds a base flow, constructed from input parameters, + * to given velocity and temperature fluctuations. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include +#include +#include +#include + +#include "channelflow/flowfield.h" +#include "modules/ilc/ilc.h" + +using namespace std; +using namespace chflow; + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + string purpose( + "Compute and add the base flow to velocity and " + "temerature fields based on the ILCFlags"); + + ArgList args(argc, argv, purpose); + + ILCFlags flags(args); + TimeStep dt(flags); + + const string uname = args.getstr(4, "", "input field of velocity fluctuations"); + const string tname = args.getstr(3, "", "input field of temperature fluctuations"); + const string ubfname = args.getstr(2, "", "output field of total velocity"); + const string tbfname = args.getstr(1, "", "output field of total temperature"); + + // put all fields into a vector + FlowField u(uname); + FlowField temp(tname); + vector fields = {u, temp, FlowField(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b())}; + + // construct OBE object + OBE obe(fields, flags); + + // + ChebyCoeff Ubase(obe.Ubase()); + ChebyCoeff Wbase(obe.Wbase()); + ChebyCoeff Tbase(obe.Tbase()); + + // add base flow; + for (int ny = 0; ny < u.Ny(); ++ny) { + u.cmplx(0, ny, 0, 0) += Complex(Ubase(ny), 0.0); + u.cmplx(0, ny, 0, 2) += Complex(Wbase(ny), 0.0); + temp.cmplx(0, ny, 0, 0) += Complex(Tbase(ny), 0.0); + } + u.cmplx(0, 0, 0, 1) -= Complex(flags.Vsuck, 0.); + + u.save(ubfname); + temp.save(tbfname); + } + cfMPI_Finalize(); +} diff --git a/modules/ilc/programs/ilc_continuesoln.cpp b/modules/ilc/programs/ilc_continuesoln.cpp new file mode 100644 index 00000000..b554ede5 --- /dev/null +++ b/modules/ilc/programs/ilc_continuesoln.cpp @@ -0,0 +1,308 @@ +/** + * Continuation program, like the program for pure shear flows. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include "cfbasics/cfbasics.h" +#include "channelflow/flowfield.h" +#include "modules/ilc/ilcdsi.h" +#include "nsolver/nsolver.h" + +using namespace std; +using namespace chflow; + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + ArgList args(argc, argv, "parametric continuation of invariant solutions in ILC"); + + ContinuationFlags cflags(args); + cflags.save(); + + unique_ptr N; + bool Rxsearch, Rzsearch, Tsearch; + NewtonSearchFlags searchflags(args); + searchflags.save(); + N = unique_ptr(new NewtonAlgorithm(searchflags)); + + Rxsearch = searchflags.xrelative; + Rzsearch = searchflags.zrelative; + Tsearch = searchflags.solntype == PeriodicOrbit ? true : false; + const bool Tnormalize = Tsearch ? false : true; + + Real diagonal = 0.0; // sqrt(Lx^2 + Lz^2) + int W = 24; + const Real EPSILON = 1e-13; + + ILCFlags ilcflags(args, searchflags.laurette); + TimeStep dt(ilcflags); + ilcflags.save(); + + args.section("Program options"); + const string muname = args.getstr("-cont", "--continuation", "", + "continuation parameter, one of [Ra, Pr, gx, gz, Grav, Tref, P, Uw, UwGrav, " + "Rot, Theta, ThArc, ThLx, ThLz, Lx, Lz, Aspect, Diag, Vs, VsNu, VsH, H]"); + const string sigmastr = + args.getstr("-sigma", "--sigma", "", "file containing sigma of sigma f^T(u) - u = 0 (default == identity)"); + const Real Unormalize = args.getreal("-un", "--unormalize", 0.0, "lower bound in energy for search"); + Real Lxtarg = args.getreal("-Lxtarg", "--LxTarget", 2 * pi, "with -cont Ltarg, aim for this value of Lx"); + Real Lztarg = args.getreal("-Lztarg", "--LzTarget", pi, "with -cont Ltarg, aim for this value of Lz"); + const bool xphasehack = + args.getflag("-xphhack", "--xphasehack", + "fix x phase so that u(x,0,Lz/2,0) - mean = 0 at x=Lx/2 (HACK for u UNSYMM in x)"); + const bool zphasehack = + args.getflag("-zphhack", "--zphasehack", + "fix z phase so that d/dz _{xy}(z) = 0 at z=Lz/2 (HACK for u UNSYMM in z)"); + const bool uUbasehack = + args.getflag("-uUbhack", "--uUbasehack", "fix u/Ubase decomposition so that == 0 at walls."); + const int nproc0 = + args.getint("-np0", "--nproc0", 0, "number of MPI-processes for transpose/number of parallel ffts"); + const int nproc1 = args.getint("-np1", "--nproc1", 0, "number of MPI-processes for one fft"); + + // check if invariant solution is relative + FieldSymmetry givenSigma; + if (sigmastr.length() != 0) + givenSigma = FieldSymmetry(sigmastr); + bool relative = Rxsearch || Rzsearch || !givenSigma.isIdentity(); + + bool restart = cflags.restartMode; + + string uname(""), tname(""), restartdir[3]; + if (restart) { + bool solutionsAvail = readContinuationInfo(restartdir, cflags); + + if (!solutionsAvail) { + restartdir[0] = args.getpath(1, "", "directory containing solution 1"); + restartdir[1] = args.getpath(2, "", "directory containing solution 2"); + restartdir[2] = args.getpath(3, "", "directory containing solution 3"); + } + + } else { + uname = + args.getstr(2, "", "velocity field of initial solution from which to start continuation"); + tname = + args.getstr(1, "", "temperature field of initial solution from which to start continuation"); + } + args.check(); + + if (muname == "") { + cerr << "Please choose --continuation with one of [Ra, Pr, gx, gz, Grav, Tref, P, Uw, UwGrav, Rot, Theta, " + "ThArc, ThLx, ThLz, Lx, Lz, Aspect, Diag, Vs, VsNu, VsH, H]" + << endl; + exit(1); + } + args.save(); + WriteProcessInfo(argc, argv); + + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + + FlowField u[3]; + FlowField temp[3]; + FieldSymmetry sigma[3]; + cfarray mu(3); + Real T[3]; + unique_ptr dsi; + dsi = unique_ptr(new ilcDSI()); + + if (restart) { + cout << "Restarting from previous solutions. Please be aware that DNSFlags and ILCFlags " + << "from the corresponding directories will overwrite any specified command line parameters!" << endl; + for (int i = 0; i < 3; ++i) { + u[i] = FlowField(restartdir[i] + "ubest", cfmpi); + temp[i] = FlowField(restartdir[i] + "tbest", cfmpi); + if (i == 0) { + printout("Optimizing FFTW...", false); + u[i].optimizeFFTW(FFTW_PATIENT); + fftw_savewisdom(); + printout("done"); + u[i] = FlowField(restartdir[i] + "ubest", cfmpi); + } + if (relative) { + sigma[i] = FieldSymmetry(restartdir[i] + "sigmabest.asc"); + } + load(mu[i], restartdir[i] + "mu.asc"); + if (Tsearch) { + load(T[i], restartdir[i] + "Tbest.asc"); + } else { + T[i] = ilcflags.T; + } + } + cout << "Loading flags from " + restartdir[0] + "dnsflags.txt and " + restartdir[0] + + "ilcflags.txt, neglecting command line switches." + << endl; + ilcflags.load(cfmpi->taskid(), restartdir[0]); + dt = TimeStep(ilcflags); + + // Some consistency checks of the initial solutions for a continuation in restart mode + if (!relative && (sigma[0] != sigma[1] || sigma[1] != sigma[2])) + cferror("loadSolutions error : initial symmetries should be equal for non-relative searchs"); + if (muname == "Aspect") { + // Check that diagonal of u[2]'s is same as u[0]. Rescale Lx,Lz keeping aspect ratio + // constant if diagonals are different. Normally reloaded aspect-continued data will + // have constant diagonals, but do check and fix in order to restart data with some + // small fuzz in diagonal. + diagonal = pythag(u[1].Lx(), u[1].Lz()); + for (int n = 2; n >= 0; n -= 2) { + if (abs(diagonal - pythag(u[n].Lx(), u[n].Lz())) >= EPSILON) { + cout << "Diagonal of u[" << n << "] != diagonal of u[1]. Rescaling Lx,Lz...." << endl; + Real alpha_n = atan(1 / mu[n]); + u[n].rescale(diagonal * cos(alpha_n), diagonal * sin(alpha_n)); + temp[n].rescale(diagonal * cos(alpha_n), diagonal * sin(alpha_n)); + } + } + } else if (muname == "Diag") { + // Check that aspect ratio of u[2]'s is same u[0]. Rescale Lx,Lz keeping diagonal + // constant if aspect ratios differ. Normally reloaded diagonl-continued data will + // have constant aspect ratio, but do check and fix in order to restart data with some + // small fuzz in aspect ratio. + Real aspect = u[1].Lx() / u[1].Lz(); + Real alpha = atan(1.0 / (u[1].Lx() / u[1].Lz())); + for (int n = 2; n >= 0; n -= 2) { + if (abs(aspect - u[n].Lx() / u[n].Lz()) >= EPSILON) { + cout << "Aspect ratio of u[" << n << "] != aspect ratio of u[1]. Rescaling Lx,Lz...." << endl; + Real diagonal_n = pythag(u[n].Lx(), u[n].Lz()); + u[n].rescale(diagonal_n * cos(alpha), diagonal_n * sin(alpha)); + temp[n].rescale(diagonal_n * cos(alpha), diagonal_n * sin(alpha)); + } + } + } else if (muname == "Lt") { + // Check that all solutions are colinear with (u[1].Lx, u[1].Lz) and (Lxtarg,Lztarg). + // If not, rescale u[2] and u[0] Lx,Lz so that they are + Real phi = atan((Lztarg - u[1].Lz()) / (Lxtarg - u[1].Lx())); + for (int i = 2; i >= 0; i -= 2) { + if (abs(atan((Lztarg - u[i].Lz()) / (Lxtarg - u[i].Lx())) - phi) >= EPSILON) { + cout << "u[" << i + << "] is not colinear with (u[1].Lx, u[1].Lz) and (Lxtarg,Lztarg). Rescaling to fix." + << endl; + u[i].rescale(Lxtarg - mu[i] * cos(phi), Lztarg - mu[i] * sin(phi)); + temp[i].rescale(Lxtarg - mu[i] * cos(phi), Lztarg - mu[i] * sin(phi)); + } + } + } + cout << endl << "loaded the following data..." << endl; + } else { // not a restart + // Compute initial data points for extrapolation from perturbations of given solution + u[1] = FlowField(uname, cfmpi); + temp[1] = FlowField(tname, cfmpi); + project(ilcflags.symmetries, u[1], "initial value u", cout); + project(ilcflags.tempsymmetries, temp[1], "initial value temp", cout); + fixdivnoslip(u[1]); + + u[2] = u[1]; + printout("Optimizing FFTW...", false); + u[2].optimizeFFTW(FFTW_PATIENT); // Overwrites u[2] + printout("done"); + fftw_savewisdom(); + u[2] = u[1]; + temp[2] = temp[1]; + u[0] = u[1]; + temp[0] = temp[1]; + + if (sigmastr.length() != 0) + sigma[0] = sigma[1] = sigma[2] = givenSigma; + T[0] = T[1] = T[2] = ilcflags.T; + + // begin superfluous output + Real phi, Ltarget; + if (muname == "Lt") { + phi = atan((Lztarg - u[1].Lz()) / (Lxtarg - u[1].Lx())); + Ltarget = spythag(u[1].Lx(), Lxtarg, u[1].Lz(), Lztarg, phi); + } else { + Lxtarg = u[1].Lx(); + Lztarg = u[1].Lz(); + phi = 0.0; + Ltarget = 0.0; + } + cout << "Some geometrical parameters. Phi and Ltarg are used for continuation along line in 2d (Lx,Lz) " + "plane." + << endl; + cout << " phi == " << phi << endl; + cout << "Ltarg == " << Ltarget << endl; + cout << "distance == " << pythag(Lxtarg - u[1].Lx(), Lztarg - u[1].Lz()) << endl; + cout << "u[1].Lx() == " << u[1].Lx() << endl; + cout << " Lxtarg == " << Lxtarg << endl; + cout << "u[1].Lz() == " << u[1].Lz() << endl; + cout << " Lztarg == " << Lztarg << endl; + + cout << "Lxtarg - Ltarg cos(phi) == " << (Lxtarg - Ltarget * cos(phi)) << endl; + cout << "Lztarg - Ltarg sin(phi) == " << (Lztarg - Ltarget * sin(phi)) << endl; + cout << "dt_Lx == " << dt.dt() / u[1].Lx() << endl; + + cout << "set up the following initial data..." << endl; + } + cout << setw(4) << "i" << setw(W) << "T" << setw(W) << setw(W) << "L2Norm(u)" << setw(W) << "sigma" << endl; + for (int i = 2; i >= 0; --i) { + Real l2normui = L2Norm(u[i]); + cout << setw(4) << i << setw(W) << T[i] << setw(W) << l2normui << setw(W) << sigma[i] << setw(W) << endl; + } + // end superfluous output + + dsi = unique_ptr(new ilcDSI(ilcflags, sigma[0], 0, dt, Tsearch, Rxsearch, Rzsearch, Tnormalize, + Unormalize, u[0], temp[0], N->getLogstream())); + + cout << setprecision(8); + printout("Working directory == " + pwd()); + printout("Command-line args == "); + ilcflags.verbosity = Silent; + for (int i = 0; i < argc; ++i) + cout << argv[i] << ' '; + cout << endl; + cout << "sigma == " << sigma[0] << endl; + cout << "T == " << T[0] << endl; + cout << "dPdx == " << ilcflags.dPdx << endl; + cout << "dt == " << dt.dt() << endl; + cout << "ILCFlags == " << ilcflags << endl << endl; + + dsi->setPhaseShifts(xphasehack, zphasehack, uUbasehack); + dsi->chooseMuILC(muname); + if (!restart) { + mu[1] = dsi->mu(); + mu[0] = mu[1] + cflags.initialParamStep; + mu[2] = mu[1] - cflags.initialParamStep; + } + if (muname == "Aspect" || muname == "Diag") { + cout << "aspect ratio || diagonal continuation : " << endl; + cout << setprecision(15); + cout.setf(std::ios::left); + cout << setw(4) << "n" << setw(20) << muname << setw(20) << "aspect" << setw(20) << "diagonal" << endl; + for (int n = 2; n >= 0; --n) + cout << setw(4) << n << setw(20) << mu[n] << setw(20) << u[n].Lx() / u[n].Lz() << setw(20) + << pythag(u[n].Lx(), u[n].Lz()) << endl; + cout.unsetf(std::ios::left); + } + + if (muname == "Lt") { + cout << "Lx,Lz target continuation : " << endl; + cout << setprecision(15); + cout.setf(std::ios::left); + cout << setw(4) << "n" << setw(20) << muname << setw(20) << "Lx" << setw(20) << "Lz" << endl; + for (int i = 2; i >= 0; --i) + cout << setw(4) << i << setw(20) << mu[i] << setw(20) << u[i].Lx() << setw(20) << u[i].Lz() << endl; + cout.unsetf(std::ios::left); + } + + cfarray x(3); + for (int i = 0; i <= 2; ++i) { + dsi->updateMu(mu[i]); + dsi->makeVectorILC(u[i], temp[i], sigma[i], T[i], x[i]); + } + + int Nunk = x[0].rows(); + int Nunk_total = Nunk; +#ifdef HAVE_MPI + MPI_Allreduce(&Nunk, &Nunk_total, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#endif + cout << Nunk_total << " unknowns" << endl; + + Real muFinal = continuation(*dsi, *N, x, mu, cflags); + cout << "Final mu is " << muFinal << endl; + } + cfMPI_Finalize(); + + return 0; +} diff --git a/modules/ilc/programs/ilc_edgetracking.cpp b/modules/ilc/programs/ilc_edgetracking.cpp new file mode 100644 index 00000000..976d7a04 --- /dev/null +++ b/modules/ilc/programs/ilc_edgetracking.cpp @@ -0,0 +1,994 @@ +/** + * Edgetracking program, like the program for pure shear flows. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include +#include +#include "channelflow/dns.h" +#include "channelflow/flowfield.h" +#include "channelflow/poissonsolver.h" +#include "channelflow/utilfuncs.h" +#include "modules/ilc/ilcdsi.h" + +#include + +using namespace std; +using namespace chflow; + +enum Attractor { NoAttractor, Laminar, Turbulent }; + +const string tab = "\t"; + +template +void push(vector& a, T b) { + int l = a.size(); + for (int i = 0; i < l - 1; ++i) { + a[i] = a[i + 1]; + } + a[l - 1] = b; +} + +class AttractorFlags { + public: + AttractorFlags(Real EcfH = 0, Real EcfL = 0, Real L2H = 0, Real L2L = 0, + // Real EnstH = 0, + // Real EnstL = 0, + bool normToHeight = 0, int tMaxDecay = 0); + + // Threshold values + Real EcfH; + Real EcfL; + Real L2H; + Real L2L; + // Real EnstH; + // Real EnstL; + bool normToHeight; + int tMaxDecay; +}; + +class EdgetrackingFlags { + public: + EdgetrackingFlags(); + Real tLastWrite; + Real epsilonAdvance; + Real epsilonBisection; + Real lambdaStep; + Real nBisectTrajectories; + Real tMinAttractor; + Real tMaxAttractor; + bool verifyLH; + bool saveMinima; + Real saveInterval; + bool saveUH; + bool saveUL; + bool keepUL; + bool checkConvergence; + bool returnFieldAtMinEcf; + string savedir; + string directoryEnergy; + string directoryBisections; + string directoryFF; + string directoryMinima; + string filenameEnergyL; + string filenameEnergyH; + string filenameTB; +}; + +class ConvergenceData { + public: + ConvergenceData(); + vector EcfMax; + vector EcfMin; + vector tEcfMin; +}; + +int EdgeStateTracking(FlowField& uL, FlowField& uH, FlowField& tempL, FlowField& tempH, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os); + +void CreateFiles(const EdgetrackingFlags& etflags); +Real Bisection(FlowField& uL, FlowField& uH, FlowField& tempL, FlowField& tempH, Real t0, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os = std::cout); // Return: lambdaH +bool AdvanceLH(FlowField& uL0, FlowField& uH0, FlowField& tempL0, FlowField& tempH0, int& t0, + const EdgetrackingFlags& etflags, const TimeStep& dt, const ILCFlags& flags, + const AttractorFlags& aflags, ostream& os, ConvergenceData& cd); +void ChooseNewLH(const FlowField& u0, const FlowField& temp0, Real& lambdaL, Real& lambdaH, const Real t, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os = std::cout); +void SaveForContinuation(const FlowField& uL, const FlowField& uH, const FlowField& tempL, const FlowField& tempH, + const string dirBisections, Real t); + +Attractor CalcAttractor(const FlowField& u0, const FlowField& temp0, FlowField& returnUField, FlowField& returnTField, + const AttractorFlags& aflags, const TimeStep& dt, const ILCFlags& flags, + const string energyfile, const string saveFinalFF = "", const int tMin = 0, + const int tMax = 10000, bool verbose = false, const int t0 = 0, std::ostream& os = std::cout); +Attractor CalcAttractor(const FlowField& u0, const FlowField& temp0, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& flags, const string energyfile, const string saveFinalFF = "", + const int tMin = 0, const int tMax = 10000, bool verbose = false, const int t0 = 0, + std::ostream& os = std::cout); +Attractor attractor(const FlowField& u, const AttractorFlags& aflags); // attractor is only defined on velocity field +Attractor attractor(const FlowField& u, int& tDecay, Real lastEcf, const AttractorFlags& aflags); +std::string a2s(Attractor a); + +/** Relative L2-Difference between two flowfields */ +Real RelL2Diff(const FlowField& u1, const FlowField& u2) { return L2Dist(u1, u2) / (L2Norm(u1) + L2Norm(u2)); } + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + WriteProcessInfo(argc, argv); + ArgList args(argc, argv, "Find edge state in ILC"); + + ILCFlags ilcflags(args); + TimeStep dt(ilcflags); + ilcflags.verbosity = Silent; + + args.section("Program options"); + const bool cont = args.getflag("-c", "--continue", "continue from a previous calculation"); + + AttractorFlags aflags; + aflags.EcfH = + args.getreal("-EcfH", "--EcfHigh", 0, "threshold in cross flow energy to consider flowfield turbulent"); + aflags.EcfL = + args.getreal("-EcfL", "--EcfLow", 0, "threshold in cross flow energy to consider flowfield laminar"); + aflags.L2H = args.getreal("-L2H", "--L2High", 0, "threshold in L2Norm to consider flowfield turbulent"); + aflags.L2L = args.getreal("-L2L", "--L2Low", 0, "threshold in L2Norm to consider flowfield laminar"); + aflags.normToHeight = + args.getflag("-n2h", "--normToHeight", "norm threshold values to height of box (instead of 1)"); + aflags.tMaxDecay = + args.getint("-maxDcy", "--tMaxDecay", 0, + "Consider a flowfield as laminar if it monotonically decays for that many timeuints"); + + EdgetrackingFlags etflags; + etflags.tMaxAttractor = args.getreal("-tMaxAtt", "--tMaxAttractor", 10000, + "maximum time for waiting that flowfield reaches attractor"); + etflags.tMinAttractor = + args.getreal("-tMinAtt", "--tMinAttractor", 0, "minimum time for waiting that flowfield reaches attractor"); + etflags.epsilonAdvance = args.getreal("-epsA", "--epsilonAdvance", 1e-4, + "maximal separation of flowfields before starting new bisectioning"); + etflags.epsilonBisection = + args.getreal("-epsB", "--epsilonBisection", 1e-6, "maximal separation of flowfields after bisectioning"); + etflags.nBisectTrajectories = + args.getint("-nBis", "--nBisecTrajectories", 1, "Minimium number of trajectories for each bisection"); + const bool chooseLHfirst = args.getflag( + "-cLH", "--chooseLHfirst", + "start edgetracking by choosing uH and uL, i.e. if you are not sure to have a turbulent field"); + etflags.lambdaStep = + args.getreal("-lS", "--lambdaStep", 0.1, "adjust lambda by this value when looking for new laminar state"); + Real lambdaH = args.getreal("-lH", "--lambdaH", 1, "initial value of lambdaH (e.g. for continuation)"); + Real lambdaL = args.getreal("-lL", "--lambdaL", 0, "initial value of lambdaL (e.g. for continuation)"); + etflags.keepUL = args.getflag("-keepUL", "--keepUL", "keep UL instead of UH"); + etflags.verifyLH = + args.getflag("-vrfyLH", "--verifyLH", + "keep on iterating uH and uL after they are separated to verify they are laminar/turbulent"); + etflags.saveInterval = args.getint("-s", "--saveInterval", 0, "save flowfield(s) to disk every s time units"); + etflags.saveUL = args.getbool("-saveUL", "--saveUL", false, "save uL"); + etflags.saveUH = args.getbool("-saveUH", "--saveUH", true, "save uH"); + etflags.saveMinima = args.getbool("-saveMin", "--saveMinima", false, "save at minima of cross flow energy"); + const int nproc0 = + args.getint("-np0", "--nproc0", 0, "number of MPI-processes for transpose/number of parallel ffts"); + const int nproc1 = args.getint("-np1", "--nproc1", 0, "number of MPI-processes for one fft"); + const string logfile = args.getstr("-log", "--logfile", "stdout", "output log (filename or \"stdout\")"); + etflags.savedir = args.getstr("-sd", "--savedir", "./", "path to save directory"); + const string u_ifname = (cont) ? "" : args.getstr(2, "initial flowfield", "initial velocity flowfield"); + const string temp_ifname = (cont) ? "" : args.getstr(1, "initial flowfield", "initial temperature flowfield"); + etflags.directoryEnergy = etflags.savedir + "energy/"; + etflags.directoryBisections = etflags.savedir + "bisections/"; + etflags.directoryFF = etflags.savedir + "flowfields/"; + etflags.directoryMinima = etflags.savedir + "minima/"; + etflags.filenameEnergyL = etflags.savedir + "energyL"; + etflags.filenameEnergyH = etflags.savedir + "energyH"; + etflags.filenameTB = etflags.savedir + "tBisections"; + + args.check(); + args.save(); + ilcflags.save(etflags.savedir); + + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + etflags.tLastWrite = 0; + + ofstream logstream; + if (logfile == "stdout" || logfile == "cout") { + ilcflags.logstream = &cout; + } else { + if (cfmpi->taskid() == 0) + logstream.open((etflags.savedir + logfile).c_str()); + ilcflags.logstream = &logstream; + } + ostream& os = (ostream&)*ilcflags.logstream; + + printout("Edge state tracking algorithm", os); + printout("Working directory == " + pwd(), os); + printout("Command-line args == ", false, os); + for (int i = 0; i < argc; ++i) + printout(string(argv[i]) + " ", false, os); + stringstream sstr; + sstr << ilcflags; + printout("\nILCFlags: " + sstr.str(), os); + + if ((ilcflags.symmetries.length() > 0) || (ilcflags.tempsymmetries.length() > 0)) { + printout("Restricting flow to invariant subspace generated by symmetries", os); + stringstream sstr2u, sstr2t; + sstr2u << "Velocity: " << ilcflags.symmetries; + sstr2t << "Temperature: " << ilcflags.tempsymmetries; + printout(sstr2u.str(), os); + printout(sstr2t.str(), os); + } + + printout("", os); // newline + + // ******** Preparations ******** + Real t = 0; + FlowField u, uL, uH; + FlowField temp, tempL, tempH; + + fftw_loadwisdom(); + if (!cont) { + CreateFiles(etflags); + u = FlowField(u_ifname, cfmpi); + temp = FlowField(temp_ifname, cfmpi); + if (temp.Nd() != 1) + cferror("Error initializing temp field: number of dimensions must be 1!"); + uH = u; // temp for calculating projection error and optimizing FFTW + if (ilcflags.symmetries.length() > 0) { + printout("Projecting Velocity FlowField to invariant subspace generated by symmetries", os); + stringstream sstr2; + sstr2 << ilcflags.symmetries; + printout(sstr2.str(), os); + u.project(ilcflags.symmetries); + printout("Projection error == " + r2s(L2Dist(u, uH))); + } + tempH = temp; // tmp for calculating projection error and optimizing FFTW + if (ilcflags.tempsymmetries.length() > 0) { + printout("Projecting Temperature FlowField to invariant subspace generated by symmetries", os); + stringstream sstr2; + sstr2 << ilcflags.tempsymmetries; + printout(sstr2.str(), os); + temp.project(ilcflags.tempsymmetries); + printout("Projection error == " + r2s(L2Dist(temp, tempH))); + } + uH.optimizeFFTW(FFTW_PATIENT); + tempH.optimizeFFTW(FFTW_PATIENT); + + fftw_savewisdom(); + if (chooseLHfirst) { + ChooseNewLH(u, temp, lambdaL, lambdaH, 0., etflags, aflags, dt, ilcflags, os); + } + uH = u; + uH *= lambdaH; + uL = u; + uL *= lambdaL; + tempH = temp; + tempH *= lambdaH; + tempL = temp; + tempL *= lambdaL; + } else { + ifstream Tlast(etflags.filenameTB); + while (Tlast >> t) + ; + ilcflags.t0 = t; + string fffileH = etflags.directoryBisections + "uH" + r2s(t); + string fffileL = etflags.directoryBisections + "uL" + r2s(t); + string tempfileH = etflags.directoryBisections + "tempH" + r2s(t); + string tempfileL = etflags.directoryBisections + "tempL" + r2s(t); + printout("Loading flowfields from " + fffileL + " and " + fffileH, os); + uH = FlowField(fffileH, cfmpi); + uH.project(ilcflags.symmetries); + uL = FlowField(fffileL, cfmpi); + uL.project(ilcflags.symmetries); + tempH = FlowField(tempfileH, cfmpi); + tempH.project(ilcflags.tempsymmetries); + tempL = FlowField(tempfileL, cfmpi); + tempL.project(ilcflags.tempsymmetries); + u = uH; + uH.optimizeFFTW(FFTW_PATIENT); + uH = u; + + // Read file energyH and try to extract last time when something was saved (it's the first number in the + // last not-empty line not beginning with #) + ifstream fin(etflags.filenameEnergyH.c_str()); + if (fin) { + string s1; + istringstream s2; + while (getline(fin, s1)) { + if (s1 != "" && s1.find("#") == string::npos) + s2.str(s1); + } + s2 >> etflags.tLastWrite; + s1 = "Found " + i2s((int)etflags.tLastWrite) + " as last time something was written to energyH"; + printout(s1, os); + } + if (etflags.tLastWrite >= ilcflags.T - 1) { + cferror("Already over time limit"); + } + } + + EdgeStateTracking(uL, uH, tempL, tempH, etflags, aflags, dt, ilcflags, os); + } + cfMPI_Finalize(); + return 0; +} + +/** +The actual edge state tracking algorithm +*/ +int EdgeStateTracking(FlowField& uL_, FlowField& uH_, FlowField& tempL_, FlowField& tempH_, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os) { + int t = ilcflags.t0; + FlowField u, uL(uL_), uH(uH_); + FlowField temp, tempL(tempL_), tempH(tempH_); + + printout("Starting edge tracking at t=" + i2s((int)t), os); + + // Main timestepping loop + bool finished = false; + ConvergenceData convergencedata; + while (!finished && t <= ilcflags.T) { + // Step 1: Bisectioning + printout("\nt = " + i2s((int)t) + ": Starting Bisection"); + // Write time to file + if (uL_.taskid() == 0) { + ofstream fT((etflags.filenameTB).c_str(), ios::app); + fT << t << endl; + fT.close(); + } + + Bisection(uL, uH, tempL, tempH, t, etflags, aflags, dt, ilcflags, os); + + // Step 2: Advance FlowFields + printout("t = " + i2s((int)t) + ": Advancing flowfields"); + bool finished = AdvanceLH(uL, uH, tempL, tempH, t, etflags, dt, ilcflags, aflags, os, convergencedata); + if (!finished) { + if (etflags.keepUL) { + u = uL; + temp = tempL; + } else { + u = uH; + temp = tempH; + } + + // Step 3: Choose new Flowfields + Real lambdaL = 1, lambdaH = 1; + ChooseNewLH(u, temp, lambdaL, lambdaH, t, etflags, aflags, dt, ilcflags, os); + uL = u; + uL *= lambdaL; + uH = u; + uH *= lambdaH; + tempL = temp; + tempL *= lambdaL; + tempH = temp; + tempH *= lambdaH; + } + } + + return 0; +} + +/** Create files for energy and bisection times, if not continuing. + */ +void CreateFiles(const EdgetrackingFlags& etflags) { + // Create files + int taskid = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &taskid); +#endif + if (taskid == 0) { + ios::openmode mode = ios::out; + mkdir(etflags.directoryEnergy); + mkdir(etflags.directoryBisections); + if (etflags.saveInterval != 0) + mkdir(etflags.directoryFF); + if (etflags.saveMinima) + mkdir(etflags.directoryMinima); + ofstream fL((etflags.filenameEnergyL).c_str(), mode); + ofstream fH((etflags.filenameEnergyH).c_str(), mode); + fL << ilcfieldstatsheader_t("t") << endl; + fH << ilcfieldstatsheader_t("t") << endl; + fL.close(); + fH.close(); + + // Bisection times + ofstream fT((etflags.filenameTB).c_str(), mode); + fT.close(); + } +} + +/** Find the edge of chaos between uL and uH. + * @param uL lower flowfield, known to decay + * @param uH upper flowfield, known to become turbulent + * @param t0 time, only for filenames + */ +Real Bisection(FlowField& uL, FlowField& uH, FlowField& tempL, FlowField& tempH, Real t0, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os) { + Real lambdaL = 0, lambdaH = 1; + Real l2d = RelL2Diff(uL, uH); + int i = 0; + FlowField u0(uH); + + printout( + "Initial L2Diff ( =L2(uH-uL)/(L2(uL)+L2(uH)) ) is " + r2s(l2d) + " and eps is " + r2s(etflags.epsilonBisection), + os); + int nTraj = etflags.nBisectTrajectories; + + SaveForContinuation(uL, uH, tempL, tempH, etflags.directoryBisections, t0); + while (l2d > etflags.epsilonBisection) { + i++; + if (nTraj == 1) { + Real lambda = (lambdaL + lambdaH) / 2.; + FlowField uT = uH; + uT += uL; + uT *= 0.5; + FlowField tempT = tempH; + tempT += tempL; + tempT *= 0.5; + + // Calculate attractor of ulambda and update lambaL or lambdaH accordingly + string f = ""; + Attractor a = CalcAttractor(uT, tempT, aflags, dt, ilcflags, + etflags.directoryEnergy + "e_bisection_t" + i2s(int(t0)) + "_i" + i2s(i), "", + etflags.tMinAttractor, etflags.tMaxAttractor, false, t0); + if (a == Turbulent) { + lambdaH = lambda; + uH = uT; + tempH = tempT; + f = "turbulent"; + } else if (a == Laminar) { + lambdaL = lambda; + uL = uT; + tempL = tempT; + f = "laminar"; + } else { + cferror("Bisection did not converge"); + } + l2d = RelL2Diff(uL, uH); + printout( + "Bisection t=" + i2s(t0) + ", i=" + FillSpaces(i, 3) + ": Found " + f + ", new L2Diff is " + r2s(l2d), + os); + } else { + // Take several steps for each "bisection" + // i.e. u = (1-alpha) uL + alpha uH + // abort, once u is turbulent and set uL to the alpha-value _two_ steps lower + Real alphaStep = 1. / (nTraj + 1); + + Attractor a = Laminar; + int j = 0; + Real alpha = 0; + FlowField uT, tempT; + while (a == Laminar && j <= nTraj + 10) { + j++; + // uT = (1-alpha)*uL + alpha*uH + alpha = j * alphaStep; + uT = uL; + uT *= (1. - alpha) / alpha; + uT += uH; + uT *= alpha; + tempT = tempL; + tempT *= (1. - alpha) / alpha; + tempT += tempH; + tempT *= alpha; + a = CalcAttractor( + uT, tempT, aflags, dt, ilcflags, + etflags.directoryEnergy + "e_bisection_t" + i2s(int(t0)) + "_i" + i2s(i) + "_j" + i2s(j), "", + etflags.tMinAttractor, etflags.tMaxAttractor, false, t0); + string f = "laminar"; + if (a == Turbulent) + f = "turbulent"; + printout("Bisection t=" + i2s(t0) + ", i=" + FillSpaces(i, 3) + ", j=" + FillSpaces(j, 3) + + ", alpha=" + r2s(alpha) + ": Found " + f, + os); + } + if (a == Laminar) + cferror("Did not find turbulent state, aborting"); + + if (j >= 3) { + // uL = (1-alpha[j-2])*uL + alpha[j-2]*uH + Real alphaL = (j - 2) * alphaStep; + uL *= (1. - alphaL) / alphaL; + uL += uH; + uL *= alphaL; + tempL *= (1. - alphaL) / alphaL; + tempL += tempH; + tempL *= alphaL; + } + + uH = uT; // because uT is the FlowField we know to become turbulent + tempH = tempT; + + l2d = RelL2Diff(uL, uH); + printout("Bisection t=" + i2s(t0) + ", i=" + FillSpaces(i, 3) + ": lL=" + r2s(lambdaL) + + ", lH=" + r2s(lambdaH) + ", new L2diff is " + r2s(l2d) + "\n", + os); + } + SaveForContinuation(uL, uH, tempL, tempH, etflags.directoryBisections, t0); + } + + printout("Bisection converged, L2Diff = " + r2s(l2d), os); + printout(" lambdaL = " + r2s(lambdaL) + ", lambdaH = " + r2s(lambdaH), os); + printout( + " l2(uH)/l2(u0) = " + r2s(L2Norm(uH) / L2Norm(u0)) + ", l2(uL)/l2(u0) = " + r2s(L2Norm(uL) / L2Norm(u0)), + os); + return lambdaH; +} + +/** Advance uL and uH. + * Advance till L2Diff(uL,uH)/L2Norm(uH)>c_epsilon or t > c_T. Check if they ar really turbulent/laminar if verifyLH is + * active + * + * @param uL0 lower flowfield, iterated while not stateIsOnEdge + * @param uH0 higher flowfield, iterated while not stateIsOnEdge + * @param t time when advancing starts + * + * @return t, final time + */ +bool AdvanceLH(FlowField& uL0, FlowField& uH0, FlowField& tempL0, FlowField& tempH0, int& t, + const EdgetrackingFlags& etflags, const TimeStep& dt, const ILCFlags& flags, + const AttractorFlags& aflags, ostream& os, ConvergenceData& cd) { + // Create TimeStep, pressureField and DNS objects + // Real t = t0; // Actual time + int t0 = t; + + TimeStep dtL = dt; + TimeStep dtH = dt; + vector fieldsH(3); + vector fieldsL(3); + fieldsL[0] = uL0; // is uL + fieldsH[0] = uH0; // is uH + fieldsL[1] = tempL0; // is tempL + fieldsH[1] = tempH0; // is tempH + fieldsL[2] = + FlowField(uH0.Nx(), uH0.Ny(), uH0.Nz(), 1, uH0.Lx(), uH0.Lz(), uH0.a(), uH0.b(), fieldsL[0].cfmpi()); // is qL + fieldsH[2] = + FlowField(uH0.Nx(), uH0.Ny(), uH0.Nz(), 1, uH0.Lx(), uH0.Lz(), uH0.a(), uH0.b(), fieldsH[0].cfmpi()); // is qH + ILC dnsH(fieldsH, flags); + ILC dnsL(fieldsL, flags); + ChebyCoeff ubase(laminarProfile(flags, uH0.a(), uH0.b(), uH0.Ny())); + ChebyCoeff wbase(uH0.Ny(), uH0.a(), uH0.b()); + PressureSolver psolver(uH0, ubase, wbase, flags.nu, flags.Vsuck, flags.nonlinearity); + + psolver.solve(fieldsH[2], fieldsH[0]); + psolver.solve(fieldsL[2], fieldsL[0]); + + // Open Files + ofstream fL, fH; + if (uH0.taskid() == 0) { + fL.open((etflags.filenameEnergyL).c_str(), ios::app); + fH.open((etflags.filenameEnergyH).c_str(), ios::app); + } + + // Advance uH and uL + Real t2 = t; // Internal time while verifying + Real tAttractor = etflags.tMaxAttractor + t; // Maximum time to verify + bool turbulent = false, laminar = false, stateIsOnEdge = true; + Real ecf0 = 0, ecf1 = 0, ecf2 = 0; // ecf(uH, t - 0,1,2) + vector ecf; + for (int i = 0; i < 100; ++i) + ecf.push_back(0); + while ((!turbulent || !laminar) && t <= flags.T && (stateIsOnEdge || (etflags.verifyLH && t2 <= tAttractor))) { + t2++; + // Advance uH + if (!turbulent) { + dnsH.advance(fieldsH, dtH.n()); + Real cflH = dnsH.CFL(fieldsH[0]); + if (dtH.adjust(cflH, false)) { + dnsH.reset_dt(dtH); + } + } + + // Advance uL + if (!laminar) { + dnsL.advance(fieldsL, dtL.n()); + Real cflL = dnsL.CFL(fieldsL[0]); + if (dtL.adjust(cflL, false)) { + dnsL.reset_dt(dtL); + } + } + + if (stateIsOnEdge) { + Real l2d = RelL2Diff(fieldsL[0], fieldsH[0]); + if (l2d < etflags.epsilonAdvance) { + t++; + if (t > etflags.tLastWrite) { // prevent writing the same stuff twice + printout(ilcfieldstats_t(fieldsH[0], fieldsH[1], t), fH); + printout(ilcfieldstats_t(fieldsL[0], fieldsL[1], t), fL); + } + ecf2 = ecf1; + ecf1 = ecf0; + ecf0 = Ecf(fieldsH[0]); + bool ecfMinimum = (ecf1 < ecf2) && (ecf1 <= ecf0); + bool ecfMaximum = (ecf1 > ecf2) && (ecf1 >= ecf0); + + // Check whether we are converged to a (relative) equilibrium where ecf doesn't change anymore + if (etflags.checkConvergence) { + bool eq = true; + int l = ecf.size(); + ecf[t % l] = ecf0; + for (int i = 0; i < l; ++i) { + if (ecf[i] < 1e-16 || abs(ecf[i] - ecf[0]) / ecf[0] > etflags.epsilonBisection) { + eq = false; + break; + } + } + if (eq) { + printout("Ecf is constant -- converged to a (relative) equilibrium", os); + return true; + } + + // Check whether we are converged to a periodic orbit (somewhat more subtle) + // We consider a state as periodic if either 4 minima and 4 maxima in a row occur with same distance + // and same values or when maxima,minima 1,3 and 2,4 are the same and occur at the same distance + if (ecfMaximum) + push(cd.EcfMax, ecf0); + if (ecfMinimum) { + push(cd.EcfMin, ecf0); + push(cd.tEcfMin, t); + + Real eps = 1e-3; + if (abs(cd.EcfMax[0] - cd.EcfMax[2]) / ecf0 < eps && + abs(cd.EcfMax[1] - cd.EcfMax[3]) / ecf0 < eps && + abs(cd.EcfMin[0] - cd.EcfMin[2]) / ecf0 < eps && + abs(cd.EcfMin[1] - cd.EcfMin[3]) / ecf0 < eps && + abs((cd.tEcfMin[3] - cd.tEcfMin[1]) - (cd.tEcfMin[2] - cd.tEcfMin[0])) < 5) { + printout("Converged to a periodic orbit", os); + return true; + } + } + } + + // Save fields + bool saveIntervalReached = etflags.saveInterval != 0 && t % (int)etflags.saveInterval == 0; + if (saveIntervalReached) { + if (etflags.saveUH) { + fieldsH[0].save(etflags.directoryFF + "uH" + i2s((int)t)); + fieldsH[1].save(etflags.directoryFF + "tH" + i2s((int)t)); + } + if (etflags.saveUL) { + fieldsL[0].save(etflags.directoryFF + "uL" + i2s((int)t)); + fieldsL[1].save(etflags.directoryFF + "tL" + i2s((int)t)); + } + } + if (ecfMinimum && etflags.saveMinima) { + if (etflags.saveUH) { + fieldsH[0].save(etflags.directoryMinima + "uH" + i2s((int)t)); + fieldsH[1].save(etflags.directoryMinima + "tH" + i2s((int)t)); + } + if (etflags.saveUL) { + fieldsL[0].save(etflags.directoryMinima + "uL" + i2s((int)t)); + fieldsL[1].save(etflags.directoryMinima + "tL" + i2s((int)t)); + } + } + uL0 = fieldsL[0]; + uH0 = fieldsH[0]; + tempL0 = fieldsL[1]; + tempH0 = fieldsH[1]; + } else { + // The flowfields are separated too far + // Use state when this occurs for the first time for the next bisectioning + // and keep on iterating to verify they become turbulent/laminar + printout("L2(uL) " + r2s(L2Norm(fieldsL[0])), os); + printout("L2(uH) " + r2s(L2Norm(fieldsH[0])), os); + stateIsOnEdge = false; + printout("Flowfields separated by " + r2s(RelL2Diff(fieldsL[0], fieldsH[0])), os); + printout("Advanced flowfields to t = " + i2s((int)t), os); + printout("", os); // newline + } + } + + // Check whether flowfields are laminar/turbulent + if (!turbulent && attractor(fieldsH[0], aflags) == Turbulent) { + turbulent = true; + printout("uH is turbulent at t = " + i2s((int)t2), os); + if (!etflags.verifyLH) { + cferror("uH is turbulent, check thresholds"); + } + if (t2 - t0 <= 2) { + cferror("t2 - t0 = " + r2s(t2 - t0) + ", aborting\nuH is turbulent, check your thresholds!"); + } + } + if (!laminar && attractor(fieldsL[0], aflags) == Laminar) { + laminar = true; + printout("uL is laminar at t = " + i2s((int)t2), os); + if (!etflags.verifyLH) { + cferror("uL is laminar, check thresholds"); + } + if (t2 - t0 <= 2) { + cferror("t2 - t0 = " + r2s(t2 - t0) + ", aborting\nuL is laminar, check your thresholds!"); + } + } + } + + // Abort if the flowfields did not reach the expected state + if (!laminar && etflags.verifyLH) { + cferror("verifyLH: uL does not become laminar"); + } + if (!turbulent && etflags.verifyLH) { + cferror("verifyLH: uH does not become turbulent"); + } + + // Close files + if (fieldsH[0].taskid() == 0) { + fH.close(); + fL.close(); + } + + if (t == flags.T) + return true; + return false; +} + +/** Choose two flowfields such that uL becomes laminar and uH turbulent. + * @param u0 flowfield to scale up/down + * @param lambdaL return reference for lambdaL + * @param lambdaH return reference for lambdaH + */ +void ChooseNewLH(const FlowField& u0, const FlowField& temp0, Real& lambdaL, Real& lambdaH, const Real t, + const EdgetrackingFlags& etflags, const AttractorFlags& aflags, const TimeStep& dt, + const ILCFlags& ilcflags, ostream& os) { + assert(lambdaL <= lambdaH); + printout("Choosing new laminar/turbulent flowfields", os); + bool equalLambda = (lambdaH == lambdaL) ? true : false; + + FlowField u = u0; + FlowField temp = temp0; + Attractor AH = NoAttractor, AL = NoAttractor; + // Choose turbulent field + while (AH != Turbulent) { + u = u0; + u *= lambdaH; + temp = temp0; + temp *= lambdaH; + AH = CalcAttractor(u, temp, aflags, dt, ilcflags, + etflags.directoryEnergy + "e_chooseLH_t" + i2s(int(t)) + "_lambdaH" + r2s(lambdaH), "", + etflags.tMinAttractor, etflags.tMaxAttractor, false, t); + + printout("FlowField with lambda = " + r2s(lambdaH) + " is " + a2s(AH), os); + if (AH == Laminar) { + lambdaL = lambdaH; + lambdaH += etflags.lambdaStep; + AL = Laminar; + } else if (AH == NoAttractor) { + cferror("flowfield is ambigous"); + } + } + // Choose laminar field + while (AL != Laminar) { + if (equalLambda) { + lambdaL -= etflags.lambdaStep; + equalLambda = false; + } + u = u0; + u *= lambdaL; + temp = temp0; + temp *= lambdaL; + // AL = Attractor ( u, "_chooseLH_t" + i2s(int(t)) + "_lambdaL" + r2s(lambdaL) ); + AL = CalcAttractor(u, temp, aflags, dt, ilcflags, + etflags.directoryEnergy + "e_chooseLH_t" + i2s(int(t)) + +"_lambdaL" + r2s(lambdaL), "", + etflags.tMinAttractor, etflags.tMaxAttractor, false, t); + + printout("FlowField with lambda = " + r2s(lambdaL) + " is " + a2s(AL), os); + if (AL == Turbulent) { + lambdaH = lambdaL; + lambdaL -= etflags.lambdaStep; + } else if (AL == NoAttractor) { + cferror("flowfield is ambigous"); + } + } + printout("Finished choosing uH, uL -- lambdaH = " + r2s(lambdaH) + ", lambdaL = " + r2s(lambdaL), os); + FlowField uL(u0); + uL *= lambdaL; + FlowField uH(u0); + uH *= lambdaH; + FlowField tempL(temp0); + tempL *= lambdaL; + FlowField tempH(temp0); + tempH *= lambdaH; + SaveForContinuation(uL, uH, tempL, tempH, etflags.directoryBisections, t); +} + +/** Save all data needed for continuation to file 'continue' + */ +void SaveForContinuation(const FlowField& uL, const FlowField& uH, const FlowField& tempL, const FlowField& tempH, + const string directoryBisections, Real t) { + string ffnameH = directoryBisections + "uH" + i2s((int)t); + string ffnameL = directoryBisections + "uL" + i2s((int)t); + uH.save(ffnameH); + uL.save(ffnameL); + + string tffnameH = directoryBisections + "tH" + i2s((int)t); + string tffnameL = directoryBisections + "tL" + i2s((int)t); + tempH.save(tffnameH); + tempL.save(tffnameL); + + fftw_savewisdom(); +} + +EdgetrackingFlags::EdgetrackingFlags() + : tLastWrite(0), + epsilonAdvance(1e-4), + epsilonBisection(1e-6), + lambdaStep(0.1), + nBisectTrajectories(1), + tMinAttractor(0), + tMaxAttractor(10000), + verifyLH(false), + saveMinima(false), + saveInterval(0), + saveUH(true), + saveUL(false), + keepUL(false), + checkConvergence(false), + returnFieldAtMinEcf(false), + savedir("./"), + directoryEnergy("energy/"), + directoryBisections("bisections/"), + directoryFF("flowfields/"), + directoryMinima("minima/"), + filenameEnergyL("energyL"), + filenameEnergyH("energyH"), + filenameTB("tBisections") {} + +ConvergenceData::ConvergenceData() { + for (int i = 0; i < 4; ++i) { + EcfMin.push_back(0); + EcfMax.push_back(0); + tEcfMin.push_back(0); + } +} + +/** Calculate the attractor of a given flowfield. + * The decicision whether a flowfield is laminar or turbulent is made based upon the parameters in the configfile. + */ +Attractor CalcAttractor(const FlowField& u0, const FlowField& temp0, const AttractorFlags& aflags, const TimeStep& dt_, + const ILCFlags& flags, const string energyfile, const string saveFinalFF, const int tMin, + const int tMax, bool verbose, const int t0, ostream& os) { + FlowField u(u0); + FlowField temp(temp0); + return CalcAttractor(u0, temp0, u, temp, aflags, dt_, flags, energyfile, saveFinalFF, tMin, tMax, verbose, t0, os); +} + +Attractor CalcAttractor(const FlowField& u0, const FlowField& temp0, FlowField& u, FlowField& temp, + const AttractorFlags& aflags, const TimeStep& dt_, const ILCFlags& flags, + const string energyfile, const string saveFinalFF, const int tMin, const int tMax, bool verbose, + const int t0, ostream& os) { + u = u0; + temp = temp0; + // Check what data to save and create files + + // Advance flowfield till it reaches an attractor or time limit + TimeStep dt = dt_; + FlowField q(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b(), u.cfmpi()); + vector fields = {u, temp, q}; + ILC dns(fields, flags); + ChebyCoeff ubase = laminarProfile(flags, u.a(), u.b(), u.Ny()); + ChebyCoeff wbase = ChebyCoeff(u.Ny(), u.a(), u.b()); + PressureSolver psolver(u, ubase, wbase, flags.nu, flags.Vsuck, flags.nonlinearity); + psolver.solve(fields[2], fields[0]); + ofstream f; + bool saveEnergy = (energyfile != ""); + if (saveEnergy) { + if (u0.taskid() == 0) + f.open(energyfile.c_str(), ios::out); + printout(ilcfieldstatsheader_t("t"), f); + printout(ilcfieldstats_t(u, temp, t0), f); + } + Attractor result = NoAttractor; + int t = 0; + int tDecay = 0; // timeunits that flowfield is monotonically decaying + while (result == NoAttractor && t < tMax) { + t++; + Real lastEcf = Ecf(fields[0]); // for calculating decay time + dns.advance(fields, dt.n()); + + Real cfl = dns.CFL(fields[0]); + if (dt.adjust(cfl, verbose, os)) { + dns.reset_dt(dt); + if (verbose) + printout("Adjusting dt to " + r2s(dt.dt()) + ", CFL is " + r2s(cfl), os); + } + + if (saveEnergy) { + printout(ilcfieldstats_t(fields[0], fields[1], t + t0), f); + } + + // Check if FF is turbulent or laminar + if (t > tMin) { + result = attractor(fields[0], tDecay, lastEcf, aflags); + } + } + // printout("", os); + if (saveFinalFF != "") + fields[0].save(saveFinalFF.c_str()); + if (saveEnergy && fields[0].taskid() == 0) + f.close(); + return result; +} + +AttractorFlags::AttractorFlags(Real EcfH_, Real EcfL_, Real L2H_, Real L2L_, bool n2h, int tMaxDecay_) + : EcfH(EcfH_), + EcfL(EcfL_), + L2H(L2H_), + L2L(L2L_), + // EnstH( EnstH_ ), + // EnstL( EnstL_ ), + normToHeight(n2h), + tMaxDecay(tMaxDecay_) {} + +Attractor attractor(const FlowField& u, const AttractorFlags& a) { + int t = 0; + return attractor(u, t, 0., a); +} + +Attractor attractor(const FlowField& u, int& tDecay, Real lastEcf, const AttractorFlags& a) { + const bool useEcf = (a.EcfH > 1e-12 && a.EcfL > 1e-12); + // const bool useEnst = ( a.EnstH > 1e-12 && a.EnstL > 1e-12 ); + const bool useL2 = (a.L2H > 1e-12 && a.L2L > 1e-12); + + if (not(useEcf || useL2)) { + cferror("No valid attractor settings specified\n"); + } + + // A FlowField has reached an attractor if all conditions are true + bool isLaminar = true, isTurbulent = true; + Real val = 0; + if (useEcf) { + val = Ecf(u); + if (a.normToHeight) + val /= (u.b() - u.a()); + if (val < a.EcfH) + isTurbulent = false; + if (val > a.EcfL) + isLaminar = false; + } + if (useL2) { + val = L2Norm(u); + if (a.normToHeight) + val /= sqrt((u.b() - u.a())); + if (val < a.L2H) + isTurbulent = false; + if (val > a.L2L) + isLaminar = false; + } + // if ( useEnst ) { + // val = Enstrophy(u); + // if ( a.normToHeight ) val /= sqrt((u.b() - u.a())); + // if ( val < a.EnstH ) + // isTurbulent = false; + // if ( val > a.EnstL ) + // isLaminar = false; + // } + + if (a.tMaxDecay > 0) { + if (Ecf(u) <= lastEcf) + tDecay++; + else + tDecay = 0; + + if (tDecay > a.tMaxDecay) + isLaminar = true; + } + + if (isTurbulent && isLaminar) { + cferror("FlowField is both laminar and turbulent, exiting"); + } else if (isTurbulent) { + return Turbulent; + } else if (isLaminar) { + return Laminar; + } + return NoAttractor; +} + +string a2s(Attractor a) { + if (a == Turbulent) + return "turbulent"; + if (a == Laminar) + return "laminar"; + return "neither laminar nor turbulent"; +} diff --git a/modules/ilc/programs/ilc_findeigenvals.cpp b/modules/ilc/programs/ilc_findeigenvals.cpp new file mode 100644 index 00000000..d838b630 --- /dev/null +++ b/modules/ilc/programs/ilc_findeigenvals.cpp @@ -0,0 +1,273 @@ +/** + * Eigenvalue program, like the program for pure shear flows. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include +#include +#include +#include +#include +#include "cfbasics/cfvector.h" +#include "channelflow/chebyshev.h" +#include "channelflow/flowfield.h" +#include "channelflow/symmetry.h" +#include "channelflow/tausolver.h" +#include "modules/ilc/ilcdsi.h" +#include "nsolver/nsolver.h" + +using namespace std; +using namespace chflow; + +// This program calculates eigenvalues of fixed point of plane Couette flow +// using Arnoldi iteration. The ideas and algorithm are based on Divakar +// Viswanath, "Recurrent motions within plane Couette turbulence", +// J. Fluid Mech. 580 (2007), http://arxiv.org/abs/physics/0604062. + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + WriteProcessInfo(argc, argv); + int taskid = 0; +#ifdef HAVE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &taskid); +#endif + + string purpose( + "compute spectrum of eigenvalues of equilibria, traveling waves, or periodic orbit using Arnoldi " + "iteration"); + ArgList args(argc, argv, purpose); + + ILCFlags ilcflags(args); + TimeStep dt(ilcflags); + ilcflags.verbosity = Silent; + + // The Eigenvals class is utilized to solve the eigenvalue problem. + // This class requires Arnoldi class. + unique_ptr E; + EigenvalsFlags eigenflags(args); + E = unique_ptr(new Eigenvals(eigenflags)); + + args.section("Program options"); + const bool poincare = + args.getflag("-poinc", "--poincare", "computing eigenvalues of map on I-D=0 Poincare section"); + + const string sigstr = + args.getstr("-sigma", "--sigma", "", "file containing sigma of sigma f^T(u) - u = 0 (default == identity)"); + + const int seed = args.getint("-sd", "--seed", 1, "seed for random number generator"); + const Real smooth = args.getreal("-s", "--smoothness", 0.4, "smoothness of initial perturb, 0 < s < 1"); + const Real EPS_du = + args.getreal("-edu", "--epsdu", 1e-7, "magnitude of Arnoldi perturbation to velocity AND temperature"); + const string duname = + args.getstr("-du", "--perturbu", "", "initial perturbation velocity field, random if unset"); + const string dtname = + args.getstr("-dtemp", "--perturbt", "", "initial perturbation temperature field, random if unset"); + + const int nproc0 = + args.getint("-np0", "--nproc0", 0, "number of MPI-processes for transpose/number of parallel ffts"); + const int nproc1 = args.getint("-np1", "--nproc1", 0, "number of MPI-processes for one fft"); + + const string uname = args.getstr(2, "", "filename of EQB, TW, or PO velocity solution"); + const string tname = args.getstr(1, "", "filename of EQB, TW, or PO temperature solution"); + + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + + args.check(); + args.save("./"); + args.save(eigenflags.outdir); + + fftw_loadwisdom(); + + srand48(seed); + const Real decay = 1.0 - smooth; + + PoincareCondition* h = poincare ? new DragDissipation() : 0; + + FlowField u(uname, cfmpi); // u*, the solution of sigma f^T(u*,temp*) - u*,temp* = 0 + FlowField temp(tname, cfmpi); // temp*, the solution of sigma f^T(u*,temp*) - u*,temp* = 0 + + const int Nx = u.Nx(); + const int Ny = u.Ny(); + const int Nz = u.Nz(); + + const int kxmin = -u.kxmaxDealiased(); + const int kxmax = u.kxmaxDealiased(); + const int kzmin = 0; + const int kzmax = u.kzmaxDealiased(); + + if (taskid == 0) { + cout << setprecision(17); + cout << " Nx == " << Nx << endl; + cout << " Ny == " << Ny << endl; + cout << " Nz == " << Nz << endl; + cout << "kxmin == " << kxmin << endl; + cout << "kxmax == " << kxmax << endl; + cout << "kzmin == " << kzmin << endl; + cout << "kzmax == " << kzmax << endl; + + cout << "dt == " << dt.dt() << endl; + cout << "dtmin == " << dt.dtmin() << endl; + cout << "dtmax == " << dt.dtmax() << endl; + cout << "CFLmin == " << dt.CFLmin() << endl; + cout << "CFLmax == " << dt.CFLmax() << endl; + } + if (!poincare) + dt.adjust_for_T(ilcflags.T); + + FieldSymmetry sigma; // defaults to identity + if (sigstr.length() != 0) + sigma = FieldSymmetry(sigstr); + + if (ilcflags.symmetries.length() > 0) { + if (taskid == 0) { + cout << "Restricting flow to invariant subspace generated by symmetries" << endl; + cout << "VelocityField: " << ilcflags.symmetries << endl; + cout << "TemperatureField: " << ilcflags.tempsymmetries << endl; + } + } + + if (taskid == 0) + cout << "DNS flags = " << ilcflags << endl; + ilcflags.save(eigenflags.outdir); + + // Set up DNS operator ("A" in Arnoldi A*b terms) + if (taskid == 0) + cout << "setting up DNS and initial fields..." << endl; + ChebyTransform trans(Ny); + + const Real l2u = L2Norm(u); + const Real eps = + (l2u < EPS_du) + ? EPS_du + : EPS_du / l2u; // for choice of epsilon, see eq. (15) in reference + // C.J. Mack, P.J. Schmid/Journal of Computational Physics 229 (2010) 541-560 + FlowField Gu(u); + FlowField Gtemp(temp); + if (taskid == 0) + cout << "computing sigma f^T(u)..." << endl; + + // Construct the dynamical-systems interface object depending on the given parameters. + unique_ptr dsi; + dsi = unique_ptr( + new ilcDSI(ilcflags, sigma, h, dt, false, false, false, false, 0.0, u, temp, E->getLogstream())); + + // Check if sigma f^T(u) - u = 0 + Eigen::VectorXd x; + field2vector(u, temp, x); + + Eigen::VectorXd Gx = dsi->eval(x); + vector2field(Gx, Gu, Gtemp); + + if (taskid == 0) + cout << "\nCFL == " << dsi->getCFL() << endl; + + Real l2normGu = L2Norm(Gu); + if (taskid == 0) + cout << "L2Norm(u - sigma f^T(u)) = " << l2normGu << endl; + + if (l2normGu > 1e-6) + cferror("error: (u, sigma, T) is not a solution such as sigma f^T(u) - u = 0"); + + FlowField du(u); + FlowField dtemp(temp); + du.setToZero(); + dtemp.setToZero(); + // Set du = EPS_du (random unit perturbation, "b" in Arnoldi A*b terms) + if (duname.length() == 0) { + cout << "Constructing du..." << endl; + + bool meanflow_perturb = true; + du.addPerturbations(kxmax, kzmax, 1.0, decay, meanflow_perturb); + + if (ilcflags.constraint == PressureGradient) { + // Modify du so that (du/dy|a + du/dy|b) == (dw/dy|a + dw/dy|b) == 0 + // i.e. mean wall shear == 0 + Real h = (du.b() - du.a()) / 2; + if (du.taskid() == du.task_coeff(0, 0)) { + ChebyCoeff du00 = Re(du.profile(0, 0, 0)); + ChebyCoeff dw00 = Re(du.profile(0, 0, 2)); + ChebyCoeff du00y = diff(du00); + ChebyCoeff dw00y = diff(dw00); + Real duy = (du00y.eval_a() + du00y.eval_b()) / 2; + Real dwy = (dw00y.eval_a() + dw00y.eval_b()) / 2; + + cout << "Modifying du so that it doesn't change mean pressure balance..." << endl; + cout << "pre-mod : " << endl; + cout << "(duya + duyb)/2 == " << duy << endl; + cout << "(dwya + dwyb)/2 == " << dwy << endl; + du.cmplx(0, 1, 0, 0) -= + h * Complex(duy, 0); // modify coeff of 1st chebyshev function T_1(y/h) = y/h for u + du.cmplx(0, 1, 0, 2) -= + h * Complex(dwy, 0); // modify coeff of 1st chebyshev function T_1(y/h) = y/h for w + + du00 = Re(du.profile(0, 0, 0)); + dw00 = Re(du.profile(0, 0, 2)); + du00y = diff(du00); + dw00y = diff(dw00); + cout << "post-mod : " << endl; + cout << "(duya + duyb)/2 == " << (du00y.eval_a() + du00y.eval_b()) / 2 << endl; + cout << "(dwya + dwyb)/2 == " << (dw00y.eval_a() + dw00y.eval_b()) / 2 << endl; + } + } else { // (ilcflags.constraint == BulkVelocity) + // modify du to have zero mean value + if (du.taskid() == du.task_coeff(0, 0)) { + cout << "Modifying du so that it doesn't change mean flow..." << endl; + ChebyCoeff du00 = Re(du.profile(0, 0, 0)); + ChebyCoeff dw00 = Re(du.profile(0, 0, 2)); + Real umean = du00.mean(); + Real wmean = dw00.mean(); + cout << "pre-mod : " << endl; + cout << "u mean == " << du00.mean() << endl; + cout << "w mean == " << dw00.mean() << endl; + du.cmplx(0, 0, 0, 0) -= + Complex(umean, 0); // modify coeff of 0th chebyshev function T_0(y/h) = 1 for u + du.cmplx(0, 0, 0, 2) -= + Complex(wmean, 0); // modify coeff of 0th chebyshev function T_0(y/h) = 1 for w + du00 = Re(du.profile(0, 0, 0)); + dw00 = Re(du.profile(0, 0, 2)); + cout << "post-mod : " << endl; + cout << "u mean == " << du00.mean() << endl; + cout << "w mean == " << dw00.mean() << endl; + } + } + } else { + du = FlowField(duname, cfmpi); + } + if (dtname.length() == 0) { + cout << "Constructing dtemp..." << endl; + bool meanflow_perturb = true; + dtemp.addPerturbations(kxmax, kzmax, 1.0, decay, meanflow_perturb); + } else { + dtemp = FlowField(dtname, cfmpi); + } + + if (ilcflags.symmetries.length() != 0) + project(ilcflags.symmetries, du); + if (ilcflags.tempsymmetries.length() != 0) + project(ilcflags.tempsymmetries, dtemp); + + printout("L2Norm(du) = " + r2s(L2Norm(du))); + printout("rescaling du by eps_du = " + r2s(EPS_du)); + du *= EPS_du / L2Norm(du); + printout("L2Norm(du) = " + r2s(L2Norm(du))); + + printout("L2Norm(dtemp) = " + r2s(L2Norm(dtemp))); + printout("rescaling dtemp by eps_du = " + r2s(EPS_du)); + dtemp *= EPS_du / L2Norm(dtemp); + printout("L2Norm(dtemp) = " + r2s(L2Norm(dtemp))); + + Eigen::VectorXd dx; + field2vector(du, dtemp, dx); + + E->solve(*dsi, x, dx, ilcflags.T, eps); + + fftw_savewisdom(); + } + cfMPI_Finalize(); +} diff --git a/modules/ilc/programs/ilc_findsoln.cpp b/modules/ilc/programs/ilc_findsoln.cpp new file mode 100644 index 00000000..4e7fd02b --- /dev/null +++ b/modules/ilc/programs/ilc_findsoln.cpp @@ -0,0 +1,127 @@ +/** + * Newton search program, like the program for pure shear flows. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include "cfbasics/cfbasics.h" +#include "channelflow/flowfield.h" +#include "modules/ilc/ilcdsi.h" +#include "nsolver/nsolver.h" + +using namespace std; +using namespace chflow; + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + ArgList args(argc, argv, "find an invariant solution in ILC using Newton-Krylov-hookstep algorithm"); + + /** Choose the Newton algorithm to be used. Currently, two options are available: simple Newton without any + * trust region optimization, and Newton with Hookstep (default). For the simple Newton, you can choose either a + * full-space algorithm to solve the Newton equations (-solver "eigen") or between the two iterative algorithms + * GMRES and BiCGStab. Newton-Hookstep requires GMRES. Note that the available parameters depend on your choice + * of the algorithm. + */ + + unique_ptr N; + NewtonSearchFlags searchflags(args); + searchflags.save(searchflags.outdir); + N = unique_ptr(new NewtonAlgorithm(searchflags)); + + ILCFlags ilcflags(args, searchflags.laurette); + TimeStep dt(ilcflags); + + bool Rxsearch, Rzsearch, Tsearch; + Rxsearch = searchflags.xrelative; + Rzsearch = searchflags.zrelative; + Tsearch = searchflags.solntype == PeriodicOrbit ? true : false; + + const bool Tnormalize = Tsearch ? false : true; + + /** Read in remaining arguments */ + + args.section("Program options"); + const string sigmastr = + args.getstr("-sigma", "--sigma", "", "file containing sigma of sigma f^T(u) - u = 0 (default == identity)"); + const Real unormalize = args.getreal("-un", "--unormalize", 0.0, "lower bound in energy for search"); + const int nproc0 = + args.getint("-np0", "--nproc0", 0, "number of MPI-processes for transpose/number of parallel ffts"); + const int nproc1 = args.getint("-np1", "--nproc1", 0, "number of MPI-processes for one fft"); + const bool msinit = + args.getflag("-MSinit", "--MSinitials", "read different files as the initial guesses for different shoots"); + const string uname = args.getstr(2, "", "initial guess for the velocity solution"); + const string tname = args.getstr(1, "", "initial guess for the temperature solution"); + + args.check(); + args.save(); + WriteProcessInfo(argc, argv); + ilcflags.save(); + cout << ilcflags << endl; + + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + + FlowField u(uname, cfmpi); + FlowField temp(tname, cfmpi); + + FieldSymmetry sigma; + if (sigmastr.length() != 0) + sigma = FieldSymmetry(sigmastr); + + /** Construct the dynamical-systems interface object depending on the given parameters. Current options are + * either standard (f(u) via forward time integration) or Laurette (f(u) via Laurettes method) + */ + unique_ptr dsi; + dsi = unique_ptr(new ilcDSI(ilcflags, sigma, 0, dt, Tsearch, Rxsearch, Rzsearch, Tnormalize, unormalize, + u, temp, N->getLogstream())); + + Eigen::VectorXd x_singleShot; + Eigen::VectorXd x; + Eigen::VectorXd yvec; + Eigen::MatrixXd y; + MultishootingDSI* msDSI = N->getMultishootingDSI(); + dsi->makeVectorILC(u, temp, sigma, ilcflags.T, x_singleShot); + msDSI->setDSI(*dsi, x_singleShot.size()); + if (msinit) { + int nSh = msDSI->nShot(); + y.resize(x_singleShot.size(), nSh); + Real Tms = ilcflags.T / nSh; + vector u_ms(nSh); + vector t_ms(nSh); + u_ms[0] = u; + t_ms[0] = temp; + for (int i = 1; i < nSh; i++) { + string uname_ms = "./Multishooting/" + uname + i2s(i); + string tname_ms = "./Multishooting/" + tname + i2s(i); + FlowField ui(uname_ms, cfmpi); + FlowField ti(tname_ms, cfmpi); + u_ms[i] = ui; + t_ms[i] = ti; + } + for (int i = 0; i < nSh; i++) { + dsi->makeVectorILC(u_ms[i], t_ms[i], sigma, Tms, yvec); + y.col(i) = yvec; + } + x = msDSI->toVector(y); + } else { + x = msDSI->makeMSVector(x_singleShot); + } + + int Nunk = x.size(); + int Nunk_total = Nunk; +#ifdef HAVE_MPI + MPI_Allreduce(&Nunk, &Nunk_total, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#endif + cout << Nunk_total << " unknowns" << endl; + + Real residual = 0; + N->solve(*dsi, x, residual); + } + + cfMPI_Finalize(); + + return 0; +} diff --git a/modules/ilc/programs/ilc_simulateflow.cpp b/modules/ilc/programs/ilc_simulateflow.cpp new file mode 100644 index 00000000..1adeb566 --- /dev/null +++ b/modules/ilc/programs/ilc_simulateflow.cpp @@ -0,0 +1,222 @@ +/** + * DNS program, like the program for pure shear flows. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include +#include +#include +#include + +#include "cfbasics/cfvector.h" +#include "cfbasics/mathdefs.h" +#include "channelflow/chebyshev.h" +#include "channelflow/dns.h" +#include "channelflow/flowfield.h" +#include "channelflow/poissonsolver.h" +#include "channelflow/symmetry.h" +#include "channelflow/tausolver.h" +#include "channelflow/utilfuncs.h" +#include "modules/ilc/ilcdsi.h" + +using namespace std; +using namespace chflow; + +string printdiagnostics(FlowField& u, const DNS& dns, Real t, const TimeStep& dt, Real nu, Real umin, bool vardt, + bool pl2norm, bool pchnorm, bool pdissip, bool pshear, bool pdiverge, bool pUbulk, bool pubulk, + bool pdPdx, bool pcfl); + +int main(int argc, char* argv[]) { + cfMPI_Init(&argc, &argv); + { + WriteProcessInfo(argc, argv); + string purpose( + "integrate inclined layer convection (ILC) in the channel flow domain from a given " + "initial condition and save velocity and temperature fields to disk."); + + ArgList args(argc, argv, purpose); + + ILCFlags flags(args); + TimeStep dt(flags); + + args.section("Program options"); + const string outdir = args.getpath("-o", "--outdir", "data/", "output directory"); + const string ulabel = args.getstr("-ul", "--ulabel", "u", "output velocity field prefix"); + const string tlabel = args.getstr("-tl", "--tlabel", "t", "output temperature field prefix"); + const bool savep = args.getflag("-sp", "--savepressure", "save pressure fields"); + + const bool pcfl = args.getflag("-cfl", "--cfl", "print CFL number each dT"); + const bool pl2norm = args.getflag("-l2", "--l2norm", "print L2Norm(u) each dT"); + const bool pchnorm = args.getbool("-ch", "--chebyNorm", true, "print chebyNorm(u) each dT"); + const bool pdissip = args.getflag("-D", "--dissipation", "print dissipation each dT"); + const bool pshear = args.getflag("-I", "--input", "print wall shear power input each dT"); + const bool pdiverge = args.getflag("-dv", "--divergence", "print divergence each dT"); + const bool pubulk = args.getflag("-u", "--ubulk", "print ubulk each dT"); + const bool pUbulk = args.getflag("-Up", "--Ubulk-print", "print Ubulk each dT"); + const bool pdPdx = args.getflag("-p", "--pressure", "print pressure gradient each dT"); + const Real umin = args.getreal("-u", "--umin", 0.0, "stop if chebyNorm(u) < umin"); + + const Real ecfmin = args.getreal("-e", "--ecfmin", 0.0, "stop if Ecf(u) < ecfmin"); + const int saveint = args.getint("-s", "--saveinterval", 1, "save fields every s dT"); + + const int nproc0 = + args.getint("-np0", "--nproc0", 0, "number of MPI-processes for transpose/number of parallel ffts"); + const int nproc1 = args.getint("-np1", "--nproc1", 0, "number of MPI-processes for one fft"); + const string uname = args.getstr(2, "", "initial condition of velocity field"); + const string tname = args.getstr(1, "", "initial condition of temperature field"); + + args.check(); + args.save("./"); + mkdir(outdir); + args.save(outdir); + flags.save(outdir); + + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + + printout("Constructing u,q, and optimizing FFTW..."); + FlowField u(uname, cfmpi); + FlowField temp(tname, cfmpi); + + const int Nx = u.Nx(); + const int Ny = u.Ny(); + const int Nz = u.Nz(); + const Real Lx = u.Lx(); + const Real Lz = u.Lz(); + const Real a = u.a(); + const Real b = u.b(); + + FlowField q(Nx, Ny, Nz, 1, Lx, Lz, a, b, cfmpi); + const bool inttime = + (abs(saveint * dt.dT() - int(saveint * dt.dT())) < 1e-12) && (abs(flags.t0 - int(flags.t0)) < 1e-12) + ? true + : false; + + cout << "Uwall == " << flags.uupperwall << endl; + cout << "Wwall == " << flags.wupperwall << endl; + cout << "DeltaT == " << flags.tlowerwall - flags.tupperwall << endl; + cout << "ilcflags == " << flags << endl; + cout << "constructing ILC DNS..." << endl; + ILC dns({u, temp, q}, flags); + // u.setnu (flags.nu); + + dns.Ubase().save(outdir + "Ubase"); + dns.Wbase().save(outdir + "Wbase"); + dns.Tbase().save(outdir + "Tbase"); + + // ChebyCoeff Ubase = laminarProfile (flags, u.a(), u.b(), u.Ny()); + // Ubase.save("Ubase2"); + + PressureSolver psolver(u, dns.Ubase(), dns.Wbase(), flags.nu, flags.Vsuck, + flags.nonlinearity); // NOT CORRECT FOR ILC + psolver.solve(q, u); + vector fields = {u, temp, q}; + + ios::openmode openflag = (flags.t0 > 0) ? ios::app : ios::out; + + ofstream eout, x0out; + openfile(eout, outdir + "energy.asc", openflag); + eout << ilcfieldstatsheader_t("t", flags) << endl; + + FlowField u0, du, tmp; + + int i = 0; + for (Real t = flags.t0; t <= flags.T; t += dt.dT()) { + string s; + s = printdiagnostics(fields[0], dns, t, dt, flags.nu, umin, dt.variable(), pl2norm, pchnorm, pdissip, + pshear, pdiverge, pUbulk, pubulk, pdPdx, pcfl); + if (ecfmin > 0 && Ecf(fields[0]) < ecfmin) { + cferror("Ecf < ecfmin == " + r2s(ecfmin) + ", exiting"); + } + + cout << s; + s = ilcfieldstats_t(fields[0], fields[1], t, flags); + eout << s << endl; + + if (saveint != 0 && i % saveint == 0) { + fields[0].save(outdir + ulabel + t2s(t, inttime)); + fields[1].save(outdir + tlabel + t2s(t, inttime)); + if (savep) + fields[1].save(outdir + "p" + t2s(t, inttime)); + } + i++; + + dns.advance(fields, dt.n()); + + if (dt.variable() && + dt.adjust(dns.CFL(fields[0]))) // TODO: dt.variable()==true is checked twice here, remove it. + dns.reset_dt(dt); + } + cout << "done!" << endl; + } + cfMPI_Finalize(); +} + +string printdiagnostics(FlowField& u, const DNS& dns, Real t, const TimeStep& dt, Real nu, Real umin, bool vardt, + bool pl2norm, bool pchnorm, bool pdissip, bool pshear, bool pdiverge, bool pUbulk, bool pubulk, + bool pdPdx, bool pcfl) { + // Printing diagnostics + stringstream sout; + sout << " t == " << t << endl; + if (vardt) + sout << " dt == " << Real(dt) << endl; + if (pl2norm) + sout << " L2Norm(u) == " << L2Norm(u) << endl; + + if (pchnorm || umin != 0.0) { + Real chnorm = chebyNorm(u); + sout << "chebyNorm(u) == " << chnorm << endl; + if (chnorm < umin) { + cout << "Exiting: chebyNorm(u) < umin." << endl; + exit(0); + } + } + Real h = 0.5 * (u.b() - u.a()); + u += dns.Ubase(); + if (pl2norm) + sout << " energy(u+U) == " << 0.5 * L2Norm(u) << endl; + if (pdissip) + sout << " dissip(u+U) == " << dissipation(u) << endl; + if (pshear) + sout << "wallshear(u+U) == " << abs(wallshearLower(u)) + abs(wallshearUpper(u)) << endl; + if (pdiverge) + sout << " divNorm(u+U) == " << divNorm(u) << endl; + if (pUbulk) + sout << "mean u+U Ubulk == " << dns.Ubulk() << endl; + u -= dns.Ubase(); + if (u.taskid() == u.task_coeff(0, 0)) { + if (pubulk) + sout << " ubulk == " << Re(u.profile(0, 0, 0)).mean() << endl; + } + if (pdPdx) + sout << " dPdx == " << dns.dPdx() << endl; + if (pl2norm) + sout << " L2Norm(u) == " << L2Norm(u) << endl; + if (pl2norm) + sout << " L2Norm3d(u) == " << L2Norm3d(u) << endl; + + Real cfl = dns.CFL(u); + if (u.taskid() == u.task_coeff(0, 0)) { + ChebyCoeff U = dns.Ubase(); + ChebyCoeff W = dns.Wbase(); + + U.makeSpectral(); + U += Re(u.profile(0, 0, 0)); + Real Ucenter = U.eval(0.5 * (u.a() + u.b())); + Real Uwall = pythag(0.5 * (U.eval_b() - U.eval_a()), 0.5 * (W.eval_b() - W.eval_a())); + Real Umean = U.mean(); + sout << " 1/nu == " << 1 / nu << endl; + sout << " Uwall h/nu == " << Uwall * h / nu << endl; + sout << " Ubulk h/nu == " << dns.Ubulk() * h / nu << endl; + sout << " Umean h/nu == " << Umean << " * " << h << " / " << nu << endl; + sout << " Umean h/nu == " << Umean * h / nu << endl; + sout << " Uparab h/nu == " << 1.5 * dns.Ubulk() * h / nu << endl; + sout << "Ucenter h/nu == " << Ucenter * h / nu << endl; + } + sout << " CFL == " << cfl << endl; + return sout.str(); +} diff --git a/modules/ilc/tests/CMakeLists.txt b/modules/ilc/tests/CMakeLists.txt new file mode 100644 index 00000000..79249907 --- /dev/null +++ b/modules/ilc/tests/CMakeLists.txt @@ -0,0 +1,26 @@ +set(ilc_TESTS ilc_timeIntegrationTest) + +file( + COPY + ${CMAKE_CURRENT_SOURCE_DIR}/data/uinit.nc ${CMAKE_CURRENT_SOURCE_DIR}/data/ufinal.nc + DESTINATION + ${CMAKE_CURRENT_BINARY_DIR}/data +) + +file( + COPY + ${CMAKE_CURRENT_SOURCE_DIR}/data/tinit.nc ${CMAKE_CURRENT_SOURCE_DIR}/data/tfinal.nc + DESTINATION + ${CMAKE_CURRENT_BINARY_DIR}/data +) + +foreach (program ${ilc_TESTS}) + install_channelflow_application(${program} OFF) + target_link_libraries(${program}_app PUBLIC ilc) +endforeach (program) + +add_serial_test(ilc_timeIntegration ilc_timeIntegrationTest) + +if (USE_MPI) + add_mpi_test(mpi_ilc_timeIntegration ilc_timeIntegrationTest) +endif () diff --git a/modules/ilc/tests/data/tfinal.nc b/modules/ilc/tests/data/tfinal.nc new file mode 100644 index 00000000..9bef355d Binary files /dev/null and b/modules/ilc/tests/data/tfinal.nc differ diff --git a/modules/ilc/tests/data/tinit.nc b/modules/ilc/tests/data/tinit.nc new file mode 100644 index 00000000..664c63f6 Binary files /dev/null and b/modules/ilc/tests/data/tinit.nc differ diff --git a/modules/ilc/tests/data/ufinal.nc b/modules/ilc/tests/data/ufinal.nc new file mode 100644 index 00000000..728bdacd Binary files /dev/null and b/modules/ilc/tests/data/ufinal.nc differ diff --git a/modules/ilc/tests/data/uinit.nc b/modules/ilc/tests/data/uinit.nc new file mode 100644 index 00000000..605fe09d Binary files /dev/null and b/modules/ilc/tests/data/uinit.nc differ diff --git a/modules/ilc/tests/ilc_timeIntegrationTest.cpp b/modules/ilc/tests/ilc_timeIntegrationTest.cpp new file mode 100644 index 00000000..e019b993 --- /dev/null +++ b/modules/ilc/tests/ilc_timeIntegrationTest.cpp @@ -0,0 +1,155 @@ +/** + * This time integration test ensures compatibility to previous versions. + * + * This file is a part of channelflow version 2.0, https://channelflow.ch . + * License is GNU GPL version 2 or later: ./LICENSE + * + * Original author: Florian Reetz + */ + +#include +#include "channelflow/flowfield.h" +#include "modules/ilc/ilc.h" +#ifdef HAVE_MPI +#include +#endif +#include "channelflow/cfmpi.h" + +using namespace std; +using namespace chflow; + +int main(int argc, char* argv[]) { + int failure = 0; + bool verbose = true; + cfMPI_Init(&argc, &argv); + { + int taskid = mpirank(); + if (taskid == 0) { + cerr << "time_integrationTest: " << flush; + if (verbose) { + cout << "\n====================================================" << endl; + cout << "time_integrationTest\n\n"; + } + } + + string purpose( + "This program loads initial conditions from the harddisk and integrates\nit for 20 time units. The " + "resulting fields are compared to the final fields \nufinal and tfinal. L2diff decides on success of " + "test."); + ArgList args(argc, argv, purpose); + int nproc0 = args.getint("-np0", "-nproc0", 0, "number of processes for transpose/parallel ffts"); + int nproc1 = args.getint("-np1", "-nproc1", 0, "number of processes for slice fft"); + bool fftwmeasure = args.getflag("-fftwmeasure", "--fftwmeasure", "use fftw_measure instead of fftw_patient"); + Real tol = args.getreal("-t", "--tolerance", 1.0e-14, "max distance allowed for the test to pass"); + string dir = "data"; + args.check(); + + if (taskid == 0 && verbose) + cout << "Creating CfMPI object..." << flush; + CfMPI* cfmpi = &CfMPI::getInstance(nproc0, nproc1); + if (taskid == 0 && verbose) + cout << "done" << endl; + + if (taskid == 0 && verbose) + cout << "Loading FlowFields..." << endl; + FlowField u(dir + "/uinit", cfmpi); + FlowField temp(dir + "/tinit", cfmpi); + if (taskid == 0 && verbose) + cout << "done" << endl; + + if (u.taskid() == 0 && verbose) { + cout << "================================================================\n"; + cout << purpose << endl << endl; + cout << "Distribution of processes is " << u.nproc0() << "x" << u.nproc1() << endl; + } + FlowField utmp(u); + if (fftwmeasure) + utmp.optimizeFFTW(FFTW_MEASURE); + else + utmp.optimizeFFTW(FFTW_PATIENT); + fftw_savewisdom(); + + // Define integration parameters + const int n = 40; // take n steps between printouts + const Real dt = 1.0 / n; // integration timestep + + // Define DNS parameters + ILCFlags flags; + flags.baseflow = + ZeroBase; // should be Laminar base for a better test, but this test was callibrated with ZeroBase + + // Run at default flags. If you change them, recreate the test files. + flags.dt = dt; + flags.verbosity = Silent; + + if (u.taskid() == 0 && verbose) + cout << "Building FlowField q..." << flush; + vector fields = { + u, temp, FlowField(u.Nx(), u.Ny(), u.Nz(), 1, u.Lx(), u.Lz(), u.a(), u.b(), u.cfmpi(), Spectral, Spectral)}; + if (u.taskid() == 0 && verbose) + cout << "done" << endl; + if (u.taskid() == 0 && verbose) + cout << "Building dns..." << flush; + ILC ilc(fields, flags); + if (u.taskid() == 0 && verbose) + cout << "done" << endl; + + Real avtime = 0; + int i = 0; + Real T = 20; + for (Real t = 0; t <= T; t += 1) { + timeval start, end; + gettimeofday(&start, 0); + Real cfl = ilc.CFL(fields[0]); + if (fields[0].taskid() == 0 && verbose) + cout << " t == " << t << endl; + if (fields[0].taskid() == 0 && verbose) + cout << " CFL == " << cfl << endl; + Real l2n = L2Norm(fields[0]); + if (fields[0].taskid() == 0 && verbose) + cout << " L2Norm(u) == " << l2n << endl; + + // Take n steps of length dt + ilc.advance(fields, n); + if (verbose) { + gettimeofday(&end, 0); + Real sec = (Real)(end.tv_sec - start.tv_sec); + Real ms = (((Real)end.tv_usec) - ((Real)start.tv_usec)); + Real timeused = sec + ms / 1000000.; + if (fields[0].taskid() == 0) + cout << "duration for this timeunit: " << timeused << endl; + if (t != 0) { + avtime += timeused; + i++; + } + if (fields[0].taskid() == 0) + cout << endl; + } + } + + FlowField uf(dir + "/ufinal", cfmpi); + FlowField tf(dir + "/tfinal", cfmpi); + Real l2d = L2Dist(uf, fields[0]) + L2Dist(tf, fields[1]); + if (l2d > tol) { + if (fields[0].taskid() == 0) { + if (verbose) + cout << endl << "Final L2Dist: " << l2d << endl; + cerr << "\t** FAIL **" << endl; + cout << "\t** FAIL **" << endl; + } + failure = 1; + } else { + if (fields[0].taskid() == 0) { + if (verbose) + cout << endl << "Final L2Dist: " << l2d << endl; + cerr << "\t pass " << endl; + cout << "\t pass " << endl; + } + } + if (fields[0].taskid() == 0 && verbose) { + cout << "Average time/timeunit: " << avtime / i << "s" << endl; + } + } + cfMPI_Finalize(); + return failure; +}