diff --git a/examples/spectralLinearElastic2Materials.i b/examples/spectralLinearElastic2Materials.i new file mode 100644 index 00000000..7803b3e0 --- /dev/null +++ b/examples/spectralLinearElastic2Materials.i @@ -0,0 +1,173 @@ +[Mesh] + type = MyTRIMMesh + dim = 3 + xmin = -5 + xmax = 5 + ymin = -5 + ymax = 5 + zmin = -5 + zmax = 5 + nx = 32 + ny = 32 + nz = 32 +[] + +[Problem] + type = FFTProblem +[] + +[AuxVariables] + [./epsilon_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./elastic_aux_var] + order = CONSTANT + family = MONOMIAL + + [../] + + [./index_buffer_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./stiffness_ratio_aux] + order = CONSTANT + family = MONOMIAL + [./InitialCondition] + type = SmoothSuperellipsoidIC + x1 = 0 + y1 = 0 + z1 = 0 + a = 2 + b = 2 + c = 2 + n = 2 + int_width = 0 + invalue = 1 + outvalue = 4 + [../] + [../] +[] + +[UserObjects] + # Buffers + [./epsilon] + type = RankTwoTensorFFTWBuffer + [../] + [./stress] + type = RankTwoTensorFFTWBuffer + [../] + # Reciprocal space: Elastic tensor + [./gamma] + type = RankFourTensorFFTWBuffer + [../] + [./elastic] + type = RankFourTensorFFTWBuffer + [../] + [./stiffness_ratio] + type = RealFFTWBuffer + moose_variable = stiffness_ratio_aux + [] + [./index_buffer] + type = RealFFTWBuffer + [] +[] + +[AuxKernels] + [./epsilon_aux] + type = FFTBufferAux + variable = epsilon_aux_var + fft_buffer = epsilon + execute_on = final + component = '0 1' + [../] + + [./stress_aux] + type = FFTBufferAux + variable = stress_aux_var + fft_buffer = stress + execute_on = final + component = '0 1' + [../] + + [./elastic_aux] + type = FFTBufferAux + variable = elastic_aux_var + fft_buffer = elastic + execute_on = final + component = '0 1 0 1' + [../] + + [./stiffness_aux] + type = FFTBufferAux + variable = stiffness_ratio_aux + fft_buffer = stiffness_ratio + execute_on = final + [../] + + [./index_aux] + type = FFTBufferAux + variable = index_buffer_aux_var + fft_buffer = index_buffer + execute_on = final + [../] +[] + +[Executioner] + type = SpectralExecutionerLinearElastic + + time_step = 1.0 + number_iterations = 300 + solver_error = 1.0e-6 + global_strain_tensor = '0.0 0.0 0.0 0.0 0.0 0.01' + young_modulus = 1e4 + poisson_ratio = 0.3 + average_material_factor = 2 +[] + +[VectorPostprocessors] + [./linevaluex] + type = LineValueSampler + variable = 'epsilon_aux_var stress_aux_var' + start_point = '-4.9999999999 0 0' + end_point = '4.9999999999 0 0' + num_points = 101 + sort_by = x + execute_on = final + [../] + [./linevaluey] + type = LineValueSampler + variable = 'epsilon_aux_var stress_aux_var' + start_point = '0 -4.9999999999 0' + end_point = '0 4.9999999999 0' + num_points = 101 + sort_by = y + execute_on = final + [../] + [./linevaluez] + type = LineValueSampler + variable = 'epsilon_aux_var stress_aux_var' + start_point = '0 0 -4.9999999999' + end_point = '0 0 4.9999999999' + num_points = 101 + sort_by = z + execute_on = final + [../] +[] + +[Outputs] + exodus = true + execute_on = 'INITIAL FINAL' + perf_graph = true + [./comp] + type = CSV + [../] +[] diff --git a/examples/spectralLinearElastic2MaterialsTM.i b/examples/spectralLinearElastic2MaterialsTM.i new file mode 100644 index 00000000..50b4e6e9 --- /dev/null +++ b/examples/spectralLinearElastic2MaterialsTM.i @@ -0,0 +1,245 @@ +[GlobalParams] + displacements = 'u_x u_y u_z' +[] + +[Mesh] + [generated_mesh] + type = GeneratedMeshGenerator + dim = 3 + xmin = -5 + xmax = 5 + ymin = -5 + ymax = 5 + zmin = -5 + zmax = 5 + nx = 32 + ny = 32 + nz = 32 + [] + + [cnode] + type = ExtraNodesetGenerator + coord = '0 0 0' + new_boundary = 100 + input = generated_mesh + [] +[] + +[AuxVariables] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] + +[./stress_xy] + order = CONSTANT + family = MONOMIAL +[../] +[./elastic_strain_xy] + order = CONSTANT + family = MONOMIAL +[../] + +[] + +[Variables] + [./u_x] + [../] + [./u_y] + [../] + [./u_z] + [../] + [./global_strain] + order = SIXTH + family = SCALAR + [../] +[] + +[Kernels] + [./TensorMechanics] + [../] +[] + +[Functions] + [./young_modulus_function] + type = ParsedFunction + vars = 'radius ym_min ym_max' + vals = '2.0 1.0 4.0' + value = 'if(sqrt(x*x+y*y+z*z) + +// Forward declarations +class InputParameters; + +/** + * Executioner for linear elastic deformation spectral solver. + */ +class SpectralExecutionerLinearElastic : public SpectralExecutionerBase +{ +public: + static InputParameters validParams(); + + SpectralExecutionerLinearElastic(const InputParameters & parameters); + /** + * Algorithm for incremental solution using forward/backward transforms of Green's function. + */ + virtual void execute() final; + +protected: + /// Time step + const Real _dt; + + /// Number of steps + const unsigned int _nsteps; + + /// Reference Young's modulus + const Real _young_modulus; + + /// Poisson ratio + const Real _poisson_ratio; + + /// Initial strain tensor + RankTwoTensor _initial_strain_tensor; + + /// Average factor to obtain homogeneous material + const Real _average_factor; + + /// User-prescribed error for fixed iteration solver + const Real _solver_error; + + /// Current time + Real _t_current; + +private: + /** + * Helper function to get the linear elastic problem Green's function. + * @param gamma_hat Linear elastic Green operator in Fourier space + * @param elasticity_tensor Elasticity tensor used to replace Gamma for some frequencies + */ + void getGreensFunction(FFTBufferBase & gamma_hat, + const RankFourTensor & elasticity_tensor); + + /** + * Compute initial stress based on homogeneous strain and space-varying elasticity tensor. + * @param epsilon_buffer Small strain FFT buffer + * @param elastic_buffer Space-varying elasticity tensor FFT buffer + */ + FFTBufferBase & + getInitialStress(FFTBufferBase & epsilon_buffer, + const FFTBufferBase & elastic_buffer); + /** + * Initialize epsilon buffer with homogeneous/global strain. + * @param epsilon_buffer Small strain FFT buffer + */ + void fillOutEpsilonBuffer(FFTBufferBase & epsilon_buffer); + + /** + * Initialize epsilon buffer with homogeneous/global strain. + * @param epsilon_buffer Small strain FFT buffer + */ + void advanceReciprocalEpsilon(FFTBufferBase & epsilon_buffer, + const FFTBufferBase & stress_buffer, + const FFTBufferBase & gamma_hat); + /** + * Update real stress buffer for current iteration. + * @param epsilon_buffer Small strain FFT buffer + * @param stress_buffer Stress FFT buffer + * @param elastic_tensor Small strain FFT buffer + */ + void updateRealSigma(const FFTBufferBase & epsilon_buffer, + FFTBufferBase & stress_buffer, + const FFTBufferBase & elastic_tensor); + + /** + * Initialize space-varying elastic coefficient tensor. + * @param ratio_buffer Initial condition ratio to distribute material's stiffnesses + * @param index_buffer Real or integer buffer assigning index to each cell + * @param elastic_tensor_buffer Space-varying elastic coefficient tensor + */ + void filloutElasticTensor(const FFTBufferBase & ratio_buffer, + FFTBufferBase & index_buffer, + FFTBufferBase & elastic_tensor_buffer); + /** + * Convergence check base on global equilibrium. + * @param stress Stress tensor for convergence check + */ + bool hasStressConverged(const FFTBufferBase & stress); +}; diff --git a/include/utils/FFTData.h b/include/utils/FFTData.h index 1f17bab0..0d011cd2 100644 --- a/include/utils/FFTData.h +++ b/include/utils/FFTData.h @@ -48,8 +48,17 @@ class FFTData FFTData & operator*=(Real rhs); FFTData & operator/=(Real rhs); FFTData & operator=(FFTData const & rhs); + FFTData & operator=(const T & rhs); ///@} + /// Templated product of FFTData + template + void setToProductRealSpace(const FFTData & m1, const FFTData & m2); + + // Apply a lambda to the reciprocal space of + template + void applyLambdaReciprocalSpace(T1 lambda); + /// return the number of proper grid cells std::size_t size() const { return _buffer.size(); } @@ -63,3 +72,12 @@ class FFTData /// FFT data buffer std::vector _buffer; }; + +template +template +void +FFTData::applyLambdaReciprocalSpace(T1 lambda) +{ + for (size_t index = 0; index < _buffer.size(); index++) + _buffer[index] = lambda(index); +} diff --git a/src/executioners/SpectralExecutionerLinearElastic.C b/src/executioners/SpectralExecutionerLinearElastic.C new file mode 100644 index 00000000..695ee569 --- /dev/null +++ b/src/executioners/SpectralExecutionerLinearElastic.C @@ -0,0 +1,318 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#include "SpectralExecutionerLinearElastic.h" +#include "InputParameters.h" +#include "RankFourTensor.h" + +registerMooseObject("MagpieApp", SpectralExecutionerLinearElastic); + +InputParameters +SpectralExecutionerLinearElastic::validParams() +{ + InputParameters params = SpectralExecutionerBase::validParams(); + params.addClassDescription("Executioner for spectral solve of diffusion equation"); + params.addParam("time_step", 1.0, "Time step for ODE integration"); + params.addParam("number_iterations", 0, "Maximum number of iterations for solver"); + params.addParam("young_modulus", 1.0, "First parameter for isotropic materials"); + params.addParam("poisson_ratio", 1.0, "Second parameter for isotropic materials"); + params.addParam("average_material_factor", 1.0, "Homogeneized factor for multiphase"); + params.addRequiredParam>( + "global_strain_tensor", + "Vector of values defining the constant applied global strain " + "to add. Components are XX, YY, ZZ, YZ, XZ, XY"); + params.addParam("solver_error", 1.0e-4, "Error for fixed iteration solver"); + + return params; +} + +SpectralExecutionerLinearElastic::SpectralExecutionerLinearElastic( + const InputParameters & parameters) + : SpectralExecutionerBase(parameters), + _dt(getParam("time_step")), + _nsteps(getParam("number_iterations")), + _young_modulus(getParam("young_modulus")), + _poisson_ratio(getParam("poisson_ratio")), + _average_factor(getParam("average_material_factor")), + _solver_error(getParam("solver_error")), + _t_current(0.0) +{ + _initial_strain_tensor.fillFromInputVector(getParam>("global_strain_tensor")); +} + +void +SpectralExecutionerLinearElastic::fillOutEpsilonBuffer( + FFTBufferBase & epsilon_buffer) +{ + epsilon_buffer.realSpace() = _initial_strain_tensor; +} + +void +SpectralExecutionerLinearElastic::getGreensFunction(FFTBufferBase & gamma_hat, + const RankFourTensor & elasticity_tensor) +{ + const int ndim = 3; + std::size_t index = 0; + + const Complex I(0.0, 1.0); + auto & gamma_hat_F = gamma_hat.reciprocalSpace(); + + // Fill fourth-order Green operator based on homogeneous material properties + for (auto ivec : gamma_hat.kTable(0)) + for (auto jvec : gamma_hat.kTable(1)) + for (auto kvec : gamma_hat.kTable(2)) + { + const std::array freq{ivec * I, jvec * I, kvec * I}; + + Real lambda0 = _young_modulus * _average_factor * _poisson_ratio / + ((1 + _poisson_ratio) * (1 - 2 * _poisson_ratio)); + Real nu0 = _young_modulus * _average_factor / (2 * (1 + _poisson_ratio)); + Real constant = (lambda0 + nu0) / (nu0 * (lambda0 + 2.0 * nu0)); + + for (int i = 0; i < ndim; i++) + for (int j = 0; j < ndim; j++) + for (int k = 0; k < ndim; k++) + for (int l = 0; l < ndim; l++) + { + Complex q_square = freq[0] * freq[0] + freq[1] * freq[1] + freq[2] * freq[2]; + if (std::abs(q_square) > 1.0e-12) + { + gamma_hat_F[index](i, j, k, l) = + -1.0 * constant * (freq[i] * freq[j] * freq[k] * freq[l]) / + (q_square * q_square) + + ((i == k) * freq[j] * freq[l] + (j == k) * freq[i] * freq[l] + + (i == l) * freq[j] * freq[k] + (j == l) * freq[i] * freq[k]) / + (4.0 * nu0 * q_square); + } + else + gamma_hat_F[index] = elasticity_tensor.invSymm(); + } + + index++; + } +} + +FFTBufferBase & +SpectralExecutionerLinearElastic::getInitialStress( + FFTBufferBase & epsilon_buffer, + const FFTBufferBase & elastic_tensor) +{ + auto & stress_buffer = getFFTBuffer("stress"); + + // Homogeneous initial state of strain + fillOutEpsilonBuffer(epsilon_buffer); + + // Set real stress buffer to product E * epsilon + stress_buffer.realSpace().setToProductRealSpace(elastic_tensor.realSpace(), + epsilon_buffer.realSpace()); + return stress_buffer; +} + +void +SpectralExecutionerLinearElastic::advanceReciprocalEpsilon( + FFTBufferBase & epsilon_buffer, + const FFTBufferBase & stress_buffer, + const FFTBufferBase & gamma_hat) +{ + Complex I(1.0, 0.0); + + epsilon_buffer.reciprocalSpace().applyLambdaReciprocalSpace( + [&gamma_hat, &stress_buffer, &epsilon_buffer](std::size_t index) { + return epsilon_buffer.reciprocalSpace()[index] - + gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index]; + }); + + // Avoid divide by zero + epsilon_buffer.reciprocalSpace()[0] = _initial_strain_tensor * I; +} + +void +SpectralExecutionerLinearElastic::updateRealSigma( + const FFTBufferBase & epsilon_buffer, + FFTBufferBase & stress_buffer, + const FFTBufferBase & elastic_tensor) +{ + // Set real stress buffer to product E * epsilon + stress_buffer.realSpace().setToProductRealSpace(elastic_tensor.realSpace(), + epsilon_buffer.realSpace()); +} + +void +SpectralExecutionerLinearElastic::filloutElasticTensor( + const FFTBufferBase & ratio_buffer, + FFTBufferBase & index_buffer, + FFTBufferBase & elastic_tensor_buffer) +{ + const auto & grid = elastic_tensor_buffer.grid(); + + int ni = grid[0]; + int nj = grid[1]; + int nk = grid[2]; + + size_t index = 0; + std::vector iso_const(2); + + for (int freq_x = 0; freq_x < ni; ++freq_x) + for (int freq_y = 0; freq_y < nj; ++freq_y) + for (int freq_z = 0; freq_z < nk; ++freq_z) + { + // Define elastic tensor + iso_const[0] = _young_modulus * ratio_buffer.realSpace()[index]; + iso_const[1] = _poisson_ratio; + + elastic_tensor_buffer.realSpace()[index].fillFromInputVector( + iso_const, RankFourTensor::symmetric_isotropic_E_nu); + index_buffer.realSpace()[index] = index; + index++; + } +} + +bool +SpectralExecutionerLinearElastic::hasStressConverged(const FFTBufferBase & stress) +{ + + const auto & grid = stress.grid(); + + const int ni = grid[0]; + const int nj = grid[1]; + const int nk = (grid[2] >> 1) + 1; + + std::size_t index = 0; + + const auto & ivec = stress.kTable(0); + const auto & jvec = stress.kTable(1); + const auto & kvec = stress.kTable(2); + + const std::vector grid_vector = stress.grid(); + + const Complex I(0.0, 1.0); + + Complex error_n(0.0, 0.0); + Complex error_0(0.0, 0.0); + + // Error: Ensure overall equilibrium + for (int freq_x = 0; freq_x < ni; ++freq_x) + for (int freq_y = 0; freq_y < nj; ++freq_y) + for (int freq_z = 0; freq_z < nk; ++freq_z) + { + const ComplexVectorValue freq{ivec[freq_x] * I, jvec[freq_y] * I, kvec[freq_z] * I}; + + ComplexVectorValue kvector_stress; + kvector_stress = stress.reciprocalSpace()[index] * freq; + + error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) + + kvector_stress(2) * kvector_stress(2); + + if (freq_x == 0 && freq_y == 0 && freq_z == 0) + error_0 = stress.reciprocalSpace()[0].L2norm() * stress.reciprocalSpace()[0].L2norm(); + + index++; + } + + Real iteration_error = std::sqrt(std::norm(error_n)) / std::sqrt(std::norm(error_0)); + Moose::out << "Iteration error: " << iteration_error << "\n"; + + if (iteration_error > _solver_error) + return false; + else + return true; +} + +void +SpectralExecutionerLinearElastic::execute() +{ + unsigned int thisStep = 0; + + _time_step = thisStep; + _time = _time_step; + _fe_problem.outputStep(EXEC_INITIAL); + _fe_problem.advanceState(); + + auto & epsilon_buffer = getFFTBuffer("epsilon"); + if (epsilon_buffer.dim() != 3) + mooseError("Error: Problem dimension not implemented in SpectralExecutionerLinearElastic."); + + auto & ratio_buffer = getFFTBuffer("stiffness_ratio"); + auto & elastic_tensor_buffer = getFFTBuffer("elastic"); + auto & index_buffer = getFFTBuffer("index_buffer"); + + // Fill out space-varying elasticity tensor + filloutElasticTensor(ratio_buffer, index_buffer, elastic_tensor_buffer); + + // Get corresponding initial stress (also fill out epsilon_buffer with initial strain) + auto & stress_buffer = getInitialStress(epsilon_buffer, elastic_tensor_buffer); + + // Get specific Green's function + auto & gamma_hat = getFFTBuffer("gamma"); + + thisStep++; + _t_current += _dt; + _time_step = thisStep; + + _fe_problem.execute(EXEC_FINAL); + _time = _t_current; + + Moose::out << "_t_current: " << _t_current << ". \n"; + + _fe_problem.outputStep(EXEC_FINAL); + + // Fill out a homogeneous elasticity tensor with some average properties + std::vector iso_const(2); + iso_const[0] = _young_modulus * _average_factor; + iso_const[1] = _poisson_ratio; + RankFourTensor elasticity_homo; + elasticity_homo.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic_E_nu); + getGreensFunction(gamma_hat, elasticity_homo); + + // Our FFTW "many" plans do not preserve the input, so explicit copies are made + FFTData stress_buffer_backup_real = stress_buffer.realSpace(); + epsilon_buffer.forward(); + + FFTData epsilon_buffer_backup_reciprocal = epsilon_buffer.reciprocalSpace(); + + bool is_converged = false; + + for (unsigned int step_no = 0; step_no < _nsteps; step_no++) + { + // Update sigma in the real space + updateRealSigma(epsilon_buffer, stress_buffer, elastic_tensor_buffer); + + stress_buffer_backup_real = stress_buffer.realSpace(); + stress_buffer.forwardRaw(); + stress_buffer.reciprocalSpace() *= stress_buffer.backwardScale(); + stress_buffer.realSpace() = stress_buffer_backup_real; + + // Convergence check: Ensure global equilibrium + is_converged = hasStressConverged(stress_buffer); + + // Compute new strain tensor in Fourier space + epsilon_buffer.reciprocalSpace() = epsilon_buffer_backup_reciprocal; + advanceReciprocalEpsilon(epsilon_buffer, stress_buffer, gamma_hat); + + // Cache reciprocal epsilon to avoid being overwritten upon backward (inverse) operation + epsilon_buffer_backup_reciprocal = epsilon_buffer.reciprocalSpace(); + epsilon_buffer.backwardRaw(); + + thisStep++; + _t_current += _dt; + _time_step = thisStep; + + _fe_problem.execute(EXEC_FINAL); + _time = _t_current; + + Moose::out << "_t_current: " << _t_current << ". \n"; + + _fe_problem.outputStep(EXEC_FINAL); + + if (is_converged) + break; + + if (step_no != _nsteps - 1) + _fe_problem.advanceState(); + } +} diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C index 33a16e4d..fe605e7c 100644 --- a/src/userobjects/FFTBufferBase.C +++ b/src/userobjects/FFTBufferBase.C @@ -120,6 +120,9 @@ FFTBufferBase::FFTBufferBase(const InputParameters & parameters) for (int j = 0; j < _grid[i]; ++j) _k_table[i][j] = 2.0 * libMesh::pi * ((j < (_grid[i] >> 1) + 1) ? Real(j) : (Real(j) - _grid[i])) / _box_size(i); + + if (i == _dim - 1) + _k_table[i].resize((_grid[i] >> 1) + 1); } _real_space_data.resize(real_space_buffer_size); _reciprocal_space_data.resize(reciprocal_space_buffer_size); diff --git a/src/userobjects/FFTWBufferBase.C b/src/userobjects/FFTWBufferBase.C index 7ed4cb60..4fc25b10 100644 --- a/src/userobjects/FFTWBufferBase.C +++ b/src/userobjects/FFTWBufferBase.C @@ -17,13 +17,14 @@ template FFTWBufferBase::FFTWBufferBase(const InputParameters & parameters) : FFTBufferBase(parameters), - _perf_plan(this->registerTimedSection("fftw_plan_r2r", 2)), + _perf_plan(this->registerTimedSection("fftw_plan_r2c", 2)), _perf_fft(this->registerTimedSection("fftw_execute", 2)) { // create plans { TIME_SECTION(_perf_plan); + // Note: These plans do not preserve the input. Additional caching may be required. _forward_plan = fftw_plan_many_dft_r2c(_dim, _grid.data(), diff --git a/src/utils/ComplexTypes.C b/src/utils/ComplexTypes.C index 9556ff51..e309f8d7 100644 --- a/src/utils/ComplexTypes.C +++ b/src/utils/ComplexTypes.C @@ -45,3 +45,6 @@ template class VectorValue; template class RankTwoTensorTempl; template class RankThreeTensorTempl; template class RankFourTensorTempl; + +template RankTwoTensorTempl RankFourTensorTempl:: +operator*(const RankTwoTensorTempl & a) const; diff --git a/src/utils/FFTData.C b/src/utils/FFTData.C index 0f162069..0f6981b8 100644 --- a/src/utils/FFTData.C +++ b/src/utils/FFTData.C @@ -72,6 +72,25 @@ FFTData::operator=(FFTData const & rhs) return *this; } +template +FFTData & +FFTData::operator=(const T & rhs) +{ + std::fill(_buffer.begin(), _buffer.end(), rhs); + return *this; +} + +template +template +void +FFTData::setToProductRealSpace(const FFTData & m1, const FFTData & m2) +{ + mooseAssert(grid.size() == 3, "Product defined for 3 dimensions"); + + for (unsigned int index = 0; index < _buffer.size(); index++) + (*this)[index] = m1[index] * m2[index]; +} + template <> void * FFTData::start(std::size_t i) @@ -191,3 +210,6 @@ template class FFTData; template class FFTData; template class FFTData; template class FFTData; + +template void FFTData::setToProductRealSpace(const FFTData & m1, + const FFTData & m2); diff --git a/tests/spectral/gold/spectralLinearElastic2Materials_out.e b/tests/spectral/gold/spectralLinearElastic2Materials_out.e new file mode 100644 index 00000000..84b33602 Binary files /dev/null and b/tests/spectral/gold/spectralLinearElastic2Materials_out.e differ diff --git a/tests/spectral/spectralLinearElastic2Materials.i b/tests/spectral/spectralLinearElastic2Materials.i new file mode 100644 index 00000000..21a76429 --- /dev/null +++ b/tests/spectral/spectralLinearElastic2Materials.i @@ -0,0 +1,165 @@ +[Mesh] + type = MyTRIMMesh + dim = 3 + xmin = -5 + xmax = 5 + ymin = -5 + ymax = 5 + zmin = -5 + zmax = 5 + nx = 32 + ny = 32 + nz = 32 +[] + +[Problem] + type = FFTProblem +[] + +[AuxVariables] + [./epsilon_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./elastic_aux_var] + order = CONSTANT + family = MONOMIAL + + [../] + + [./index_buffer_aux_var] + order = CONSTANT + family = MONOMIAL + [../] + + [./stiffness_ratio_aux] + order = CONSTANT + family = MONOMIAL + [./InitialCondition] + type = SmoothSuperellipsoidIC + x1 = 0 + y1 = 0 + z1 = 0 + a = 2 + b = 2 + c = 2 + n = 2 + int_width = 0 + invalue = 1 + outvalue = 4 + [../] + [../] +[] + +[UserObjects] + # Buffers + [./epsilon] + type = RankTwoTensorFFTWBuffer + [../] + [./stress] + type = RankTwoTensorFFTWBuffer + [../] + # Reciprocal space: Elastic tensor + [./gamma] + type = RankFourTensorFFTWBuffer + [../] + [./elastic] + type = RankFourTensorFFTWBuffer + [../] + [./stiffness_ratio] + type = RealFFTWBuffer + moose_variable = stiffness_ratio_aux + [] + [./index_buffer] + type = RealFFTWBuffer + [] +[] + +[AuxKernels] + [./epsilon_aux] + type = FFTBufferAux + variable = epsilon_aux_var + fft_buffer = epsilon + execute_on = final + component = '0 1' + [../] + + [./stress_aux] + type = FFTBufferAux + variable = stress_aux_var + fft_buffer = stress + execute_on = final + component = '0 1' + [../] + + [./elastic_aux] + type = FFTBufferAux + variable = elastic_aux_var + fft_buffer = elastic + execute_on = final + component = '0 1 0 1' + [../] + + [./stiffness_aux] + type = FFTBufferAux + variable = stiffness_ratio_aux + fft_buffer = stiffness_ratio + execute_on = final + [../] + + [./index_aux] + type = FFTBufferAux + variable = index_buffer_aux_var + fft_buffer = index_buffer + execute_on = final + [../] +[] + +[Executioner] + type = SpectralExecutionerLinearElastic + + number_iterations = 300 + solver_error = 3.0e-5 + global_strain_tensor = '0.0 0.0 0.0 0.0 0.0 0.01' + + time_step = 1.0 + young_modulus = 1e4 + poisson_ratio = 0.3 + average_material_factor = 2 +[] + +[VectorPostprocessors] + [./linevaluex] + type = LineValueSampler + variable = 'epsilon_aux_var stress_aux_var' + start_point = '-4.9999999999 0 0' + end_point = '4.9999999999 0 0' + num_points = 11 + sort_by = x + execute_on = final + [../] + [./linevaluez] + type = LineValueSampler + variable = 'epsilon_aux_var stress_aux_var' + start_point = '0 0 -4.9999999999' + end_point = '0 0 4.9999999999' + num_points = 11 + sort_by = z + execute_on = final + [../] +[] + +[Outputs] + exodus = true + execute_on = 'INITIAL FINAL' + perf_graph = true + [./comp] + type = CSV + [../] +[] diff --git a/tests/spectral/tests b/tests/spectral/tests index 5c73cda2..55a1a9ee 100644 --- a/tests/spectral/tests +++ b/tests/spectral/tests @@ -22,7 +22,15 @@ abs_zero = 1e-6 max_parallel = 1 max_threads = 1 - # required_objects = NEMLStress allow_test_objects = true [../] + [./spectralLinearElastic2Materials] + input = 'spectralLinearElastic2Materials.i' + type = 'Exodiff' + exodiff = 'spectralLinearElastic2Materials_out.e' + abs_zero = 1e-6 + max_parallel = 1 + max_threads = 1 + valgrind = 'HEAVY' # timing out + [../] []