Compare commits
7 Commits
cjy-oneapi
...
cjy-oneapi
| Author | SHA1 | Date | |
|---|---|---|---|
| 6796384bf4 | |||
| c974a88d6d | |||
| 09ffdb553d | |||
| 699e443c7a | |||
| 24bfa44911 | |||
| 6738854a9d | |||
| 223ec17a54 |
1
.gitignore
vendored
1
.gitignore
vendored
@@ -1,7 +1,6 @@
|
||||
__pycache__
|
||||
GW150914
|
||||
GW150914-origin
|
||||
GW150914-mini
|
||||
docs
|
||||
*.tmp
|
||||
|
||||
|
||||
@@ -16,14 +16,12 @@ import numpy
|
||||
File_directory = "GW150914" ## output file directory
|
||||
Output_directory = "binary_output" ## binary data file directory
|
||||
## The file directory name should not be too long
|
||||
MPI_processes = 8 ## number of mpi processes used in the simulation
|
||||
MPI_processes = 64 ## number of mpi processes used in the simulation
|
||||
|
||||
GPU_Calculation = "no" ## Use GPU or not
|
||||
## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
|
||||
CPU_Part = 1.0
|
||||
GPU_Part = 0.0
|
||||
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
@@ -1,233 +0,0 @@
|
||||
|
||||
#################################################
|
||||
##
|
||||
## This file provides the input parameters required for numerical relativity.
|
||||
## XIAOQU
|
||||
## 2024/03/19 --- 2025/09/14
|
||||
## Modified for GW150914-mini test case
|
||||
##
|
||||
#################################################
|
||||
|
||||
import numpy
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting MPI processes and the output file directory
|
||||
|
||||
File_directory = "GW150914-mini" ## output file directory
|
||||
Output_directory = "binary_output" ## binary data file directory
|
||||
## The file directory name should not be too long
|
||||
MPI_processes = 4 ## number of mpi processes used in the simulation (Reduced for laptop)
|
||||
|
||||
GPU_Calculation = "no" ## Use GPU or not
|
||||
## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
|
||||
CPU_Part = 1.0
|
||||
GPU_Part = 0.0
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the physical system and numerical method
|
||||
|
||||
Symmetry = "equatorial-symmetry" ## Symmetry of System: choose equatorial-symmetry、no-symmetry、octant-symmetry
|
||||
Equation_Class = "BSSN" ## Evolution Equation: choose "BSSN", "BSSN-EScalar", "BSSN-EM", "Z4C"
|
||||
## If "BSSN-EScalar" is chosen, it is necessary to set other parameters below
|
||||
Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical"
|
||||
Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45"
|
||||
Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order"
|
||||
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the time evolutionary information
|
||||
|
||||
Start_Evolution_Time = 0.0 ## start evolution time t0
|
||||
Final_Evolution_Time = 100.0 ## final evolution time t1 (Reduced for quick test)
|
||||
Check_Time = 10.0
|
||||
Dump_Time = 10.0 ## time inteval dT for dumping binary data
|
||||
D2_Dump_Time = 10.0 ## dump the ascii data for 2d surface after dT'
|
||||
Analysis_Time = 1.0 ## dump the puncture position and GW psi4 after dT"
|
||||
Evolution_Step_Number = 10000000 ## stop the calculation after the maximal step number
|
||||
Courant_Factor = 0.5 ## Courant Factor
|
||||
Dissipation = 0.15 ## Kreiss-Oliger Dissipation Strength
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the grid structure
|
||||
|
||||
basic_grid_set = "Patch" ## grid structure: choose "Patch" or "Shell-Patch"
|
||||
grid_center_set = "Cell" ## grid center: chose "Cell" or "Vertex"
|
||||
|
||||
grid_level = 7 ## total number of AMR grid levels (Reduced from 9)
|
||||
static_grid_level = 4 ## number of AMR static grid levels (Reduced from 5)
|
||||
moving_grid_level = grid_level - static_grid_level ## number of AMR moving grid levels
|
||||
|
||||
analysis_level = 0
|
||||
refinement_level = 3 ## time refinement start from this grid level
|
||||
|
||||
largest_box_xyz_max = [320.0, 320.0, 320.0] ## scale of the largest box
|
||||
## not ne cess ary to be cubic for "Patch" grid s tructure
|
||||
## need to be a cubic box for "Shell-Patch" grid structure
|
||||
largest_box_xyz_min = - numpy.array(largest_box_xyz_max)
|
||||
|
||||
static_grid_number = 48 ## grid points of each static AMR grid (in x direction) (Reduced from 96)
|
||||
## (grid points in y and z directions are automatically adjusted)
|
||||
moving_grid_number = 24 ## grid points of each moving AMR grid (Reduced from 48)
|
||||
shell_grid_number = [32, 32, 100] ## grid points of Shell-Patch grid
|
||||
## in (phi, theta, r) direction
|
||||
devide_factor = 2.0 ## resolution between different grid levels dh0/dh1, only support 2.0 now
|
||||
|
||||
|
||||
static_grid_type = 'Linear' ## AMR static grid structure , only supports "Linear"
|
||||
moving_grid_type = 'Linear' ## AMR moving grid structure , only supports "Linear"
|
||||
|
||||
quarter_sphere_number = 48 ## grid number of 1/4 s pher ical surface (Reduced from 96)
|
||||
## (which is needed for evaluating the spherical surface integral)
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the puncture information
|
||||
|
||||
puncture_number = 2
|
||||
|
||||
position_BH = numpy.zeros( (puncture_number, 3) )
|
||||
parameter_BH = numpy.zeros( (puncture_number, 3) )
|
||||
dimensionless_spin_BH = numpy.zeros( (puncture_number, 3) )
|
||||
momentum_BH = numpy.zeros( (puncture_number, 3) )
|
||||
|
||||
puncture_data_set = "Manually" ## Method to give Puncture’s positions and momentum
|
||||
## choose "Manually" or "Automatically-BBH"
|
||||
## Prefer to choose "Manually", because "Automatically-BBH" is developing now
|
||||
|
||||
## initial orbital distance and ellipticity for BBHs system
|
||||
## ( needed for "Automatically-BBH" case , not affect the "Manually" case )
|
||||
Distance = 10.0
|
||||
e0 = 0.0
|
||||
|
||||
## black hole parameter (M Q* a*)
|
||||
parameter_BH[0] = [ 36.0/(36.0+29.0), 0.0, +0.31 ]
|
||||
parameter_BH[1] = [ 29.0/(36.0+29.0), 0.0, -0.46 ]
|
||||
## dimensionless spin in each direction
|
||||
dimensionless_spin_BH[0] = [ 0.0, 0.0, +0.31 ]
|
||||
dimensionless_spin_BH[1] = [ 0.0, 0.0, -0.46 ]
|
||||
|
||||
## use Brugmann's convention
|
||||
## -----0-----> y
|
||||
## - +
|
||||
|
||||
#---------------------------------------------
|
||||
|
||||
## If puncture_data_set is chosen to be "Manually", it is necessary to set the position and momentum of each puncture manually
|
||||
|
||||
## initial position for each puncture
|
||||
position_BH[0] = [ 0.0, 10.0*29.0/(36.0+29.0), 0.0 ]
|
||||
position_BH[1] = [ 0.0, -10.0*36.0/(36.0+29.0), 0.0 ]
|
||||
|
||||
## initial mumentum for each puncture
|
||||
## (needed for "Manually" case, does not affect the "Automatically-BBH" case)
|
||||
momentum_BH[0] = [ -0.09530152296974252, -0.00084541526517121, 0.0 ]
|
||||
momentum_BH[1] = [ +0.09530152296974252, +0.00084541526517121, 0.0 ]
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the gravitational wave information
|
||||
|
||||
GW_L_max = 4 ## maximal L number in gravitational wave
|
||||
GW_M_max = 4 ## maximal M number in gravitational wave
|
||||
Detector_Number = 12 ## number of dector
|
||||
Detector_Rmin = 50.0 ## nearest dector distance
|
||||
Detector_Rmax = 160.0 ## farest dector distance
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Setting the apprent horizon
|
||||
|
||||
AHF_Find = "no" ## whether to find the apparent horizon: choose "yes" or "no"
|
||||
|
||||
AHF_Find_Every = 24
|
||||
AHF_Dump_Time = 20.0
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Other parameters (testing)
|
||||
## Only influence the Equation_Class = "BSSN-EScalar" case
|
||||
|
||||
FR_a2 = 3.0 ## f(R) = R + a2 * R^2
|
||||
FR_l2 = 10000.0
|
||||
FR_phi0 = 0.00005
|
||||
FR_r0 = 120.0
|
||||
FR_sigma0 = 8.0
|
||||
FR_Choice = 2 ## Choice options: 1 2 3 4 5
|
||||
## 1: phi(r) = phi0 * Exp(-(r-r0)**2/sigma0)
|
||||
## V(r) = 0
|
||||
## 2: phi(r) = phi0 * a2^2/(1+a2^2)
|
||||
## V(r) = Exp(-8*Sqrt(PI/3)*phi(r)) * (1-Exp(4*Sqrt(PI/3)*phi(r)))**2 / (32*PI*a2)
|
||||
## 3: Schrodinger-Newton gived by system phi(r)
|
||||
## V(r) = Exp(-8*Sqrt(PI/3)*phi(r)) * (1-Exp(4*Sqrt(PI/3)*phi(r)))**2 / (32*PI*a2)
|
||||
## 4: phi(r) = phi0 * 0.5 * ( tanh((r+r0)/sigma0) - tanh((r-r0)/sigma0) )
|
||||
## V(r) = 0
|
||||
## f(R) = R + a2*R^2 with a2 = +oo
|
||||
## 5: phi(r) = phi0 * Exp(-(r-r0)**2/sigma)
|
||||
## V(r) = 0
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
#################################################
|
||||
|
||||
## Other parameters (testing)
|
||||
## (please do not change if not necessary)
|
||||
|
||||
boundary_choice = "BAM-choice" ## Sommerfeld boundary condition : choose "BAM-choice" or "Shibata-choice"
|
||||
## prefer "BAM-choice"
|
||||
|
||||
gauge_choice = 0 ## gauge choice
|
||||
## 0: B^i gauge
|
||||
## 1: David's puncture gauge
|
||||
## 2: MB B^i gauge
|
||||
## 3: RIT B^i gauge
|
||||
## 4: MB beta gauge
|
||||
## 5: RIT beta gauge
|
||||
## 6: MGB1 B^i gauge
|
||||
## 7: MGB2 B^i gauge
|
||||
## prefer 0 or 1
|
||||
|
||||
tetrad_type = 2 ## tetradtype
|
||||
## v:r; u: phi; w: theta
|
||||
## v^a = (x,y,z)
|
||||
## 0: orthonormal order: v,u,w
|
||||
## v^a = (x,y,z)
|
||||
## m = (phi - i theta)/sqrt(2)
|
||||
## following Frans, Eq.(8) of PRD 75, 124018(2007)
|
||||
## 1: orthonormal order: w,u,v
|
||||
## m = (theta + i phi)/sqrt(2)
|
||||
## following Sperhake, Eq.(3.2) of PRD 85, 124062(2012)
|
||||
## 2: orthonormal order: v,u,w
|
||||
## v_a = (x,y,z)
|
||||
## m = (phi - i theta)/sqrt(2)
|
||||
## following Frans, Eq.(8) of PRD 75, 124018(2007)
|
||||
## this version recommend set to 2
|
||||
## prefer 2
|
||||
|
||||
#################################################
|
||||
@@ -1,224 +0,0 @@
|
||||
##################################################################
|
||||
##
|
||||
## AMSS-NCKU Numerical Relativity Mini Test Program
|
||||
## Author: Assistant (based on Xiaoqu's code)
|
||||
## 2026/01/20
|
||||
##
|
||||
## This script runs a scaled-down version of the GW150914 test case
|
||||
## suitable for laptop testing.
|
||||
##
|
||||
##################################################################
|
||||
|
||||
import os
|
||||
import shutil
|
||||
import sys
|
||||
import time
|
||||
|
||||
# --- Context Manager for Input File Swapping ---
|
||||
class InputFileSwapper:
|
||||
def __init__(self, mini_file="AMSS_NCKU_Input_Mini.py", target_file="AMSS_NCKU_Input.py"):
|
||||
self.mini_file = mini_file
|
||||
self.target_file = target_file
|
||||
self.backup_file = target_file + ".bak"
|
||||
self.swapped = False
|
||||
|
||||
def __enter__(self):
|
||||
print(f"[MiniProgram] Swapping {self.target_file} with {self.mini_file}...")
|
||||
if os.path.exists(self.target_file):
|
||||
shutil.move(self.target_file, self.backup_file)
|
||||
shutil.copy(self.mini_file, self.target_file)
|
||||
self.swapped = True
|
||||
return self
|
||||
|
||||
def __exit__(self, exc_type, exc_value, traceback):
|
||||
if self.swapped:
|
||||
print(f"[MiniProgram] Restoring original {self.target_file}...")
|
||||
os.remove(self.target_file)
|
||||
if os.path.exists(self.backup_file):
|
||||
shutil.move(self.backup_file, self.target_file)
|
||||
|
||||
def main():
|
||||
# Use the swapper to ensure all imported modules see the mini configuration
|
||||
with InputFileSwapper():
|
||||
|
||||
# Import modules AFTER swapping input file
|
||||
try:
|
||||
import AMSS_NCKU_Input as input_data
|
||||
import print_information
|
||||
import setup
|
||||
import numerical_grid
|
||||
import generate_macrodef
|
||||
import makefile_and_run
|
||||
import generate_TwoPuncture_input
|
||||
import renew_puncture_parameter
|
||||
import plot_xiaoqu
|
||||
import plot_GW_strain_amplitude_xiaoqu
|
||||
except ImportError as e:
|
||||
print(f"Error importing modules: {e}")
|
||||
return
|
||||
|
||||
print_information.print_program_introduction()
|
||||
|
||||
print("\n" + "#"*60)
|
||||
print(" RUNNING MINI TEST CASE: GW150914-mini")
|
||||
print("#"*60 + "\n")
|
||||
|
||||
# --- Directory Setup ---
|
||||
File_directory = os.path.join(input_data.File_directory)
|
||||
|
||||
if os.path.exists(File_directory):
|
||||
print(f" Output directory '{File_directory}' exists. Removing for mini test...")
|
||||
shutil.rmtree(File_directory, ignore_errors=True)
|
||||
|
||||
os.mkdir(File_directory)
|
||||
shutil.copy("AMSS_NCKU_Input.py", File_directory) # Copies the current (mini) input
|
||||
|
||||
output_directory = os.path.join(File_directory, "AMSS_NCKU_output")
|
||||
os.mkdir(output_directory)
|
||||
|
||||
binary_results_directory = os.path.join(output_directory, input_data.Output_directory)
|
||||
os.mkdir(binary_results_directory)
|
||||
|
||||
figure_directory = os.path.join(File_directory, "figure")
|
||||
os.mkdir(figure_directory)
|
||||
|
||||
print(" Output directories generated.\n")
|
||||
|
||||
# --- Setup and Input Generation ---
|
||||
setup.print_input_data(File_directory)
|
||||
setup.generate_AMSSNCKU_input()
|
||||
setup.print_puncture_information()
|
||||
|
||||
print("\n Generating AMSS-NCKU input parfile...")
|
||||
numerical_grid.append_AMSSNCKU_cgh_input()
|
||||
|
||||
print("\n Plotting initial grid...")
|
||||
numerical_grid.plot_initial_grid()
|
||||
|
||||
print("\n Generating macro files...")
|
||||
generate_macrodef.generate_macrodef_h()
|
||||
generate_macrodef.generate_macrodef_fh()
|
||||
|
||||
# --- Compilation Preparation ---
|
||||
print("\n Preparing to compile and run...")
|
||||
|
||||
AMSS_NCKU_source_path = "AMSS_NCKU_source"
|
||||
AMSS_NCKU_source_copy = os.path.join(File_directory, "AMSS_NCKU_source_copy")
|
||||
|
||||
if not os.path.exists(AMSS_NCKU_source_path):
|
||||
print(" Error: AMSS_NCKU_source not found! Please run in the project root.")
|
||||
return
|
||||
|
||||
shutil.copytree(AMSS_NCKU_source_path, AMSS_NCKU_source_copy)
|
||||
|
||||
macrodef_h_path = os.path.join(File_directory, "macrodef.h")
|
||||
macrodef_fh_path = os.path.join(File_directory, "macrodef.fh")
|
||||
|
||||
shutil.copy2(macrodef_h_path, AMSS_NCKU_source_copy)
|
||||
shutil.copy2(macrodef_fh_path, AMSS_NCKU_source_copy)
|
||||
|
||||
# --- Compilation ---
|
||||
cwd = os.getcwd()
|
||||
os.chdir(AMSS_NCKU_source_copy)
|
||||
|
||||
print(" Compiling ABE...")
|
||||
makefile_and_run.makefile_ABE()
|
||||
|
||||
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
|
||||
print(" Compiling TwoPunctureABE...")
|
||||
makefile_and_run.makefile_TwoPunctureABE()
|
||||
|
||||
os.chdir(cwd)
|
||||
|
||||
# --- Copy Executables ---
|
||||
if (input_data.GPU_Calculation == "no"):
|
||||
ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
|
||||
else:
|
||||
ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
|
||||
|
||||
if not os.path.exists(ABE_file):
|
||||
print(" Error: ABE executable compilation failed.")
|
||||
return
|
||||
|
||||
shutil.copy2(ABE_file, output_directory)
|
||||
|
||||
TwoPuncture_file = os.path.join(AMSS_NCKU_source_copy, "TwoPunctureABE")
|
||||
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
|
||||
if not os.path.exists(TwoPuncture_file):
|
||||
print(" Error: TwoPunctureABE compilation failed.")
|
||||
return
|
||||
shutil.copy2(TwoPuncture_file, output_directory)
|
||||
|
||||
# --- Execution ---
|
||||
start_time = time.time()
|
||||
|
||||
if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
|
||||
print("\n Generating TwoPuncture input...")
|
||||
generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
|
||||
|
||||
AMSS_NCKU_TwoPuncture_inputfile = 'AMSS-NCKU-TwoPuncture.input'
|
||||
AMSS_NCKU_TwoPuncture_inputfile_path = os.path.join( File_directory, AMSS_NCKU_TwoPuncture_inputfile )
|
||||
shutil.copy2( AMSS_NCKU_TwoPuncture_inputfile_path, os.path.join(output_directory, 'TwoPunctureinput.par') )
|
||||
|
||||
print(" Running TwoPunctureABE...")
|
||||
os.chdir(output_directory)
|
||||
makefile_and_run.run_TwoPunctureABE()
|
||||
os.chdir(cwd)
|
||||
|
||||
# Update Puncture Parameter
|
||||
renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
|
||||
|
||||
AMSS_NCKU_inputfile = 'AMSS-NCKU.input'
|
||||
AMSS_NCKU_inputfile_path = os.path.join(File_directory, AMSS_NCKU_inputfile)
|
||||
shutil.copy2( AMSS_NCKU_inputfile_path, os.path.join(output_directory, 'input.par') )
|
||||
|
||||
print("\n Input files ready. Launching ABE...")
|
||||
|
||||
os.chdir(output_directory)
|
||||
makefile_and_run.run_ABE()
|
||||
os.chdir(cwd)
|
||||
|
||||
end_time = time.time()
|
||||
elapsed_time = end_time - start_time
|
||||
|
||||
# --- Post-processing ---
|
||||
print("\n Copying output files for inspection...")
|
||||
AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "setting.par")
|
||||
if os.path.exists(AMSS_NCKU_error_file_path):
|
||||
shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "AMSSNCKU_setting_parameter") )
|
||||
|
||||
AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "Error.log")
|
||||
if os.path.exists(AMSS_NCKU_error_file_path):
|
||||
shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "Error.log") )
|
||||
|
||||
for fname in ["bssn_BH.dat", "bssn_ADMQs.dat", "bssn_psi4.dat", "bssn_constraint.dat"]:
|
||||
fpath = os.path.join(binary_results_directory, fname)
|
||||
if os.path.exists(fpath):
|
||||
shutil.copy(fpath, os.path.join(output_directory, fname))
|
||||
|
||||
# --- Plotting ---
|
||||
print("\n Plotting results...")
|
||||
try:
|
||||
plot_xiaoqu.generate_puncture_orbit_plot( binary_results_directory, figure_directory )
|
||||
plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
|
||||
plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
|
||||
|
||||
for i in range(input_data.Detector_Number):
|
||||
plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
|
||||
plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
|
||||
|
||||
for i in range(input_data.Detector_Number):
|
||||
plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
|
||||
|
||||
for i in range(input_data.grid_level):
|
||||
plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
|
||||
|
||||
plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
|
||||
except Exception as e:
|
||||
print(f"Warning: Plotting failed: {e}")
|
||||
|
||||
print(f"\n Program Cost = {elapsed_time:.2f} Seconds \n")
|
||||
print(" AMSS-NCKU-Python simulation finished (Mini Test).\n")
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -5,6 +5,7 @@
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <string>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
@@ -60,13 +61,110 @@ TwoPunctures::TwoPunctures(double mp, double mm, double b,
|
||||
F = dvector(0, ntotal - 1);
|
||||
allocate_derivs(&u, ntotal);
|
||||
allocate_derivs(&v, ntotal);
|
||||
|
||||
// Allocate workspace buffers for hot-path allocation elimination
|
||||
int N = maximum3(n1, n2, n3);
|
||||
int maxn = maximum2(n1, n2);
|
||||
|
||||
// LineRelax_be workspace (sized for n2)
|
||||
ws_diag_be = new double[n2];
|
||||
ws_e_be = new double[n2 - 1];
|
||||
ws_f_be = new double[n2 - 1];
|
||||
ws_b_be = new double[n2];
|
||||
ws_x_be = new double[n2];
|
||||
|
||||
// LineRelax_al workspace (sized for n1)
|
||||
ws_diag_al = new double[n1];
|
||||
ws_e_al = new double[n1 - 1];
|
||||
ws_f_al = new double[n1 - 1];
|
||||
ws_b_al = new double[n1];
|
||||
ws_x_al = new double[n1];
|
||||
|
||||
// ThomasAlgorithm workspace (sized for max(n1,n2))
|
||||
ws_thomas_y = new double[maxn];
|
||||
|
||||
// JFD_times_dv workspace (sized for nvar)
|
||||
ws_jfd_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_jfd_dU, nvar);
|
||||
allocate_derivs(&ws_jfd_U, nvar);
|
||||
|
||||
// chebft_Zeros workspace (sized for N+1)
|
||||
ws_cheb_c = dvector(0, N);
|
||||
|
||||
// fourft workspace (sized for N/2+1 each)
|
||||
ws_four_a = dvector(0, N / 2);
|
||||
ws_four_b = dvector(0, N / 2);
|
||||
|
||||
// Derivatives_AB3 workspace
|
||||
ws_deriv_p = dvector(0, N);
|
||||
ws_deriv_dp = dvector(0, N);
|
||||
ws_deriv_d2p = dvector(0, N);
|
||||
ws_deriv_q = dvector(0, N);
|
||||
ws_deriv_dq = dvector(0, N);
|
||||
ws_deriv_r = dvector(0, N);
|
||||
ws_deriv_dr = dvector(0, N);
|
||||
ws_deriv_indx = ivector(0, N);
|
||||
|
||||
// F_of_v workspace
|
||||
ws_fov_sources = new double[n1 * n2 * n3];
|
||||
ws_fov_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_fov_U, nvar);
|
||||
|
||||
// J_times_dv workspace
|
||||
ws_jtdv_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_jtdv_dU, nvar);
|
||||
allocate_derivs(&ws_jtdv_U, nvar);
|
||||
}
|
||||
|
||||
TwoPunctures::~TwoPunctures()
|
||||
{
|
||||
int const nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
|
||||
int N = maximum3(n1, n2, n3);
|
||||
|
||||
free_dvector(F, 0, ntotal - 1);
|
||||
free_derivs(&u, ntotal);
|
||||
free_derivs(&v, ntotal);
|
||||
|
||||
// Free workspace buffers
|
||||
delete[] ws_diag_be;
|
||||
delete[] ws_e_be;
|
||||
delete[] ws_f_be;
|
||||
delete[] ws_b_be;
|
||||
delete[] ws_x_be;
|
||||
|
||||
delete[] ws_diag_al;
|
||||
delete[] ws_e_al;
|
||||
delete[] ws_f_al;
|
||||
delete[] ws_b_al;
|
||||
delete[] ws_x_al;
|
||||
|
||||
delete[] ws_thomas_y;
|
||||
|
||||
free_dvector(ws_jfd_values, 0, nvar - 1);
|
||||
free_derivs(&ws_jfd_dU, nvar);
|
||||
free_derivs(&ws_jfd_U, nvar);
|
||||
|
||||
free_dvector(ws_cheb_c, 0, N);
|
||||
|
||||
free_dvector(ws_four_a, 0, N / 2);
|
||||
free_dvector(ws_four_b, 0, N / 2);
|
||||
|
||||
free_dvector(ws_deriv_p, 0, N);
|
||||
free_dvector(ws_deriv_dp, 0, N);
|
||||
free_dvector(ws_deriv_d2p, 0, N);
|
||||
free_dvector(ws_deriv_q, 0, N);
|
||||
free_dvector(ws_deriv_dq, 0, N);
|
||||
free_dvector(ws_deriv_r, 0, N);
|
||||
free_dvector(ws_deriv_dr, 0, N);
|
||||
free_ivector(ws_deriv_indx, 0, N);
|
||||
|
||||
delete[] ws_fov_sources;
|
||||
free_dvector(ws_fov_values, 0, nvar - 1);
|
||||
free_derivs(&ws_fov_U, nvar);
|
||||
|
||||
free_dvector(ws_jtdv_values, 0, nvar - 1);
|
||||
free_derivs(&ws_jtdv_dU, nvar);
|
||||
free_derivs(&ws_jtdv_U, nvar);
|
||||
}
|
||||
|
||||
void TwoPunctures::Solve()
|
||||
@@ -655,7 +753,7 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
|
||||
int k, j, isignum;
|
||||
double fac, sum, Pion, *c;
|
||||
|
||||
c = dvector(0, n);
|
||||
c = ws_cheb_c;
|
||||
Pion = Pi / n;
|
||||
if (inv == 0)
|
||||
{
|
||||
@@ -686,7 +784,6 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
|
||||
}
|
||||
for (j = 0; j < n; j++)
|
||||
u[j] = c[j];
|
||||
free_dvector(c, 0, n);
|
||||
}
|
||||
|
||||
/* --------------------------------------------------------------------------*/
|
||||
@@ -774,8 +871,8 @@ void TwoPunctures::fourft(double *u, int N, int inv)
|
||||
double x, x1, fac, Pi_fac, *a, *b;
|
||||
|
||||
M = N / 2;
|
||||
a = dvector(0, M);
|
||||
b = dvector(1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/
|
||||
a = ws_four_a;
|
||||
b = ws_four_b - 1; /* offset to match dvector(1,M) indexing */
|
||||
fac = 1. / M;
|
||||
Pi_fac = Pi * fac;
|
||||
if (inv == 0)
|
||||
@@ -824,8 +921,6 @@ void TwoPunctures::fourft(double *u, int N, int inv)
|
||||
iy = -iy;
|
||||
}
|
||||
}
|
||||
free_dvector(a, 0, M);
|
||||
free_dvector(b, 1, M);
|
||||
}
|
||||
|
||||
/* -----------------------------------------*/
|
||||
@@ -1118,14 +1213,14 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
|
||||
double *p, *dp, *d2p, *q, *dq, *r, *dr;
|
||||
|
||||
N = maximum3(n1, n2, n3);
|
||||
p = dvector(0, N);
|
||||
dp = dvector(0, N);
|
||||
d2p = dvector(0, N);
|
||||
q = dvector(0, N);
|
||||
dq = dvector(0, N);
|
||||
r = dvector(0, N);
|
||||
dr = dvector(0, N);
|
||||
indx = ivector(0, N);
|
||||
p = ws_deriv_p;
|
||||
dp = ws_deriv_dp;
|
||||
d2p = ws_deriv_d2p;
|
||||
q = ws_deriv_q;
|
||||
dq = ws_deriv_dq;
|
||||
r = ws_deriv_r;
|
||||
dr = ws_deriv_dr;
|
||||
indx = ws_deriv_indx;
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -1208,14 +1303,6 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
|
||||
}
|
||||
}
|
||||
}
|
||||
free_dvector(p, 0, N);
|
||||
free_dvector(dp, 0, N);
|
||||
free_dvector(d2p, 0, N);
|
||||
free_dvector(q, 0, N);
|
||||
free_dvector(dq, 0, N);
|
||||
free_dvector(r, 0, N);
|
||||
free_dvector(dr, 0, N);
|
||||
free_ivector(indx, 0, N);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
void TwoPunctures::Newton(int const nvar, int const n1, int const n2, int const n3,
|
||||
@@ -1284,10 +1371,11 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
|
||||
derivs U;
|
||||
double *sources;
|
||||
|
||||
values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&U, nvar);
|
||||
values = ws_fov_values;
|
||||
U = ws_fov_U;
|
||||
|
||||
sources = (double *)calloc(n1 * n2 * n3, sizeof(double));
|
||||
sources = ws_fov_sources;
|
||||
memset(sources, 0, n1 * n2 * n3 * sizeof(double));
|
||||
if (0)
|
||||
{
|
||||
double *s_x, *s_y, *s_z;
|
||||
@@ -1442,9 +1530,6 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
|
||||
{
|
||||
fclose(debugfile);
|
||||
}
|
||||
free(sources);
|
||||
free_dvector(values, 0, nvar - 1);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
double TwoPunctures::norm_inf(double const *F, int const ntotal)
|
||||
@@ -1850,11 +1935,12 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
|
||||
|
||||
Derivatives_AB3(nvar, n1, n2, n3, dv);
|
||||
|
||||
values = ws_jtdv_values;
|
||||
dU = ws_jtdv_dU;
|
||||
U = ws_jtdv_U;
|
||||
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&dU, nvar);
|
||||
allocate_derivs(&U, nvar);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
for (k = 0; k < n3; k++)
|
||||
@@ -1908,9 +1994,6 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
|
||||
}
|
||||
}
|
||||
}
|
||||
free_dvector(values, 0, nvar - 1);
|
||||
free_derivs(&dU, nvar);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
@@ -1957,17 +2040,11 @@ void TwoPunctures::LineRelax_be(double *dv,
|
||||
{
|
||||
int j, m, Ic, Ip, Im, col, ivar;
|
||||
|
||||
double *diag = new double[n2];
|
||||
double *e = new double[n2 - 1]; /* above diagonal */
|
||||
double *f = new double[n2 - 1]; /* below diagonal */
|
||||
double *b = new double[n2]; /* rhs */
|
||||
double *x = new double[n2]; /* solution vector */
|
||||
|
||||
// gsl_vector *diag = gsl_vector_alloc(n2);
|
||||
// gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
|
||||
// gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
|
||||
// gsl_vector *b = gsl_vector_alloc(n2); /* rhs */
|
||||
// gsl_vector *x = gsl_vector_alloc(n2); /* solution vector */
|
||||
double *diag = ws_diag_be;
|
||||
double *e = ws_e_be; /* above diagonal */
|
||||
double *f = ws_f_be; /* below diagonal */
|
||||
double *b = ws_b_be; /* rhs */
|
||||
double *x = ws_x_be; /* solution vector */
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -1977,62 +2054,35 @@ void TwoPunctures::LineRelax_be(double *dv,
|
||||
}
|
||||
diag[n2 - 1] = 0;
|
||||
|
||||
// gsl_vector_set_zero(diag);
|
||||
// gsl_vector_set_zero(e);
|
||||
// gsl_vector_set_zero(f);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
Ip = Index(ivar, i, j + 1, k, nvar, n1, n2, n3);
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
Im = Index(ivar, i, j - 1, k, nvar, n1, n2, n3);
|
||||
b[j] = rhs[Ic];
|
||||
// gsl_vector_set(b,j,rhs[Ic]);
|
||||
for (m = 0; m < ncols[Ic]; m++)
|
||||
{
|
||||
col = cols[Ic][m];
|
||||
if (col != Ip && col != Ic && col != Im)
|
||||
b[j] -= JFD[Ic][m] * dv[col];
|
||||
// *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
|
||||
else
|
||||
{
|
||||
if (col == Im && j > 0)
|
||||
f[j - 1] = JFD[Ic][m];
|
||||
// gsl_vector_set(f,j-1,JFD[Ic][m]);
|
||||
if (col == Ic)
|
||||
diag[j] = JFD[Ic][m];
|
||||
// gsl_vector_set(diag,j,JFD[Ic][m]);
|
||||
if (col == Ip && j < n2 - 1)
|
||||
e[j] = JFD[Ic][m];
|
||||
// gsl_vector_set(e,j,JFD[Ic][m]);
|
||||
}
|
||||
}
|
||||
}
|
||||
// A x = b
|
||||
// A = ( d_0 e_0 0 0 )
|
||||
// ( f_0 d_1 e_1 0 )
|
||||
// ( 0 f_1 d_2 e_2 )
|
||||
// ( 0 0 f_2 d_3 )
|
||||
//
|
||||
ThomasAlgorithm(n2, f, diag, e, x, b);
|
||||
// gsl_linalg_solve_tridiag(diag, e, f, b, x);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
dv[Ic] = x[j];
|
||||
// dv[Ic] = gsl_vector_get(x, j);
|
||||
}
|
||||
}
|
||||
|
||||
delete[] diag;
|
||||
delete[] e;
|
||||
delete[] f;
|
||||
delete[] b;
|
||||
delete[] x;
|
||||
// gsl_vector_free(diag);
|
||||
// gsl_vector_free(e);
|
||||
// gsl_vector_free(f);
|
||||
// gsl_vector_free(b);
|
||||
// gsl_vector_free(x);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
@@ -2049,8 +2099,8 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
|
||||
derivs dU, U;
|
||||
|
||||
allocate_derivs(&dU, nvar);
|
||||
allocate_derivs(&U, nvar);
|
||||
dU = ws_jfd_dU;
|
||||
U = ws_jfd_U;
|
||||
|
||||
if (k < 0)
|
||||
k = k + n3;
|
||||
@@ -2168,9 +2218,6 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
LinEquations(A, B, X, R, x, r, phi, y, z, dU, U, values);
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
values[ivar] *= FAC;
|
||||
|
||||
free_derivs(&dU, nvar);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
#undef FAC
|
||||
/*-----------------------------------------------------------*/
|
||||
@@ -2202,17 +2249,11 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
{
|
||||
int i, m, Ic, Ip, Im, col, ivar;
|
||||
|
||||
double *diag = new double[n1];
|
||||
double *e = new double[n1 - 1]; /* above diagonal */
|
||||
double *f = new double[n1 - 1]; /* below diagonal */
|
||||
double *b = new double[n1]; /* rhs */
|
||||
double *x = new double[n1]; /* solution vector */
|
||||
|
||||
// gsl_vector *diag = gsl_vector_alloc(n1);
|
||||
// gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
|
||||
// gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
|
||||
// gsl_vector *b = gsl_vector_alloc(n1); /* rhs */
|
||||
// gsl_vector *x = gsl_vector_alloc(n1); /* solution vector */
|
||||
double *diag = ws_diag_al;
|
||||
double *e = ws_e_al; /* above diagonal */
|
||||
double *f = ws_f_al; /* below diagonal */
|
||||
double *b = ws_b_al; /* rhs */
|
||||
double *x = ws_x_al; /* solution vector */
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -2222,57 +2263,35 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
}
|
||||
diag[n1 - 1] = 0;
|
||||
|
||||
// gsl_vector_set_zero(diag);
|
||||
// gsl_vector_set_zero(e);
|
||||
// gsl_vector_set_zero(f);
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
Ip = Index(ivar, i + 1, j, k, nvar, n1, n2, n3);
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
Im = Index(ivar, i - 1, j, k, nvar, n1, n2, n3);
|
||||
b[i] = rhs[Ic];
|
||||
// gsl_vector_set(b,i,rhs[Ic]);
|
||||
for (m = 0; m < ncols[Ic]; m++)
|
||||
{
|
||||
col = cols[Ic][m];
|
||||
if (col != Ip && col != Ic && col != Im)
|
||||
b[i] -= JFD[Ic][m] * dv[col];
|
||||
// *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
|
||||
else
|
||||
{
|
||||
if (col == Im && i > 0)
|
||||
f[i - 1] = JFD[Ic][m];
|
||||
// gsl_vector_set(f,i-1,JFD[Ic][m]);
|
||||
if (col == Ic)
|
||||
diag[i] = JFD[Ic][m];
|
||||
// gsl_vector_set(diag,i,JFD[Ic][m]);
|
||||
if (col == Ip && i < n1 - 1)
|
||||
e[i] = JFD[Ic][m];
|
||||
// gsl_vector_set(e,i,JFD[Ic][m]);
|
||||
}
|
||||
}
|
||||
}
|
||||
ThomasAlgorithm(n1, f, diag, e, x, b);
|
||||
// gsl_linalg_solve_tridiag(diag, e, f, b, x);
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
dv[Ic] = x[i];
|
||||
// dv[Ic] = gsl_vector_get(x, i);
|
||||
}
|
||||
}
|
||||
|
||||
delete[] diag;
|
||||
delete[] e;
|
||||
delete[] f;
|
||||
delete[] b;
|
||||
delete[] x;
|
||||
|
||||
// gsl_vector_free(diag);
|
||||
// gsl_vector_free(e);
|
||||
// gsl_vector_free(f);
|
||||
// gsl_vector_free(b);
|
||||
// gsl_vector_free(x);
|
||||
}
|
||||
/* -------------------------------------------------------------------------*/
|
||||
// a[N], b[N-1], c[N-1], x[N], q[N]
|
||||
@@ -2284,44 +2303,29 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
//"Parallel Scientific Computing in C++ and MPI" P361
|
||||
void TwoPunctures::ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q)
|
||||
{
|
||||
// In-place Thomas algorithm: uses a[] as d workspace, b[] as l workspace.
|
||||
// c[] is already u (above-diagonal). ws_thomas_y is pre-allocated workspace.
|
||||
int i;
|
||||
double *l, *u, *d, *y;
|
||||
l = new double[N - 1];
|
||||
u = new double[N - 1];
|
||||
d = new double[N];
|
||||
y = new double[N];
|
||||
|
||||
/* LU Decomposition */
|
||||
d[0] = a[0];
|
||||
u[0] = c[0];
|
||||
double *y = ws_thomas_y;
|
||||
|
||||
/* LU Decomposition (in-place: a becomes d, b becomes l) */
|
||||
for (i = 0; i < N - 2; i++)
|
||||
{
|
||||
l[i] = b[i] / d[i];
|
||||
d[i + 1] = a[i + 1] - l[i] * u[i];
|
||||
u[i + 1] = c[i + 1];
|
||||
b[i] = b[i] / a[i];
|
||||
a[i + 1] = a[i + 1] - b[i] * c[i];
|
||||
}
|
||||
|
||||
l[N - 2] = b[N - 2] / d[N - 2];
|
||||
d[N - 1] = a[N - 1] - l[N - 2] * u[N - 2];
|
||||
b[N - 2] = b[N - 2] / a[N - 2];
|
||||
a[N - 1] = a[N - 1] - b[N - 2] * c[N - 2];
|
||||
|
||||
/* Forward Substitution [L][y] = [q] */
|
||||
y[0] = q[0];
|
||||
for (i = 1; i < N; i++)
|
||||
y[i] = q[i] - l[i - 1] * y[i - 1];
|
||||
y[i] = q[i] - b[i - 1] * y[i - 1];
|
||||
|
||||
/* Backward Substitution [U][x] = [y] */
|
||||
x[N - 1] = y[N - 1] / d[N - 1];
|
||||
|
||||
x[N - 1] = y[N - 1] / a[N - 1];
|
||||
for (i = N - 2; i >= 0; i--)
|
||||
x[i] = (y[i] - u[i] * x[i + 1]) / d[i];
|
||||
|
||||
delete[] l;
|
||||
delete[] u;
|
||||
delete[] d;
|
||||
delete[] y;
|
||||
|
||||
return;
|
||||
x[i] = (y[i] - c[i] * x[i + 1]) / a[i];
|
||||
}
|
||||
// --------------------------------------------------------------------------*/
|
||||
// Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are know*/*/
|
||||
|
||||
@@ -42,6 +42,33 @@ private:
|
||||
|
||||
int ntotal;
|
||||
|
||||
// Pre-allocated workspace buffers for hot-path allocation elimination
|
||||
// LineRelax_be workspace (sized for n2)
|
||||
double *ws_diag_be, *ws_e_be, *ws_f_be, *ws_b_be, *ws_x_be;
|
||||
// LineRelax_al workspace (sized for n1)
|
||||
double *ws_diag_al, *ws_e_al, *ws_f_al, *ws_b_al, *ws_x_al;
|
||||
// ThomasAlgorithm workspace (sized for max(n1,n2))
|
||||
double *ws_thomas_y;
|
||||
// JFD_times_dv workspace (sized for nvar)
|
||||
double *ws_jfd_values;
|
||||
derivs ws_jfd_dU, ws_jfd_U;
|
||||
// chebft_Zeros workspace (sized for max(n1,n2,n3)+1)
|
||||
double *ws_cheb_c;
|
||||
// fourft workspace (sized for max(n1,n2,n3)/2+1 each)
|
||||
double *ws_four_a, *ws_four_b;
|
||||
// Derivatives_AB3 workspace
|
||||
double *ws_deriv_p, *ws_deriv_dp, *ws_deriv_d2p;
|
||||
double *ws_deriv_q, *ws_deriv_dq;
|
||||
double *ws_deriv_r, *ws_deriv_dr;
|
||||
int *ws_deriv_indx;
|
||||
// F_of_v workspace
|
||||
double *ws_fov_sources;
|
||||
double *ws_fov_values;
|
||||
derivs ws_fov_U;
|
||||
// J_times_dv workspace
|
||||
double *ws_jtdv_values;
|
||||
derivs ws_jtdv_dU, ws_jtdv_U;
|
||||
|
||||
struct parameters
|
||||
{
|
||||
int nvar, n1, n2, n3;
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1103,6 +1103,103 @@
|
||||
|
||||
end subroutine fderivs
|
||||
!-----------------------------------------------------------------------------
|
||||
! fderivs variant: reuses caller-provided fh work array (memory pool)
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine fderivs_fh(ex,f,fx,fy,fz,X,Y,Z,SYM1,SYM2,SYM3, &
|
||||
symmetry,onoff,fh)
|
||||
implicit none
|
||||
|
||||
integer, intent(in ):: ex(1:3),symmetry,onoff
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: fx,fy,fz
|
||||
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
|
||||
real*8, intent(in ):: SYM1,SYM2,SYM3
|
||||
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh
|
||||
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8, dimension(3) :: SoA
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
|
||||
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
|
||||
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
|
||||
real*8, parameter :: ZEO=0.d0,ONE=1.d0
|
||||
real*8, parameter :: TWO=2.d0,EIT=8.d0
|
||||
real*8, parameter :: F12=1.2d1
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
dZ = Z(2)-Z(1)
|
||||
|
||||
imax = ex(1)
|
||||
jmax = ex(2)
|
||||
kmax = ex(3)
|
||||
|
||||
imin = 1
|
||||
jmin = 1
|
||||
kmin = 1
|
||||
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
|
||||
|
||||
SoA(1) = SYM1
|
||||
SoA(2) = SYM2
|
||||
SoA(3) = SYM3
|
||||
|
||||
call symmetry_bd(2,ex,f,fh,SoA)
|
||||
|
||||
d12dx = ONE/F12/dX
|
||||
d12dy = ONE/F12/dY
|
||||
d12dz = ONE/F12/dZ
|
||||
|
||||
d2dx = ONE/TWO/dX
|
||||
d2dy = ONE/TWO/dY
|
||||
d2dz = ONE/TWO/dZ
|
||||
|
||||
fx = ZEO
|
||||
fy = ZEO
|
||||
fz = ZEO
|
||||
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
#if 0
|
||||
if(i+2 <= imax .and. i-2 >= imin)then
|
||||
fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin)then
|
||||
fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
|
||||
endif
|
||||
if(j+2 <= jmax .and. j-2 >= jmin)then
|
||||
fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
|
||||
elseif(j+1 <= jmax .and. j-1 >= jmin)then
|
||||
fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
|
||||
endif
|
||||
if(k+2 <= kmax .and. k-2 >= kmin)then
|
||||
fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
|
||||
elseif(k+1 <= kmax .and. k-1 >= kmin)then
|
||||
fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
|
||||
endif
|
||||
#else
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
|
||||
fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
|
||||
fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. &
|
||||
j+1 <= jmax .and. j-1 >= jmin .and. &
|
||||
k+1 <= kmax .and. k-1 >= kmin) then
|
||||
fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
|
||||
fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
|
||||
fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
|
||||
endif
|
||||
#endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine fderivs_fh
|
||||
!-----------------------------------------------------------------------------
|
||||
!
|
||||
! single derivatives dx
|
||||
!
|
||||
@@ -1939,31 +2036,30 @@
|
||||
return
|
||||
|
||||
end subroutine fddyz
|
||||
subroutine fderivs_batch4(ex,f1,f2,f3,f4, &
|
||||
f1x,f1y,f1z,f2x,f2y,f2z,f3x,f3y,f3z,f4x,f4y,f4z, &
|
||||
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
|
||||
|
||||
!-----------------------------------------------------------------------------
|
||||
! fdderivs variant: reuses caller-provided fh work array (memory pool)
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine fdderivs_fh(ex,f,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
|
||||
SYM1,SYM2,SYM3,symmetry,onoff,fh)
|
||||
implicit none
|
||||
|
||||
integer, intent(in ):: ex(1:3),symmetry,onoff
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2,f3,f4
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f3x,f3y,f3z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f4x,f4y,f4z
|
||||
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
|
||||
real*8, intent(in ):: SYM1,SYM2,SYM3
|
||||
|
||||
!~~~~~~ other variables
|
||||
integer, intent(in ):: ex(1:3),symmetry,onoff
|
||||
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ):: f
|
||||
real*8, dimension(ex(1),ex(2),ex(3)),intent(out):: fxx,fxy,fxz,fyy,fyz,fzz
|
||||
real*8, intent(in ):: X(ex(1)),Y(ex(2)),Z(ex(3))
|
||||
real*8, intent(in ):: SYM1,SYM2,SYM3
|
||||
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh
|
||||
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2,fh3,fh4
|
||||
real*8, dimension(3) :: SoA
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
|
||||
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
|
||||
real*8 :: Sdxdx,Sdydy,Sdzdz,Fdxdx,Fdydy,Fdzdz
|
||||
real*8 :: Sdxdy,Sdxdz,Sdydz,Fdxdy,Fdxdz,Fdydz
|
||||
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
|
||||
real*8, parameter :: ZEO=0.d0,ONE=1.d0
|
||||
real*8, parameter :: TWO=2.d0,EIT=8.d0
|
||||
real*8, parameter :: F12=1.2d1
|
||||
real*8, parameter :: ZEO=0.d0, ONE=1.d0, TWO=2.d0, F1o4=2.5d-1
|
||||
real*8, parameter :: F8=8.d0, F16=1.6d1, F30=3.d1
|
||||
real*8, parameter :: F1o12=ONE/1.2d1, F1o144=ONE/1.44d2
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
@@ -1984,264 +2080,118 @@
|
||||
SoA(2) = SYM2
|
||||
SoA(3) = SYM3
|
||||
|
||||
call symmetry_bd(2,ex,f1,fh1,SoA)
|
||||
call symmetry_bd(2,ex,f2,fh2,SoA)
|
||||
call symmetry_bd(2,ex,f3,fh3,SoA)
|
||||
call symmetry_bd(2,ex,f4,fh4,SoA)
|
||||
call symmetry_bd(2,ex,f,fh,SoA)
|
||||
|
||||
d12dx = ONE/F12/dX
|
||||
d12dy = ONE/F12/dY
|
||||
d12dz = ONE/F12/dZ
|
||||
Sdxdx = ONE /( dX * dX )
|
||||
Sdydy = ONE /( dY * dY )
|
||||
Sdzdz = ONE /( dZ * dZ )
|
||||
|
||||
d2dx = ONE/TWO/dX
|
||||
d2dy = ONE/TWO/dY
|
||||
d2dz = ONE/TWO/dZ
|
||||
Fdxdx = F1o12 /( dX * dX )
|
||||
Fdydy = F1o12 /( dY * dY )
|
||||
Fdzdz = F1o12 /( dZ * dZ )
|
||||
|
||||
f1x = ZEO; f1y = ZEO; f1z = ZEO
|
||||
f2x = ZEO; f2y = ZEO; f2z = ZEO
|
||||
f3x = ZEO; f3y = ZEO; f3z = ZEO
|
||||
f4x = ZEO; f4y = ZEO; f4z = ZEO
|
||||
Sdxdy = F1o4 /( dX * dY )
|
||||
Sdxdz = F1o4 /( dX * dZ )
|
||||
Sdydz = F1o4 /( dY * dZ )
|
||||
|
||||
Fdxdy = F1o144 /( dX * dY )
|
||||
Fdxdz = F1o144 /( dX * dZ )
|
||||
Fdydz = F1o144 /( dY * dZ )
|
||||
|
||||
fxx = ZEO
|
||||
fyy = ZEO
|
||||
fzz = ZEO
|
||||
fxy = ZEO
|
||||
fxz = ZEO
|
||||
fyz = ZEO
|
||||
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
#if 0
|
||||
if(i+2 <= imax .and. i-2 >= imin)then
|
||||
fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
|
||||
-fh(i+2,j,k)+F16*fh(i+1,j,k) )
|
||||
elseif(i+1 <= imax .and. i-1 >= imin)then
|
||||
fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k))
|
||||
endif
|
||||
if(j+2 <= jmax .and. j-2 >= jmin)then
|
||||
fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
|
||||
-fh(i,j+2,k)+F16*fh(i,j+1,k) )
|
||||
elseif(j+1 <= jmax .and. j-1 >= jmin)then
|
||||
fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k))
|
||||
endif
|
||||
if(k+2 <= kmax .and. k-2 >= kmin)then
|
||||
fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
|
||||
-fh(i,j,k+2)+F16*fh(i,j,k+1) )
|
||||
elseif(k+1 <= kmax .and. k-1 >= kmin)then
|
||||
fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1))
|
||||
endif
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. j+2 <= jmax .and. j-2 >= jmin)then
|
||||
fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) &
|
||||
-F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) &
|
||||
+F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) &
|
||||
- (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. j+1 <= jmax .and. j-1 >= jmin)then
|
||||
fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k))
|
||||
endif
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. k+2 <= kmax .and. k-2 >= kmin)then
|
||||
fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) &
|
||||
-F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) &
|
||||
+F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) &
|
||||
- (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. k+1 <= kmax .and. k-1 >= kmin)then
|
||||
fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
|
||||
endif
|
||||
if(j+2 <= jmax .and. j-2 >= jmin .and. k+2 <= kmax .and. k-2 >= kmin)then
|
||||
fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) &
|
||||
-F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) &
|
||||
+F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) &
|
||||
- (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
|
||||
elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then
|
||||
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
|
||||
endif
|
||||
#else
|
||||
! for bam comparison
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
|
||||
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
|
||||
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
|
||||
|
||||
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
|
||||
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
|
||||
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
|
||||
|
||||
f3x(i,j,k)=d12dx*(fh3(i-2,j,k)-EIT*fh3(i-1,j,k)+EIT*fh3(i+1,j,k)-fh3(i+2,j,k))
|
||||
f3y(i,j,k)=d12dy*(fh3(i,j-2,k)-EIT*fh3(i,j-1,k)+EIT*fh3(i,j+1,k)-fh3(i,j+2,k))
|
||||
f3z(i,j,k)=d12dz*(fh3(i,j,k-2)-EIT*fh3(i,j,k-1)+EIT*fh3(i,j,k+1)-fh3(i,j,k+2))
|
||||
|
||||
f4x(i,j,k)=d12dx*(fh4(i-2,j,k)-EIT*fh4(i-1,j,k)+EIT*fh4(i+1,j,k)-fh4(i+2,j,k))
|
||||
f4y(i,j,k)=d12dy*(fh4(i,j-2,k)-EIT*fh4(i,j-1,k)+EIT*fh4(i,j+1,k)-fh4(i,j+2,k))
|
||||
f4z(i,j,k)=d12dz*(fh4(i,j,k-2)-EIT*fh4(i,j,k-1)+EIT*fh4(i,j,k+1)-fh4(i,j,k+2))
|
||||
fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
|
||||
-fh(i+2,j,k)+F16*fh(i+1,j,k) )
|
||||
fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
|
||||
-fh(i,j+2,k)+F16*fh(i,j+1,k) )
|
||||
fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
|
||||
-fh(i,j,k+2)+F16*fh(i,j,k+1) )
|
||||
fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) &
|
||||
-F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) &
|
||||
+F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) &
|
||||
- (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
|
||||
fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) &
|
||||
-F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) &
|
||||
+F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) &
|
||||
- (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
|
||||
fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) &
|
||||
-F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) &
|
||||
+F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) &
|
||||
- (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. &
|
||||
j+1 <= jmax .and. j-1 >= jmin .and. &
|
||||
k+1 <= kmax .and. k-1 >= kmin) then
|
||||
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
|
||||
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
|
||||
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
|
||||
|
||||
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
|
||||
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
|
||||
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
|
||||
|
||||
f3x(i,j,k)=d2dx*(-fh3(i-1,j,k)+fh3(i+1,j,k))
|
||||
f3y(i,j,k)=d2dy*(-fh3(i,j-1,k)+fh3(i,j+1,k))
|
||||
f3z(i,j,k)=d2dz*(-fh3(i,j,k-1)+fh3(i,j,k+1))
|
||||
|
||||
f4x(i,j,k)=d2dx*(-fh4(i-1,j,k)+fh4(i+1,j,k))
|
||||
f4y(i,j,k)=d2dy*(-fh4(i,j-1,k)+fh4(i,j+1,k))
|
||||
f4z(i,j,k)=d2dz*(-fh4(i,j,k-1)+fh4(i,j,k+1))
|
||||
fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k))
|
||||
fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k))
|
||||
fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1))
|
||||
fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k))
|
||||
fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
|
||||
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine fderivs_batch4
|
||||
!-----------------------------------------------------------------------------
|
||||
! batch first derivatives (3 fields), same symmetry setup
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine fderivs_batch3(ex,f1,f2,f3, &
|
||||
f1x,f1y,f1z,f2x,f2y,f2z,f3x,f3y,f3z, &
|
||||
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
|
||||
implicit none
|
||||
|
||||
integer, intent(in ):: ex(1:3),symmetry,onoff
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2,f3
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f3x,f3y,f3z
|
||||
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
|
||||
real*8, intent(in ):: SYM1,SYM2,SYM3
|
||||
|
||||
!~~~~~~ other variables
|
||||
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2,fh3
|
||||
real*8, dimension(3) :: SoA
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
|
||||
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
|
||||
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
|
||||
real*8, parameter :: ZEO=0.d0,ONE=1.d0
|
||||
real*8, parameter :: TWO=2.d0,EIT=8.d0
|
||||
real*8, parameter :: F12=1.2d1
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
dZ = Z(2)-Z(1)
|
||||
|
||||
imax = ex(1)
|
||||
jmax = ex(2)
|
||||
kmax = ex(3)
|
||||
|
||||
imin = 1
|
||||
jmin = 1
|
||||
kmin = 1
|
||||
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
|
||||
|
||||
SoA(1) = SYM1
|
||||
SoA(2) = SYM2
|
||||
SoA(3) = SYM3
|
||||
|
||||
call symmetry_bd(2,ex,f1,fh1,SoA)
|
||||
call symmetry_bd(2,ex,f2,fh2,SoA)
|
||||
call symmetry_bd(2,ex,f3,fh3,SoA)
|
||||
|
||||
d12dx = ONE/F12/dX
|
||||
d12dy = ONE/F12/dY
|
||||
d12dz = ONE/F12/dZ
|
||||
|
||||
d2dx = ONE/TWO/dX
|
||||
d2dy = ONE/TWO/dY
|
||||
d2dz = ONE/TWO/dZ
|
||||
|
||||
f1x = ZEO; f1y = ZEO; f1z = ZEO
|
||||
f2x = ZEO; f2y = ZEO; f2z = ZEO
|
||||
f3x = ZEO; f3y = ZEO; f3z = ZEO
|
||||
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
|
||||
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
|
||||
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
|
||||
|
||||
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
|
||||
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
|
||||
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
|
||||
|
||||
f3x(i,j,k)=d12dx*(fh3(i-2,j,k)-EIT*fh3(i-1,j,k)+EIT*fh3(i+1,j,k)-fh3(i+2,j,k))
|
||||
f3y(i,j,k)=d12dy*(fh3(i,j-2,k)-EIT*fh3(i,j-1,k)+EIT*fh3(i,j+1,k)-fh3(i,j+2,k))
|
||||
f3z(i,j,k)=d12dz*(fh3(i,j,k-2)-EIT*fh3(i,j,k-1)+EIT*fh3(i,j,k+1)-fh3(i,j,k+2))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. &
|
||||
j+1 <= jmax .and. j-1 >= jmin .and. &
|
||||
k+1 <= kmax .and. k-1 >= kmin) then
|
||||
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
|
||||
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
|
||||
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
|
||||
|
||||
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
|
||||
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
|
||||
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
|
||||
|
||||
f3x(i,j,k)=d2dx*(-fh3(i-1,j,k)+fh3(i+1,j,k))
|
||||
f3y(i,j,k)=d2dy*(-fh3(i,j-1,k)+fh3(i,j+1,k))
|
||||
f3z(i,j,k)=d2dz*(-fh3(i,j,k-1)+fh3(i,j,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine fderivs_batch3
|
||||
!-----------------------------------------------------------------------------
|
||||
! batch first derivatives (2 fields), same symmetry setup
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine fderivs_batch2(ex,f1,f2, &
|
||||
f1x,f1y,f1z,f2x,f2y,f2z, &
|
||||
X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
|
||||
implicit none
|
||||
|
||||
integer, intent(in ):: ex(1:3),symmetry,onoff
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(in ):: f1,f2
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f1x,f1y,f1z
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: f2x,f2y,f2z
|
||||
real*8, intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
|
||||
real*8, intent(in ):: SYM1,SYM2,SYM3
|
||||
|
||||
!~~~~~~ other variables
|
||||
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh1,fh2
|
||||
real*8, dimension(3) :: SoA
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
|
||||
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
|
||||
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
|
||||
real*8, parameter :: ZEO=0.d0,ONE=1.d0
|
||||
real*8, parameter :: TWO=2.d0,EIT=8.d0
|
||||
real*8, parameter :: F12=1.2d1
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
dZ = Z(2)-Z(1)
|
||||
|
||||
imax = ex(1)
|
||||
jmax = ex(2)
|
||||
kmax = ex(3)
|
||||
|
||||
imin = 1
|
||||
jmin = 1
|
||||
kmin = 1
|
||||
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
|
||||
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
|
||||
|
||||
SoA(1) = SYM1
|
||||
SoA(2) = SYM2
|
||||
SoA(3) = SYM3
|
||||
|
||||
call symmetry_bd(2,ex,f1,fh1,SoA)
|
||||
call symmetry_bd(2,ex,f2,fh2,SoA)
|
||||
|
||||
d12dx = ONE/F12/dX
|
||||
d12dy = ONE/F12/dY
|
||||
d12dz = ONE/F12/dZ
|
||||
|
||||
d2dx = ONE/TWO/dX
|
||||
d2dy = ONE/TWO/dY
|
||||
d2dz = ONE/TWO/dZ
|
||||
|
||||
f1x = ZEO; f1y = ZEO; f1z = ZEO
|
||||
f2x = ZEO; f2y = ZEO; f2z = ZEO
|
||||
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
f1x(i,j,k)=d12dx*(fh1(i-2,j,k)-EIT*fh1(i-1,j,k)+EIT*fh1(i+1,j,k)-fh1(i+2,j,k))
|
||||
f1y(i,j,k)=d12dy*(fh1(i,j-2,k)-EIT*fh1(i,j-1,k)+EIT*fh1(i,j+1,k)-fh1(i,j+2,k))
|
||||
f1z(i,j,k)=d12dz*(fh1(i,j,k-2)-EIT*fh1(i,j,k-1)+EIT*fh1(i,j,k+1)-fh1(i,j,k+2))
|
||||
|
||||
f2x(i,j,k)=d12dx*(fh2(i-2,j,k)-EIT*fh2(i-1,j,k)+EIT*fh2(i+1,j,k)-fh2(i+2,j,k))
|
||||
f2y(i,j,k)=d12dy*(fh2(i,j-2,k)-EIT*fh2(i,j-1,k)+EIT*fh2(i,j+1,k)-fh2(i,j+2,k))
|
||||
f2z(i,j,k)=d12dz*(fh2(i,j,k-2)-EIT*fh2(i,j,k-1)+EIT*fh2(i,j,k+1)-fh2(i,j,k+2))
|
||||
elseif(i+1 <= imax .and. i-1 >= imin .and. &
|
||||
j+1 <= jmax .and. j-1 >= jmin .and. &
|
||||
k+1 <= kmax .and. k-1 >= kmin) then
|
||||
f1x(i,j,k)=d2dx*(-fh1(i-1,j,k)+fh1(i+1,j,k))
|
||||
f1y(i,j,k)=d2dy*(-fh1(i,j-1,k)+fh1(i,j+1,k))
|
||||
f1z(i,j,k)=d2dz*(-fh1(i,j,k-1)+fh1(i,j,k+1))
|
||||
|
||||
f2x(i,j,k)=d2dx*(-fh2(i-1,j,k)+fh2(i+1,j,k))
|
||||
f2y(i,j,k)=d2dy*(-fh2(i,j-1,k)+fh2(i,j+1,k))
|
||||
f2z(i,j,k)=d2dz*(-fh2(i,j,k-1)+fh2(i,j,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine fderivs_batch2
|
||||
end subroutine fdderivs_fh
|
||||
|
||||
#elif (ghost_width == 4)
|
||||
! sixth order code
|
||||
@@ -2380,9 +2330,6 @@
|
||||
|
||||
end subroutine fderivs
|
||||
!-----------------------------------------------------------------------------
|
||||
! batch first derivatives (4 fields), same symmetry setup
|
||||
!-----------------------------------------------------------------------------
|
||||
!-----------------------------------------------------------------------------
|
||||
!
|
||||
! single derivatives dx
|
||||
!
|
||||
|
||||
@@ -18,49 +18,61 @@
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz
|
||||
|
||||
!~~~~~~~> Local variable:
|
||||
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: trA,detg
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
|
||||
|
||||
integer :: i,j,k
|
||||
real*8 :: lgxx,lgyy,lgzz,ldetg
|
||||
real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
|
||||
real*8 :: ltrA,lscale
|
||||
real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
|
||||
|
||||
!~~~~~~>
|
||||
|
||||
gxx = dxx + ONE
|
||||
gyy = dyy + ONE
|
||||
gzz = dzz + ONE
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
|
||||
detg = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
|
||||
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
|
||||
gupxx = ( gyy * gzz - gyz * gyz ) / detg
|
||||
gupxy = - ( gxy * gzz - gyz * gxz ) / detg
|
||||
gupxz = ( gxy * gyz - gyy * gxz ) / detg
|
||||
gupyy = ( gxx * gzz - gxz * gxz ) / detg
|
||||
gupyz = - ( gxx * gyz - gxy * gxz ) / detg
|
||||
gupzz = ( gxx * gyy - gxy * gxy ) / detg
|
||||
lgxx = dxx(i,j,k) + ONE
|
||||
lgyy = dyy(i,j,k) + ONE
|
||||
lgzz = dzz(i,j,k) + ONE
|
||||
|
||||
trA = gupxx * Axx + gupyy * Ayy + gupzz * Azz &
|
||||
+ TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
|
||||
ldetg = lgxx * lgyy * lgzz &
|
||||
+ gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) &
|
||||
+ gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) &
|
||||
- gxz(i,j,k) * lgyy * gxz(i,j,k) &
|
||||
- gxy(i,j,k) * gxy(i,j,k) * lgzz &
|
||||
- lgxx * gyz(i,j,k) * gyz(i,j,k)
|
||||
|
||||
Axx = Axx - F1o3 * gxx * trA
|
||||
Axy = Axy - F1o3 * gxy * trA
|
||||
Axz = Axz - F1o3 * gxz * trA
|
||||
Ayy = Ayy - F1o3 * gyy * trA
|
||||
Ayz = Ayz - F1o3 * gyz * trA
|
||||
Azz = Azz - F1o3 * gzz * trA
|
||||
lgupxx = ( lgyy * lgzz - gyz(i,j,k) * gyz(i,j,k) ) / ldetg
|
||||
lgupxy = - ( gxy(i,j,k) * lgzz - gyz(i,j,k) * gxz(i,j,k) ) / ldetg
|
||||
lgupxz = ( gxy(i,j,k) * gyz(i,j,k) - lgyy * gxz(i,j,k) ) / ldetg
|
||||
lgupyy = ( lgxx * lgzz - gxz(i,j,k) * gxz(i,j,k) ) / ldetg
|
||||
lgupyz = - ( lgxx * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / ldetg
|
||||
lgupzz = ( lgxx * lgyy - gxy(i,j,k) * gxy(i,j,k) ) / ldetg
|
||||
|
||||
detg = ONE / ( detg ** F1o3 )
|
||||
|
||||
gxx = gxx * detg
|
||||
gxy = gxy * detg
|
||||
gxz = gxz * detg
|
||||
gyy = gyy * detg
|
||||
gyz = gyz * detg
|
||||
gzz = gzz * detg
|
||||
ltrA = lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
|
||||
+ lgupzz * Azz(i,j,k) &
|
||||
+ TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
|
||||
+ lgupyz * Ayz(i,j,k))
|
||||
|
||||
dxx = gxx - ONE
|
||||
dyy = gyy - ONE
|
||||
dzz = gzz - ONE
|
||||
Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
|
||||
Axy(i,j,k) = Axy(i,j,k) - F1o3 * gxy(i,j,k) * ltrA
|
||||
Axz(i,j,k) = Axz(i,j,k) - F1o3 * gxz(i,j,k) * ltrA
|
||||
Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
|
||||
Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * gyz(i,j,k) * ltrA
|
||||
Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
|
||||
|
||||
lscale = ONE / ( ldetg ** F1o3 )
|
||||
|
||||
dxx(i,j,k) = lgxx * lscale - ONE
|
||||
gxy(i,j,k) = gxy(i,j,k) * lscale
|
||||
gxz(i,j,k) = gxz(i,j,k) * lscale
|
||||
dyy(i,j,k) = lgyy * lscale - ONE
|
||||
gyz(i,j,k) = gyz(i,j,k) * lscale
|
||||
dzz(i,j,k) = lgzz * lscale - ONE
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
@@ -82,51 +94,71 @@
|
||||
real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz
|
||||
|
||||
!~~~~~~~> Local variable:
|
||||
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: trA
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz
|
||||
real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
|
||||
|
||||
integer :: i,j,k
|
||||
real*8 :: lgxx,lgyy,lgzz,lscale
|
||||
real*8 :: lgxy,lgxz,lgyz
|
||||
real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
|
||||
real*8 :: ltrA
|
||||
real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
|
||||
|
||||
!~~~~~~>
|
||||
|
||||
gxx = dxx + ONE
|
||||
gyy = dyy + ONE
|
||||
gzz = dzz + ONE
|
||||
! for g
|
||||
gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
|
||||
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
|
||||
gupzz = ONE / ( gupzz ** F1o3 )
|
||||
|
||||
gxx = gxx * gupzz
|
||||
gxy = gxy * gupzz
|
||||
gxz = gxz * gupzz
|
||||
gyy = gyy * gupzz
|
||||
gyz = gyz * gupzz
|
||||
gzz = gzz * gupzz
|
||||
! for g: normalize determinant first
|
||||
lgxx = dxx(i,j,k) + ONE
|
||||
lgyy = dyy(i,j,k) + ONE
|
||||
lgzz = dzz(i,j,k) + ONE
|
||||
lgxy = gxy(i,j,k)
|
||||
lgxz = gxz(i,j,k)
|
||||
lgyz = gyz(i,j,k)
|
||||
|
||||
dxx = gxx - ONE
|
||||
dyy = gyy - ONE
|
||||
dzz = gzz - ONE
|
||||
! for A
|
||||
lscale = lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
|
||||
+ lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
|
||||
- lgxy * lgxy * lgzz - lgxx * lgyz * lgyz
|
||||
|
||||
gupxx = ( gyy * gzz - gyz * gyz )
|
||||
gupxy = - ( gxy * gzz - gyz * gxz )
|
||||
gupxz = ( gxy * gyz - gyy * gxz )
|
||||
gupyy = ( gxx * gzz - gxz * gxz )
|
||||
gupyz = - ( gxx * gyz - gxy * gxz )
|
||||
gupzz = ( gxx * gyy - gxy * gxy )
|
||||
lscale = ONE / ( lscale ** F1o3 )
|
||||
|
||||
trA = gupxx * Axx + gupyy * Ayy + gupzz * Azz &
|
||||
+ TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
|
||||
lgxx = lgxx * lscale
|
||||
lgxy = lgxy * lscale
|
||||
lgxz = lgxz * lscale
|
||||
lgyy = lgyy * lscale
|
||||
lgyz = lgyz * lscale
|
||||
lgzz = lgzz * lscale
|
||||
|
||||
Axx = Axx - F1o3 * gxx * trA
|
||||
Axy = Axy - F1o3 * gxy * trA
|
||||
Axz = Axz - F1o3 * gxz * trA
|
||||
Ayy = Ayy - F1o3 * gyy * trA
|
||||
Ayz = Ayz - F1o3 * gyz * trA
|
||||
Azz = Azz - F1o3 * gzz * trA
|
||||
dxx(i,j,k) = lgxx - ONE
|
||||
gxy(i,j,k) = lgxy
|
||||
gxz(i,j,k) = lgxz
|
||||
dyy(i,j,k) = lgyy - ONE
|
||||
gyz(i,j,k) = lgyz
|
||||
dzz(i,j,k) = lgzz - ONE
|
||||
|
||||
! for A: trace-free using normalized metric (det=1, no division needed)
|
||||
lgupxx = ( lgyy * lgzz - lgyz * lgyz )
|
||||
lgupxy = - ( lgxy * lgzz - lgyz * lgxz )
|
||||
lgupxz = ( lgxy * lgyz - lgyy * lgxz )
|
||||
lgupyy = ( lgxx * lgzz - lgxz * lgxz )
|
||||
lgupyz = - ( lgxx * lgyz - lgxy * lgxz )
|
||||
lgupzz = ( lgxx * lgyy - lgxy * lgxy )
|
||||
|
||||
ltrA = lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
|
||||
+ lgupzz * Azz(i,j,k) &
|
||||
+ TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
|
||||
+ lgupyz * Ayz(i,j,k))
|
||||
|
||||
Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
|
||||
Axy(i,j,k) = Axy(i,j,k) - F1o3 * lgxy * ltrA
|
||||
Axz(i,j,k) = Axz(i,j,k) - F1o3 * lgxz * ltrA
|
||||
Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
|
||||
Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * lgyz * ltrA
|
||||
Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
|
||||
@@ -324,7 +324,6 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -350,7 +349,6 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -379,7 +377,6 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -886,7 +883,6 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -912,7 +908,6 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -941,7 +936,6 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
funcc = 0.d0
|
||||
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
|
||||
do i=0,ord-1
|
||||
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
|
||||
@@ -1118,64 +1112,65 @@ end subroutine d2dump
|
||||
! Lagrangian polynomial interpolation
|
||||
!------------------------------------------------------------------------------
|
||||
|
||||
subroutine polint(xa,ya,x,y,dy,ordn)
|
||||
|
||||
subroutine polint(xa, ya, x, y, dy, ordn)
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input Parameter:
|
||||
integer,intent(in) :: ordn
|
||||
real*8, dimension(ordn), intent(in) :: xa,ya
|
||||
integer, intent(in) :: ordn
|
||||
real*8, dimension(ordn), intent(in) :: xa, ya
|
||||
real*8, intent(in) :: x
|
||||
real*8, intent(out) :: y,dy
|
||||
real*8, intent(out) :: y, dy
|
||||
|
||||
!~~~~~~> Other parameter:
|
||||
integer :: i, m, ns, n_m
|
||||
real*8, dimension(ordn) :: c, d, ho
|
||||
real*8 :: dif, dift, hp, h, den_val
|
||||
|
||||
integer :: m,n,ns
|
||||
real*8, dimension(ordn) :: c,d,den,ho
|
||||
real*8 :: dif,dift
|
||||
c = ya
|
||||
d = ya
|
||||
ho = xa - x
|
||||
|
||||
!~~~~~~>
|
||||
ns = 1
|
||||
dif = abs(x - xa(1))
|
||||
|
||||
n=ordn
|
||||
m=ordn
|
||||
|
||||
c=ya
|
||||
d=ya
|
||||
ho=xa-x
|
||||
|
||||
ns=1
|
||||
dif=abs(x-xa(1))
|
||||
do m=1,n
|
||||
dift=abs(x-xa(m))
|
||||
if(dift < dif) then
|
||||
ns=m
|
||||
dif=dift
|
||||
end if
|
||||
do i = 2, ordn
|
||||
dift = abs(x - xa(i))
|
||||
if (dift < dif) then
|
||||
ns = i
|
||||
dif = dift
|
||||
end if
|
||||
end do
|
||||
|
||||
y=ya(ns)
|
||||
ns=ns-1
|
||||
do m=1,n-1
|
||||
den(1:n-m)=ho(1:n-m)-ho(1+m:n)
|
||||
if (any(den(1:n-m) == 0.0))then
|
||||
write(*,*) 'failure in polint for point',x
|
||||
write(*,*) 'with input points: ',xa
|
||||
stop
|
||||
endif
|
||||
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
|
||||
d(1:n-m)=ho(1+m:n)*den(1:n-m)
|
||||
c(1:n-m)=ho(1:n-m)*den(1:n-m)
|
||||
if (2*ns < n-m) then
|
||||
dy=c(ns+1)
|
||||
y = ya(ns)
|
||||
ns = ns - 1
|
||||
|
||||
do m = 1, ordn - 1
|
||||
n_m = ordn - m
|
||||
do i = 1, n_m
|
||||
hp = ho(i)
|
||||
h = ho(i+m)
|
||||
den_val = hp - h
|
||||
|
||||
if (den_val == 0.0d0) then
|
||||
write(*,*) 'failure in polint for point',x
|
||||
write(*,*) 'with input points: ',xa
|
||||
stop
|
||||
end if
|
||||
|
||||
den_val = (c(i+1) - d(i)) / den_val
|
||||
|
||||
d(i) = h * den_val
|
||||
c(i) = hp * den_val
|
||||
end do
|
||||
|
||||
if (2 * ns < n_m) then
|
||||
dy = c(ns + 1)
|
||||
else
|
||||
dy=d(ns)
|
||||
ns=ns-1
|
||||
dy = d(ns)
|
||||
ns = ns - 1
|
||||
end if
|
||||
y=y+dy
|
||||
y = y + dy
|
||||
end do
|
||||
|
||||
return
|
||||
|
||||
end subroutine polint
|
||||
!------------------------------------------------------------------------------
|
||||
!
|
||||
@@ -1183,35 +1178,37 @@ end subroutine d2dump
|
||||
!
|
||||
!------------------------------------------------------------------------------
|
||||
subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
|
||||
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input parameters:
|
||||
integer,intent(in) :: ordn
|
||||
real*8, dimension(1:ordn), intent(in) :: x1a,x2a
|
||||
real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
|
||||
real*8, intent(in) :: x1,x2
|
||||
real*8, intent(out) :: y,dy
|
||||
|
||||
!~~~~~~> Other parameters:
|
||||
|
||||
#ifdef POLINT_LEGACY_ORDER
|
||||
integer :: i,m
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8, dimension(ordn) :: yntmp
|
||||
|
||||
m=size(x1a)
|
||||
|
||||
do i=1,m
|
||||
|
||||
yntmp=ya(i,:)
|
||||
call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
|
||||
|
||||
end do
|
||||
|
||||
call polint(x1a,ymtmp,x1,y,dy,ordn)
|
||||
#else
|
||||
integer :: j
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8 :: dy_temp
|
||||
|
||||
do j=1,ordn
|
||||
call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
|
||||
end do
|
||||
call polint(x2a, ymtmp, x2, y, dy, ordn)
|
||||
#endif
|
||||
|
||||
return
|
||||
|
||||
end subroutine polin2
|
||||
!------------------------------------------------------------------------------
|
||||
!
|
||||
@@ -1219,18 +1216,15 @@ end subroutine d2dump
|
||||
!
|
||||
!------------------------------------------------------------------------------
|
||||
subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
|
||||
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input parameters:
|
||||
integer,intent(in) :: ordn
|
||||
real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
|
||||
real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
|
||||
real*8, intent(in) :: x1,x2,x3
|
||||
real*8, intent(out) :: y,dy
|
||||
|
||||
!~~~~~~> Other parameters:
|
||||
|
||||
#ifdef POLINT_LEGACY_ORDER
|
||||
integer :: i,j,m,n
|
||||
real*8, dimension(ordn,ordn) :: yatmp
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
@@ -1239,24 +1233,33 @@ end subroutine d2dump
|
||||
|
||||
m=size(x1a)
|
||||
n=size(x2a)
|
||||
|
||||
do i=1,m
|
||||
do j=1,n
|
||||
|
||||
yqtmp=ya(i,j,:)
|
||||
call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
|
||||
|
||||
end do
|
||||
|
||||
yntmp=yatmp(i,:)
|
||||
call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
|
||||
|
||||
end do
|
||||
|
||||
call polint(x1a,ymtmp,x1,y,dy,ordn)
|
||||
#else
|
||||
integer :: j, k
|
||||
real*8, dimension(ordn,ordn) :: yatmp
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8 :: dy_temp
|
||||
|
||||
do k=1,ordn
|
||||
do j=1,ordn
|
||||
call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
|
||||
end do
|
||||
end do
|
||||
do k=1,ordn
|
||||
call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
|
||||
end do
|
||||
call polint(x3a, ymtmp, x3, y, dy, ordn)
|
||||
#endif
|
||||
|
||||
return
|
||||
|
||||
end subroutine polin3
|
||||
!--------------------------------------------------------------------------------------
|
||||
! calculate L2norm
|
||||
|
||||
@@ -215,6 +215,99 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
|
||||
|
||||
end subroutine kodis
|
||||
|
||||
!-----------------------------------------------------------------------------
|
||||
! kodis variant: reuses caller-provided fh work array (memory pool)
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine kodis_fh(ex,X,Y,Z,f,f_rhs,SoA,Symmetry,eps,fh)
|
||||
|
||||
implicit none
|
||||
! argument variables
|
||||
integer,intent(in) :: Symmetry
|
||||
integer,dimension(3),intent(in)::ex
|
||||
real*8, dimension(1:3), intent(in) :: SoA
|
||||
double precision,intent(in),dimension(ex(1))::X
|
||||
double precision,intent(in),dimension(ex(2))::Y
|
||||
double precision,intent(in),dimension(ex(3))::Z
|
||||
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::f
|
||||
double precision,intent(inout),dimension(ex(1),ex(2),ex(3))::f_rhs
|
||||
real*8,intent(in) :: eps
|
||||
real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh
|
||||
! local variables
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax
|
||||
integer :: i,j,k
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8, parameter :: ONE=1.d0,SIX=6.d0,FIT=1.5d1,TWT=2.d1
|
||||
real*8,parameter::cof=6.4d1 ! 2^6
|
||||
integer, parameter :: NO_SYMM=0, OCTANT=2
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
dZ = Z(2)-Z(1)
|
||||
|
||||
imax = ex(1)
|
||||
jmax = ex(2)
|
||||
kmax = ex(3)
|
||||
|
||||
imin = 1
|
||||
jmin = 1
|
||||
kmin = 1
|
||||
|
||||
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2
|
||||
if(Symmetry == OCTANT .and. dabs(X(1)) < dX) imin = -2
|
||||
if(Symmetry == OCTANT .and. dabs(Y(1)) < dY) jmin = -2
|
||||
|
||||
call symmetry_bd(3,ex,f,fh,SoA)
|
||||
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
|
||||
if(i-3 >= imin .and. i+3 <= imax .and. &
|
||||
j-3 >= jmin .and. j+3 <= jmax .and. &
|
||||
k-3 >= kmin .and. k+3 <= kmax) then
|
||||
#if 0
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dX/cof * ( &
|
||||
(fh(i-3,j,k)+fh(i+3,j,k)) - &
|
||||
SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
|
||||
FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dY/cof * ( &
|
||||
(fh(i,j-3,k)+fh(i,j+3,k)) - &
|
||||
SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
|
||||
FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dZ/cof * ( &
|
||||
(fh(i,j,k-3)+fh(i,j,k+3)) - &
|
||||
SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
|
||||
FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
#else
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/cof *( ( &
|
||||
(fh(i-3,j,k)+fh(i+3,j,k)) - &
|
||||
SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
|
||||
FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
|
||||
TWT* fh(i,j,k) )/dX + &
|
||||
( &
|
||||
(fh(i,j-3,k)+fh(i,j+3,k)) - &
|
||||
SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
|
||||
FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
|
||||
TWT* fh(i,j,k) )/dY + &
|
||||
( &
|
||||
(fh(i,j,k-3)+fh(i,j,k+3)) - &
|
||||
SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
|
||||
FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
|
||||
TWT* fh(i,j,k) )/dZ )
|
||||
#endif
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine kodis_fh
|
||||
|
||||
#elif (ghost_width == 4)
|
||||
! sixth order code
|
||||
!------------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
@@ -487,6 +487,160 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
|
||||
|
||||
end subroutine lopsided
|
||||
|
||||
!-----------------------------------------------------------------------------
|
||||
! lopsided variant: reuses caller-provided fh work array (memory pool)
|
||||
!-----------------------------------------------------------------------------
|
||||
subroutine lopsided_fh(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA,fh)
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input parameters:
|
||||
|
||||
integer, intent(in) :: ex(1:3),Symmetry
|
||||
real*8, intent(in) :: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3))
|
||||
real*8,dimension(ex(1),ex(2),ex(3)),intent(in) :: f,Sfx,Sfy,Sfz
|
||||
|
||||
real*8,dimension(ex(1),ex(2),ex(3)),intent(inout):: f_rhs
|
||||
real*8,dimension(3),intent(in) ::SoA
|
||||
real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh
|
||||
|
||||
!~~~~~~> local variables:
|
||||
integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
|
||||
real*8 :: dX,dY,dZ
|
||||
real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
|
||||
real*8, parameter :: ZEO=0.d0,ONE=1.d0, F3=3.d0
|
||||
real*8, parameter :: TWO=2.d0,F6=6.0d0,F18=1.8d1
|
||||
real*8, parameter :: F12=1.2d1, F10=1.d1,EIT=8.d0
|
||||
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
|
||||
|
||||
dX = X(2)-X(1)
|
||||
dY = Y(2)-Y(1)
|
||||
dZ = Z(2)-Z(1)
|
||||
|
||||
d12dx = ONE/F12/dX
|
||||
d12dy = ONE/F12/dY
|
||||
d12dz = ONE/F12/dZ
|
||||
|
||||
d2dx = ONE/TWO/dX
|
||||
d2dy = ONE/TWO/dY
|
||||
d2dz = ONE/TWO/dZ
|
||||
|
||||
imax = ex(1)
|
||||
jmax = ex(2)
|
||||
kmax = ex(3)
|
||||
|
||||
imin = 1
|
||||
jmin = 1
|
||||
kmin = 1
|
||||
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2
|
||||
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -2
|
||||
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -2
|
||||
|
||||
call symmetry_bd(3,ex,f,fh,SoA)
|
||||
|
||||
! upper bound set ex-1 only for efficiency,
|
||||
! the loop body will set ex 0 also
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
#if 0
|
||||
!! old code - same as original lopsided
|
||||
#else
|
||||
!! new code, 2012dec27, based on bam
|
||||
! x direction
|
||||
if(Sfx(i,j,k) > ZEO)then
|
||||
if(i+3 <= imax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
|
||||
-F6*fh(i+2,j,k)+ fh(i+3,j,k))
|
||||
elseif(i+2 <= imax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
|
||||
elseif(i+1 <= imax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
|
||||
-F6*fh(i-2,j,k)+ fh(i-3,j,k))
|
||||
endif
|
||||
elseif(Sfx(i,j,k) < ZEO)then
|
||||
if(i-3 >= imin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
|
||||
-F6*fh(i-2,j,k)+ fh(i-3,j,k))
|
||||
elseif(i-2 >= imin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
|
||||
elseif(i-1 >= imin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
|
||||
-F6*fh(i+2,j,k)+ fh(i+3,j,k))
|
||||
endif
|
||||
endif
|
||||
|
||||
! y direction
|
||||
if(Sfy(i,j,k) > ZEO)then
|
||||
if(j+3 <= jmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
|
||||
-F6*fh(i,j+2,k)+ fh(i,j+3,k))
|
||||
elseif(j+2 <= jmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfy(i,j,k)*d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
|
||||
elseif(j+1 <= jmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
|
||||
-F6*fh(i,j-2,k)+ fh(i,j-3,k))
|
||||
endif
|
||||
elseif(Sfy(i,j,k) < ZEO)then
|
||||
if(j-3 >= jmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
|
||||
-F6*fh(i,j-2,k)+ fh(i,j-3,k))
|
||||
elseif(j-2 >= jmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfy(i,j,k)*d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
|
||||
elseif(j-1 >= jmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
|
||||
-F6*fh(i,j+2,k)+ fh(i,j+3,k))
|
||||
endif
|
||||
endif
|
||||
|
||||
! z direction
|
||||
if(Sfz(i,j,k) > ZEO)then
|
||||
if(k+3 <= kmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
|
||||
-F6*fh(i,j,k+2)+ fh(i,j,k+3))
|
||||
elseif(k+2 <= kmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfz(i,j,k)*d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
|
||||
elseif(k+1 <= kmax)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
|
||||
-F6*fh(i,j,k-2)+ fh(i,j,k-3))
|
||||
endif
|
||||
elseif(Sfz(i,j,k) < ZEO)then
|
||||
if(k-3 >= kmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)- &
|
||||
Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
|
||||
-F6*fh(i,j,k-2)+ fh(i,j,k-3))
|
||||
elseif(k-2 >= kmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfz(i,j,k)*d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
|
||||
elseif(k-1 >= kmin)then
|
||||
f_rhs(i,j,k)=f_rhs(i,j,k)+ &
|
||||
Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
|
||||
-F6*fh(i,j,k+2)+ fh(i,j,k+3))
|
||||
endif
|
||||
endif
|
||||
#endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
end subroutine lopsided_fh
|
||||
|
||||
#elif (ghost_width == 4)
|
||||
! sixth order code
|
||||
! Compute advection terms in right hand sides of field equations
|
||||
|
||||
@@ -16,10 +16,10 @@ LDLIBS = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore
|
||||
## -fp-model fast=2: Aggressive floating-point optimizations
|
||||
## -fma: Enable fused multiply-add instructions
|
||||
## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues
|
||||
CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma \
|
||||
CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
|
||||
-Dfortran3 -Dnewc -I${MKLROOT}/include
|
||||
f90appflags = -O3 -xHost -fp-model fast=2 -fma \
|
||||
-fpp -I${MKLROOT}/include
|
||||
f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo \
|
||||
-align array64byte -fpp -I${MKLROOT}/include
|
||||
f90 = ifx
|
||||
f77 = ifx
|
||||
CXX = icpx
|
||||
|
||||
@@ -392,17 +392,6 @@ def generate_macrodef_fh():
|
||||
print( "# Finite_Difference_Method #define ghost_width setting error!!!", file=file1 )
|
||||
print( file=file1 )
|
||||
|
||||
# Define macro DEBUG_NAN_CHECK
|
||||
# 0: off (default), 1: on
|
||||
|
||||
debug_nan_check = getattr(input_data, "Debug_NaN_Check", 0)
|
||||
if debug_nan_check:
|
||||
print( "#define DEBUG_NAN_CHECK 1", file=file1 )
|
||||
print( file=file1 )
|
||||
else:
|
||||
print( "#define DEBUG_NAN_CHECK 0", file=file1 )
|
||||
print( file=file1 )
|
||||
|
||||
# Whether to use a shell-patch grid
|
||||
# use shell or not
|
||||
|
||||
@@ -525,9 +514,6 @@ def generate_macrodef_fh():
|
||||
print( " 6th order: 4", file=file1 )
|
||||
print( " 8th order: 5", file=file1 )
|
||||
print( file=file1 )
|
||||
print( "define DEBUG_NAN_CHECK", file=file1 )
|
||||
print( " 0: off (default), 1: on", file=file1 )
|
||||
print( file=file1 )
|
||||
print( "define WithShell", file=file1 )
|
||||
print( " use shell or not", file=file1 )
|
||||
print( file=file1 )
|
||||
|
||||
@@ -35,8 +35,7 @@ Equation_Class = "BSSN" ## Evolution Equation: choose
|
||||
## If "BSSN-EScalar" is chosen, it is necessary to set other parameters below
|
||||
Initial_Data_Method = "Ansorg-TwoPuncture" ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical"
|
||||
Time_Evolution_Method = "runge-kutta-45" ## time evolution method: choose "runge-kutta-45"
|
||||
Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order"
|
||||
Debug_NaN_Check = 0 ## enable NaN checks in compute_rhs_bssn: 0 (off) or 1 (on)
|
||||
Finite_Diffenence_Method = "4th-order" ## finite-difference method: choose "2nd-order", "4th-order", "6th-order", "8th-order"
|
||||
|
||||
#################################################
|
||||
|
||||
|
||||
@@ -15,13 +15,13 @@ import subprocess
|
||||
## taskset ensures all child processes inherit the CPU affinity mask
|
||||
## This forces make and all compiler processes to use only nohz_full cores (4-55, 60-111)
|
||||
## Format: taskset -c 4-55,60-111 ensures processes only run on these cores
|
||||
#NUMACTL_CPU_BIND = "taskset -c 4-55,60-111"
|
||||
NUMACTL_CPU_BIND = ""
|
||||
NUMACTL_CPU_BIND = "taskset -c 16-47,64-95"
|
||||
#NUMACTL_CPU_BIND = "taskset -c 0-111"
|
||||
|
||||
## Build parallelism configuration
|
||||
## Use nohz_full cores (4-55, 60-111) for compilation: 52 + 52 = 104 cores
|
||||
## Set make -j to utilize available cores for faster builds
|
||||
BUILD_JOBS = 14
|
||||
BUILD_JOBS = 64
|
||||
|
||||
|
||||
##################################################################
|
||||
|
||||
Reference in New Issue
Block a user