""" Compute diffusion-limited oxidation and ageing in rigid medium. """ from __future__ import absolute_import from sfepy import data_dir from sfepy.base.base import output filename_mesh = data_dir + '/meshes/3d/cylinder.mesh' t0 = 0.0 t1 = 0.1 n_step = 11 def get_conductivity(ts, coors, problem, equations=None, mode=None, **kwargs): """ Calculates the conductivity as 2+10*T and returns it. This relation results in larger T gradients where T is small. """ if mode == 'qp': # T-field values in quadrature points coordinates given by integral i # - they are the same as in `coors` argument. T_values = problem.evaluate('ev_volume_integrate.i.Omega(T)', mode='qp', verbose=False) val = 2 + 10 * (T_values + 2) output('conductivity: min:', val.min(), 'max:', val.max()) val.shape = (val.shape[0] * val.shape[1], 1, 1) return {'val' : val} materials = { 'coef' : 'get_conductivity', } fields = { 'temperature' : ('real', 1, 'Omega', 1), } variables = { 'T' : ('unknown field', 'temperature', 0, True), 's' : ('test field', 'temperature', 'T'), } regions = { 'Omega' : 'all', 'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'), 'Gamma_Right' : ('vertices in (x > 0.099999)', 'facet'), } ebcs = { 'T1' : ('Gamma_Left', {'T.0' : 2.0}), 'T2' : ('Gamma_Right', {'T.0' : -2.0}), } def get_ic(coor, ic): """Non-constant initial condition.""" import numpy as nm # Normalize x coordinate. mi, ma = coor[:,0].min(), coor[:,0].max() nx = (coor[:,0] - mi) / (ma - mi) return nm.where( (nx > 0.25) & (nx < 0.75 ), 8.0 * (nx - 0.5), 0.0 ) functions = { 'get_ic' : (get_ic,), 'get_conductivity' : (get_conductivity,), } ics = { 'ic' : ('Omega', {'T.0' : 0.}), } integrals = { 'i' : 1, } equations = { 'Temperature' : """dw_volume_dot.i.Omega( s, dT/dt ) +dw_laplace.i.Omega( coef.val, s, T ) = 0""" } solvers = { 'ls' : ('ls.scipy_direct', {}), 'newton' : ('nls.newton', { 'i_max' : 1, 'eps_a' : 1e-10, 'eps_r' : 1.0, }), 'ts' : ('ts.simple', { 't0' : t0, 't1' : t1, 'dt' : None, 'n_step' : n_step, # has precedence over dt! 'quasistatic' : True, }), } options = { 'nls' : 'newton', 'ls' : 'ls', 'ts' : 'ts', 'save_steps' : -1, }